Shear state of freely evolving granular gases

Shear state of freely evolving granular gases J. Javier Brey, M. J. Ruiz-Montero, and A. Domínguez Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain 共Received 27 May 2008; published 1 October 2008兲 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. DOI: 10.1103/PhysRevE.78.041301 PACS number共s兲: 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 关1兴. 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 关2,3兴, among other reasons because they serve as a sensitive proving ground for kinetic theory and nonequilibrium 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 关4,5兴, 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. 关6兴. 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 关7,8兴. 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 关9,10兴. Consequently, the spatial inhomogeneities of the hydrodynamic fields are expected to be small. Although the flow is nonstationary, 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 关11–13兴. Nonstationary states in which all the time dependence occurs through the temperature, such as the homogenous cooling state 关14兴, as well as nonuniform states where all the spatial dependence also occurs through the temperature are peculiar of granular gases 关15兴. The existence of this free shear state has been suggested previously in the framework of a time-dependent Ginzburg- Landau model for granular gases 关10兴, and also as the result of a nonlinear analysis of the shearing instability 关16兴. 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 N-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 PHYSICAL REVIEW E 78, 041301 共2008兲 1539-3755/2008/78共4兲/041301共9兲©2008 The American Physical Society041301-1 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 n共r,t兲, the flow velocity u共r,t兲, and the granular temperature T共r,t兲of a system composed by smooth inelastic hard spheres 共d=3兲or disks 共d=2兲of mass m, diameter ␴ , and coefficient of normal restitution ␣ , have the form ⳵ n ⳵ t+⵱·共nu兲=0, 共1兲 冉 ⳵ ⳵ t+u·⵱ 冊 u+共mn兲−1 ⵱·P=0, 共2兲 冉 ⳵ ⳵ t+u·⵱+ ␨ 冊 T+2 nd共P:⵱u+⵱·q兲=0. 共3兲 For a dilute gas and to Navier-Stokes order, the pressure tensor P, the heat flux q, and the cooling rate ␨ are given by the constitutive relations 关9,17,18兴 Pij共r,t兲=p共r,t兲 ␦ ij − ␩ 冉 ⳵ uj ⳵ ri + ⳵ ui ⳵ rj −2 d ␦ ij ⵱·u 冊 ,共4兲 q共r,t兲=− ␬ ⵱T− ␮ ⵱n,共5兲 ␨ = ␨ 共0兲+ ␨ 共2兲.共6兲 Here p=nT 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 ␨ 共0兲denotes the contribution to the cooling rate of zeroth order in the gradients of the hydrodynamic fields, while ␨ 共2兲denotes the second order in the gradients contributions. The latter will be neglected in the following, as 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 关9兴. The transport coefficients appearing in Eqs. 共4兲and 共5兲can be expressed as ␩ 共T, ␣ 兲= ␩ 0共T兲 ␩ *共 ␣ 兲,共7兲 ␬ 共T, ␣ 兲= ␬ 0共T兲 ␬ *共 ␣ 兲,共8兲 ␮ 共n,T, ␣ 兲=T ␬ 0共T兲 n ␮ *共 ␣ 兲,共9兲 where ␩ 0共T兲and ␬ 0共T兲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 the Appendix. Finally, the expression of the zeroth order cooling rate reads ␨ 共0兲共n,T, ␣ 兲=nT ␩ 0共T兲 ␨ *共 ␣ 兲.共10兲 The function ␨ *共 ␣ 兲is also given in the Appendix. 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 nH, a vanishing velocity field uH=0, and a time dependent temperature TH共t兲obeying the equation ⳵ TH共t兲 ⳵ t+ ␨ H 共0兲共t兲TH共t兲=0. 共11兲 This is the so-called homogeneous cooling state 共HCS兲关14兴. 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 关4,5,9兴. This shear instability has been confirmed by molecular dynamics simulations 关4,5兴and also by Monte Carlo simulations of the effective dynamics associated to the Boltzmann equation 关19兴. 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 关20兴. 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 关20兴. 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 关6–8兴. Consider hydrodynamic flows verifying the following conditions. 共i兲There are gradients in only one direction taken as the Xaxis, and 共ii兲the hydrodynamic velocity field has the form ui共r,t兲= ␦ i,yu共x,t兲. 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, n=n共x兲, and the component Pxx of the pressure tensor must be uniform. In the framework of the Navier-Stokes approximation, the latter implies that the pressure is also uniform, p=p共t兲. Moreover, the equations for the velocity and temperature fields become BREY, RUIZ-MONTERO, AND DOMÍNGUEZ PHYSICAL REVIEW E 78, 041301 共2008兲 041301-2 ⳵ u ⳵ t−共mn兲−1 ⳵ ⳵ x 冉 ␩ ⳵ u ⳵ x 冊 =0, 共12兲 ⳵ T ⳵ t−2共nd兲−1 再 ␩ 冉 ⳵ u ⳵ x 冊 2 + ⳵ ⳵ x 冋 冉 ␬ −n ␮ T 冊 ⳵ T ⳵ x 册 冎 +T ␨ 共0兲=0. 共13兲 To solve the above equations, the boundary conditions must be specified. A system of Nparticles enclosed in a cubic 共d =3兲or square 共d=2兲box of side Lwill 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 ␸ 共x兲 by n共x兲=nH关1+ ␸ 共x兲兴,共14兲 with nH⬅N/Ld. The conservation of the number of particles and the periodic boundary conditions imply that 冕 0 L dx ␸ 共x兲=0, ␸ 共x+L兲= ␸ 共x兲,0艋x⬍L.共15兲 Moreover, it will be assumed that 兩 ␸ 共x兲兩Ⰶ1, 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 T共x,t兲⬇ ␪ 共t兲关1− ␸ 共x兲兴,共16兲 where ␪ 共t兲⬅p共t兲 nH =兰0 Ldxn共x兲T共x,t兲 兰0 Ldxn共x兲共17兲 is the spatial average temperature of the system at time t.It must be realized that assuming 兩 ␸ 兩Ⰶ1 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 ␪ 共t兲or the velocity field u共x,t兲. 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 l⬅ ␯ 共t兲 2 冋 m ␪ 共t兲 册 1/2 x共18兲 and s⬅1 2 冕 0 t dt⬘ ␯ 共t⬘兲,共19兲 respectively. Here ␯ 共t兲⬅nH ␪ 共t兲 ␩ 0共 ␪ 兲共20兲 is a characteristic frequency, proportional to the Boltzmann collision frequency ␯ B共t兲, namely, ␯ B共t兲=共2+d兲 ␯ 共t兲/4. Note that the prefactor 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 ⳵ u共l,s兲 ⳵ s− ␩ *共 ␣ 兲 2 ⳵ 2u共l,s兲 ⳵ l2=0, 共21兲 where terms of order u ␸ have been neglected. The solution of the above equation is a superposition of monocromatic waves of the form uk共l,s兲= ␹ k共l兲 ␻ k共s兲共22兲 with ␹ k共l兲= sin共kl + ␾ k兲共23兲 and ␻ k共s兲= ␻ k,0e−k2 ␩ *s/2,共24兲 where ␾ kand ␻ k,0 are arbitrary constants. The possible values of kare restricted by the periodic boundary conditions which imply k=2 ␲ q lM ,共25兲 qbeing a positive integer and lMthe value of lfor x=L, i.e., lM= ␯ 共t兲 2 冋 m ␪ 共t兲 册 1/2 L.共26兲 The form of Eq. 共24兲indicates that in the limit of large time s, the dynamics of the system is governed by the fundamental mode corresponding to the lowest possible value of k, i.e., k=km=2 ␲ /lM. Therefore, for large enough times, u共l,s兲⬇ukm共l,s兲= ␹ 共l兲 ␻ 共s兲共27兲 with ␹ 共l兲= ␹ km共l兲and ␻ 共s兲= ␻ km共s兲. In the same approximation under consideration, 兩 ␸ 共x兲兩Ⰶ1, the evolution equation for the temperature 共13兲leads to 1 ␪ 共s兲 d ␪ 共s兲 ds −m ␩ *共 ␣ 兲 ␻ 2共s兲 ␪ 共s兲d 冉 d ␹ 共l兲 dl 冊 2 + ␨ *共 ␣ 兲关2+ ␸ 共l兲兴 +d+2 2共d−1兲关 ␬ *共 ␣ 兲− ␮ *共 ␣ 兲兴d2 ␸ 共l兲 dl2=0. 共28兲 Using Eq. 共23兲this becomes SHEAR STATE OF FREELY EVOLVING GRANULAR GASES PHYSICAL REVIEW E 78, 041301 共2008兲 041301-3 1 ␪ 共s兲 d ␪ 共s兲 ds − ␩ *共 ␣ 兲m ␻ 2共s兲km 2 2 ␪ 共s兲d兵1 + cos关2共kml+ ␾ km兲兴其 + ␨ *共 ␣ 兲 ⫻关2+ ␸ 共l兲兴 +d+2 2共d−1兲关 ␬ *共 ␣ 兲− ␮ *共 ␣ 兲兴d2 ␸ 共l兲 dl2=0. 共29兲 By requiring the sum of the position dependent terms to cancel it is obtained that ␸ 共l兲=− 共d−1兲 ␩ *共 ␣ 兲 4d共d+2兲 ␥ 共km, ␣ 兲acos关2共kml+ ␾ km兲兴,共30兲 where ␥ 共 ␣ ,km兲⬅ ␬ *共 ␣ 兲− ␮ *共 ␣ 兲−共d−1兲 ␨ *共 ␣ 兲 2共d+2兲km 2共31兲 and a⬅m ␻ 2共s兲 ␪ 共s兲.共32兲 Consistency requires that abe actually independent of sand, moreover, that it be defined positive. The former of these conditions together with Eq. 共21兲leads to 1 ␪ 共s兲 d ␪ 共s兲 ds =2 ␻ 共s兲 d ␻ 共s兲 ds =−km 2 ␩ *共 ␣ 兲.共33兲 This equation will be used later on to determine the effective cooling rate of the average temperature ␪ 共s兲of the system in the free shear state. On the other hand, equating to zero the part of Eq. 共29兲that is position independent yields, after employing Eq. 共33兲, a=2d关2 ␨ *共 ␣ 兲−km 2 ␩ *共 ␣ 兲兴 km 2 ␩ *共 ␣ 兲.共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 km 2⬍2 ␨ * ␩ *共35兲 or, equivalently, L⬎Lc, with the critical size Lcgiven by Lc=共2+d兲⌫共d/2兲 2 ␲ 共d−3兲/2nH ␴ 共d−1兲 冉 ␩ * 2 ␨ * 冊 1/2 .共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 关9兴. 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 ␸ 共l兲=−Acos关2共kml+ ␾ km兲兴,共37兲 A=共d−1兲 ␩ *共 ␣ 兲 2共d+2兲 ␥ 共km, ␣ 兲 冋 冉 L Lc 冊 2 −1 册 .共38兲 Therefore, the condition 兩 ␸ 共x兲兩Ⰶ1 formally implies that 兩L −Lc兩/LcⰆ1. 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 ␻ ˜ 共s兲⬅ ␻ 共s兲 冋 m 2 ␪ 共s兲 册 1/2 =d1/2 冋 冉 L Lc 冊 2 −1 册 1/2 .共39兲 Therefore, the scaled macroscopic velocity field is timeindependent and all its dependence on the average density and inelasticity occurs through the critical length Lc. A similar expression can be obtained for the ratio between the effective cooling rate of the shear state ␨ S *and the cooling rate of the HCS, ␨ *. The former is identified by writing Eq. 共33兲 in the form ⳵ ␪ 共s兲 ⳵ s+2 ␨ S * ␪ 共s兲=0, 共40兲 with ␨ S *=km 2 ␩ *共 ␣ 兲 2,共41兲 that leads to ␨ S * ␨ *= 冉 Lc L 冊 2 .共42兲 Of course, in the limit L→Lc, ␪ 共s兲→TH共s兲and the law for the HCS given by Eq. 共11兲is recovered. The average energy per particle is e ¯ T⬅m 2u2+T ¯ d 2,共43兲 with the bar denoting spatial average. Using Eqs. 共17兲,共22兲, 共23兲,共32兲, and 共34兲it is found that e ¯ T共s兲= ␨ *共 ␣ 兲 ␪ 共s兲d km 2 ␩ *共 ␣ 兲= 冉 L Lc 冊 2 ␪ 共s兲d 2.共44兲 For L→Lc, the equilibrium relation e ¯ T=dT/2 is recovered as expected. This equation has been previously obtained by Wakou et al. 关10兴, 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 共1− ␣ Ⰶ1兲. 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. 关16兴have carried out a nonlinear analysis of the hydrodynamic equations of a system of inelastic hard disks close to the instability threshold. When kmis slightly BREY, RUIZ-MONTERO, AND DOMÍNGUEZ PHYSICAL REVIEW E 78, 041301 共2008兲 041301-4 smaller than kc⬅2 ␲ /Lc, vorticity modes with wave number kmexhibit 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 k=kmby 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 Ref. 关16兴do not give much 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 关16兴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 关21,22兴, 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 Lwas always larger than, but close to, the critical one Lc, predicted by Eq. 共36兲, for the considered value of the restitution coefficient ␣ . Then, according to the theory, the system is expected to generate the nonlinear free shear flow described in Sec. III, if it is stable. As usual in DSMC simulations, the mass mof the particles will be taken as the unit of mass, the average mean free path ␭⬅共2冑2nH ␴ 兲−1 as the unit of length, and 2T共0兲 ⬅2 ␪ 共0兲, where T共0兲is the initial temperature, as the unit of energy. The effect of a collision between particles rand sis to instantaneously modify their velocities according to the rule vr→vr ⬘=vr−1+ ␣ 2共 ␴ ˆ·vrs兲 ␴ ˆ,共45兲 vs→vs ⬘=vs+1+ ␣ 2共 ␴ ˆ·vrs兲 ␴ ˆ,共46兲 where vrs⬅vr−vsis the relative velocity and ␴ ˆis the unit vector pointing from the center of particle sto the center of particle rat contact. This corresponds to the scenario in which the constitutive relations 共4兲–共6兲were derived 关9,17兴. 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 关23,24兴. 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. A. 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 0 10 20 30 10 20 30 τ=241.2 (b) 0 10 20 30 10 20 30 τ=120.5 (a) 0 10 20 30 10 20 30 τ=363.2 (c) 0 10 20 30 10 20 30 τ=604.5 ( d ) FIG. 1. Time evolution of the velocity field scaled with the average thermal velocity u共m/2 ␪ 兲1/2in a system with ␣ =0.95 and L/Lc⬇1.205. The indicated times ␶ are measured as the average number of collisions per particle. The lengths of the arrows in the plots are proportional to the scaled velocities. SHEAR STATE OF FREELY EVOLVING GRANULAR GASES PHYSICAL REVIEW E 78, 041301 共2008兲 041301-5 the scaled velocity field u共m/2 ␪ 兲1/2is plotted at four different times for a system with ␣ =0.95 and L=39 共L/Lc ⬇1.205兲. 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 xaxis 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. 2the ratio between the spatial average temperature ␪ 共t兲and the HCS temperature at the same time as predicted by the Haff law 共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 共 ␶ ⱗ250兲in spite of the system having already well developed vortices 共see Fig. 1兲. Afterwards, for 250ⱗ ␶ ⱗ600, 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 ␪ 共t兲still obeys a Hafflike 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 ␣ =0.95 and L/Lc⬇1.205. Two Fourier components of the relative density field ␳ ⬅n/nHare plotted ␳ 2kmand ␳ km,km. 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, ␳ km,kmis 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. B. 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 Anis plotted in Fig. 4as a function of L/Lcfor three values of the restitution coefficient ␣ =0.97, 0.95, and 0.9. 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 ␣ .Itis seen that the agreement is quite good when L/Lcis close to 1, 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 AnⰆ1. In Fig. 5, the amplitude ATof the steady cosine profile for 1−T共x,t兲/ ␪ 共t兲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 0 200 400 600 800 100 0 τ 0 1 2 3 θ/ THCS FIG. 2. 共Color online兲Time evolution of the spatial average temperature normalized by the temperature of the reference HCS for the same system as in Fig. 1. Time ␶ is measured as the average number of collisions per particle. The several curves correspond to different simulation trajectories. 0 200 400 600 800 100 0 τ 0 0.05 0.1 ρ FIG. 3. 共Color online兲Time evolution of the second transversal mode of the relative density ␳ 2km共solid lines兲and the diagonal mode ␳ km,km共dashed lines兲. The system is the same as in Fig. 1. The several curves correspond to different trajectories, and time ␶ is measured again as the average number of collisions per particle. BREY, RUIZ-MONTERO, AND DOMÍNGUEZ PHYSICAL REVIEW E 78, 041301 共2008兲 041301-6 the density profile, i.e., AT=An=A. Nevertheless, the simulation data indicate that this is not the case, except for values of L/Lcclose to unity. The deviations of the numerical data from the theoretical predictions are in opposite directions for ATand An, 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 Pxx of the pressure tensor was found to be uniform, as required by the balance equations 共1兲–共3兲. Nevertheless, for L/Lcsmall enough, the agreement between theory and simulation can be qualified as satisfactory. Let us mention that if the elastic values of the transport coefficients 共 ␩ *= ␬ * =1, ␮ *=0兲were used, the theoretical prediction for Aas a function of L/Lcwould 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 L/Lc. This latter property was seen to be verified in the simulation within the statistical uncertainties. Moreover, the results displayed in Fig. 6show 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 共39兲. In Fig. 7, the ratio between the cooling rate of the free shear state, ␨ S *, 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 tempera- 1 1.2 1.4 1.6 L/Lc 0 0.1 0.2 0.3 A n α=0.97 α=0.95 α=0.90 FIG. 4. 共Color online兲Dimensionless amplitude of the cosine density profile in the free shear state as a function of the ratio between the size of the system Land its critical value for the shear instability Lc. Three values of the restitution coefficient are considered, as indicated. The symbols are simulation results and the lines are the predictions of Eq. 共38兲, the lowest line corresponding to the larger ␣ . 1 1.2 1.4 1.6 L/L c 0 0.1 0.2 0.3 0 . 4 A T α=0.97 α=0.95 α=0.90 FIG. 5. 共Color online兲The same as in Fig. 4but for the amplitude of the cosine temperature profile. 1 1.2 1.4 1.6 L/L c 0 0.5 1 1.5 2 2 .5 ~ ω α=0.97 α=0.95 α=0.90 FIG. 6. 共Color online兲Dimensionless amplitude of the macroscopic flow ␻ ˜ ⬅ ␻ 共t兲/关2 ␪ 共t兲/m兴1/2in the free shear state as a function of the length Lof the system normalized by its critical value Lc. The solid line is the theoretical prediction 共39兲, while the symbols are simulation results for the same systems as in Fig. 4. 1 1.2 1.4 1.6 L/L c 0.2 0.4 0.6 0.8 1 ζ∗ S/ζ∗α=0.97 α=0.95 α=0.90 FIG. 7. 共Color online兲Ratio between the cooling rates of the free shear state and of the HCS as a function of the ratio L/Lc. The symbols have the same meaning as in Fig. 4and the solid line is the theoretical prediction 共41兲. SHEAR STATE OF FREELY EVOLVING GRANULAR GASES PHYSICAL REVIEW E 78, 041301 共2008兲 041301-7 ture. 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 nonstationary 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 nonequilibrium 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 Lc.Asa 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 is the case for the steady uniform shear state of a granular gas 关25兴. 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 Anand AT. In Ref. 关6兴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 twovortex 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 Lbecomes much larger than Lc共typically, L/Lcⲏ1.7兲, 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. 关16兴. 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 关26兴. Moreover, in Ref. 关16兴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 u, and when the expressions are inserted into the equation for u, 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 ␪ 共t兲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 关27兴and/or velocity dependent restitution coefficients 关2兴, 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. ACKNOWLEDGMENT This research was supported by the Ministerio de Educación y Ciencia 共Spain兲through Grant No. FIS2008-01339 共partially financed by FEDER funds兲. APPENDIX: 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 ␩ 0=2+d 8⌫ 冉 d 2 冊 ␲ −共d−1兲/2共mT兲1/2 ␴ −共d−1兲,共A1兲 ␬ 0=d共d+2兲2 16共d−1兲⌫ 冉 d 2 冊 ␲ −共d−1兲/2 冉 T m 冊 1/2 ␴ −共d−1兲,共A2兲 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 ␩ *共 ␣ 兲= 冋 ␯ 1 *共 ␣ 兲− ␨ *共 ␣ 兲 2 册 −1 ,共A3兲 ␬ *共 ␣ 兲= 冋 ␯ 2 *共 ␣ 兲−2d d−1 ␨ *共 ␣ 兲 册 −1 关1+c*共 ␣ 兲兴,共A4兲 ␮ *共 ␣ 兲=2 ␨ *共 ␣ 兲 冋 ␬ *共 ␣ 兲+共d−1兲c*共 ␣ 兲 2d ␨ *共 ␣ 兲 册 ⫻ 冋 2共d−1兲 d ␯ 2 *共 ␣ 兲−3 ␨ *共 ␣ 兲 册 −1 ,共A5兲 BREY, RUIZ-MONTERO, AND DOMÍNGUEZ PHYSICAL REVIEW E 78, 041301 共2008兲 041301-8 ␨ *共 ␣ 兲=2+d 4d共1− ␣ 2兲 冋 1+3c*共 ␣ 兲 32 册 .共A6兲 In the above expressions, ␯ 1 *共 ␣ 兲=共3−3 ␣ +2d兲共1+ ␣ 兲 4d 冋 1−c*共 ␣ 兲 64 册 ,共A7兲 ␯ 2 *共 ␣ 兲=1+ ␣ d−1 冋 d−1 2+3共d+8兲共1− ␣ 兲 16 +4+5d−3共4−d兲 ␣ 1024 c*共 ␣ 兲 册 ,共A8兲 c*=32共1− ␣ 兲共1−2 ␣ 2兲 9+24d+共8d−41兲 ␣ +30 ␣ 2共1− ␣ 兲.共A9兲 It is easily verified that ␬ *and ␩ *tend to unity when ␣ goes to one, while ␮ *and ␨ *vanish in this limit. 关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. Comput. 8,1 共1993兲. 关5兴S. McNamara and W. R. Young, Phys. Rev. E 50,R28共1994兲; 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. M. 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兲. SHEAR STATE OF FREELY EVOLVING GRANULAR GASES PHYSICAL REVIEW E 78, 041301 共2008兲 041301-9