Attempted density blowup in a freely cooling dilute granular gas: hydrodynamics versus molecular dynamics
Andrea Puglisi, Michael Assaf, Itzhak Fouxon, Baruch Meerson
aa r X i v : . [ c ond - m a t . s o f t ] J a n Attempted density blowup in a freely cooling dilute granular gas: hydrodynamicsversus molecular dynamics
Andrea Puglisi , Michael Assaf , Itzhak Fouxon , and Baruch Meerson Dipartimento di Fisica, Universit`a La Sapienza, p.le Aldo Moro 2, 00185 Roma, Italy Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel
It has been recently shown (Fouxon et al.
PACS numbers: 45.70.Qj
I. INTRODUCTION
Spontaneous clustering of particles in granular gaseshas attracted much recent interest [1, 2, 3, 4, 5, 6, 7, 8,9, 10, 11, 12]. As other pattern-forming instabilities, theclustering instability of a freely cooling granular gas hasserved as a sensitive probe of theoretical modeling and,first of all, of the Navier-Stokes granular hydrodynam-ics (GHD). Although the formal criteria of its validitymay be quite restrictive, see below, GHD has a greatpower, sometimes far beyond its formal validity limits[13], in predicting a host of collective phenomena in gran-ular flows, such as shocks, vortices and clusters. In therecent years, GHD has been applied to a variety of non-stationary dilute granular flows [10, 14, 15, 16, 17, 18].Non-stationary flows are both appealing and challeng-ing for continuum modeling of granular dynamics. As inother areas of continuum modeling, this is especially truewhen a non-stationary flow develops a finite-time singu-laritiy [19]. Examples are provided by the finite-timeblowup of the gas density: at zero gravity [10, 16, 17](as described by the ideal
GHD discussed below, IGHDfrom now on), and at finite gravity [15] as described bythe more complete, non-ideal GHD. Of course, a den-sity blowup in a gas with finite-size particles can only bean intermediate asymptotics, as the blowup is ultimatelyarrested: either by close-packing effects [11], or by thegradient-dependent transport [18]. Still, the attemptedblowup regimes, signaling the development of high den-sity regions in the gas, are fascinating and worth a de-tailed investigation. One such regime has been recentlyaddressed by Fouxon et al. [16, 17]. They dealt witha macroscopically one-dimensional, freely cooling, dilutegranular flow in the framework of ideal GHD that ne-glects the gradient-dependent transport effects: the heat diffusion and viscosity. Fouxon et al. derived a class ofexact solutions to the ideal equations, for which an ini-tially smooth flow develops a finite-time density blowup.Close to the blowup time τ , the maximum gas densityexhibits a power law behavior ∼ ( τ − t ) − . The velocitygradient blows up as ∼ − ( τ − t ) − , whereas the veloc-ity itself remains continuous and forms a cusp, ratherthan a shock discontinuity, at the singularity. The gastemperature vanishes at the singularity, but the pressureremains finite and almost constant. Extensive numeri-cal simulations with the ideal hydrodynamic equationsshowed that the singularity, exhibited by the exact solu-tions, is universal, as it develops for quite general initialconditions [16, 17]. The reason behind this universalityis in the fact that the sound travel time through the re-gion of the developing singularity is much shorter thanthe characteristic inelastic cooling time of the gas in thatregion. As a result, the pressure gradient (almost) van-ishes in the singularity region, and the local features ofthis isobaric singularity become essentially independentof how the flow was initiated and how it behaves at largedistances from the singularity. This singularity is of thesame type as the one that develops, in the framework ofthe IGHD equations, in a general low Mach number flowof a freely cooling granular gas [18].Here we perform molecular dynamics (MD) simula-tions of a freely cooling gas of nearly elastically collidinghard disks, aimed at identifying the “attempted” densityblowup regime, predicted by the ideal analytical solu-tions [16, 17]. We simulate a freely evolving dilute gas ofnearly elastically colliding hard disks in a narrow channelwith perfectly elastic side walls. In this geometry boththe clustering mode in the transverse directions, and theshear mode are suppressed (see Refs. [2, 3, 10, 11] fordetailed criteria). As a result, the coarse-grained, or hy-drodynamic flow depends only on the longitudinal coor-dinate along the channel (and time), as it was assumedin Refs. [16, 17]. We choose the initial conditions of theMD simulations so that the coarse-grained density, ve-locity and temperature fields are those producing one ofthe exact blowup solutions of the ideal GHD equations.Then we follow the time history of the hydrodynamicfields in the MD simulations and compare it with thatpredicted by the ideal exact solution and with numericalsolutions of non-ideal GHD equations.The remainder of this paper is organized as follows.In Section II we briefly summarize the ideal GHD anal-ysis [16, 17] of the density blowup: we present the idealGHD equations, one of their exact solutions, its mainfeatures and expected limits of its validity. In SectionIII we describe our MD simulations and compare the hy-drodynamic quantities, computed from the simulations,with the exact solution of the ideal GHD equations. Wefind that the exact solution is in remarkable quantitativeagreement with the MD simulations over an extendedtime interval, proving the existence of the attempted den-sity blowup regime. As the attempted singularity is ap-proached, the hydrodynamic fields, as observed in theMD simulations, deviate from the predictions of the ex-act solution. To investigate the mechanism of breakdownof the ideal solution, we extend the hydrodynamic the-ory, in Section IV, in two separate ways. In the first onewe take into account the gradient-dependent transport:the heat diffusion and viscosity, but continue to assumethat the gas is dilute. In the second one we neglect thegradient-dependent transport but take into account, ina semi-phenomenological way, finite density corrections.Section V summarizes our results and puts them into amore general context of hydrodynamic scenarios of clus-tering in freely evolving granular gases. II. HYDRODYNAMIC THEORY AND DENSITYBLOWUPA. Ideal granular hydrodynamics and exactsolution
We consider a two-dimensional granular gas of identi-cal hard and smooth disks with diameter σ and mass setto unity, and adopt a simple model where the inelasticparticle collisions are characterized by a constant coeffi-cient of normal restitution r . Throughout this paper wewill only deal with nearly elastic collisions,1 − r ≪ , (1)and assume a very small Knudsen number: l free /L ≪ . (2)Here l free is the mean free path of the particles, and L isthe characteristic length scale of the hydrodynamic fieldsthat may depend on time. In addition, we will assume in this Section that the local gas density ρ is much smallerthan the close-packing density of disks ρ c = 2 / ( √ σ ): ρσ ≪ . (3)The strong inequalities (2) and (3) need to be verified aposteriori , once the hydrodynamic problem in questionis solved. The strong inequalities (1)-(3) enable one toemploy the well-established equations of Navier-Stokesgranular hydrodynamics (see, e.g. [13, 20]) that dealwith three coarse-grained fields: the mass density ρ ( x , t ),the mean flow velocity v ( x , t ) and the granular tempera-ture T ( x , t ). In a sufficiently narrow channel these fieldsdepend only on the longitudinal coordinate x , and thehydrodynamic equations become ∂ρ∂t + ∂ ( ρv ) ∂x = 0 , (4) ρ (cid:18) ∂v∂t + v ∂v∂x (cid:19) = − ∂ ( ρT ) ∂x + ν ∂∂x (cid:18) √ T ∂v∂x (cid:19) , (5) ∂T∂t + v ∂T∂x = − T ∂v∂x − Λ ρT / + κ ρ ∂∂x (cid:18) √ T ∂T∂x (cid:19) + ν √ Tρ (cid:18) ∂v∂x (cid:19) , (6)where Λ = √ π (1 − r ) σ , ν = (2 √ πσ ) − and κ =2 / ( √ πσ ), see e.g. Refs. [6, 20]. Equations (4)-(6) differfrom the hydrodynamic equations for a gas of elastically colliding disks only by the presence in Eq. (6) of the in-elastic loss rate term − Λ ρT / , that describes the pro-portionality of the energy loss per particle to the num-ber of particle collisions per unit time (proportional to ρT / ) and to the energy loss per collision (proportionalto T ). This additional term, however, brings a whole newphysics (and mathematics) into the problem.Let us rewrite Eqs. (4)-(6) in dimensionless variables.We will measure the gas density, temperature and thevelocity in the units of ρ , T / p T /
2, respectively,where ρ and T are some characteristic values of theinitial density and temperature. The time and distancewill be measured in the units of τ = 4Λ ρ √ T and l = τ r T , (7)respectively. As one can see from Eq. (6), τ is the char-acteristic cooling time of the gas due to the collisionalenergy loss, while l is the characteristic distance a soundwave travels during time τ . The numerical factors inEqs. (7) are chosen for convenience. We will keep identi-cal notation for the rescaled and physical quantities, andtake care that no confusion arises. Using the rescaledquantities, we rewrite Eqs. (4)-(6) as ∂ρ∂t + ∂ ( ρv ) ∂x = 0 , (8) ρ ∂v∂t + ρv ∂v∂x = − ∂ ( ρT ) ∂x + 1 − r √ ∂∂x (cid:18) √ T ∂v∂x (cid:19) , (9) ∂T∂t + v ∂T∂x = − T ∂v∂x − √ ρT / + 1 − r √ ρ ∂∂x (cid:18) √ T ∂T∂x (cid:19) + (1 − r ) √ T √ ρ (cid:18) ∂v∂x (cid:19) , (10)Let us assume that the characteristic magnitudes of the rescaled hydrodynamic fields, and of their spatial andtemporal derivatives, are of order unity (this assumptionneeds to be checked a posteriori ). Then we can neglectthe viscous and thermal conduction terms, as they scaleas 1 − r ≪
1, and arrive at the IGHD equations: ∂ρ∂t + ∂ ( ρv ) ∂x = 0 , ρ (cid:18) ∂v∂t + v ∂v∂x (cid:19) = − ∂ ( ρT ) ∂x , (11) ∂T∂t + v ∂T∂x = − T ∂v∂x − √ ρT / . (12)These equations were investigated in Refs.[16, 17], wherea family of exact analytic solutions was derived. Herewe will consider a representative and simple particularsolution that evolves from the following initial conditions: ρ ( x, t = 0) = cosh − x, T ( x, t = 0) = 2 , (13) v ( x, t = 0) = − x ) . (14)To remind the reader, we are using rescaled variableshere. Back in the physical variables, the initial profilesare ρ ( x,
0) = ρ cosh( x/l ) , (15) T ( x,
0) = T , (16) v ( x,
0) = − p T arcsin h tanh (cid:16) xl (cid:17)i , (17)where l is defined in Eq. (7). That is, at t = 0 the densityprofile has a maximum ρ and width l , the temperature T is uniform, and there is an inflow of the gas towardsthe origin with v ( x → ±∞ , t = 0) = ∓ π p T /
2. Theinitial scale of variation of the fields, l , is by a factor1 / (1 − r ) greater than the mean free path of the gas,justifying the use of the hydrodynamic description. Fur-thermore, both the magnitudes, and the spatial scales ofthe rescaled fields are of order unity which justifies, atleast for finite times, the use of the ideal equations (11)and (12). Figure 1 shows the initial density and velocityfields of the flow in the Eulerian coordinate.Now we go over from the Eulerian coordinate x to theLagrangian mass coordinate m = R x ρ ( x ′ , dx ′ , see e.g. [21]. For the initial density profile (13), the Euleriancoordinate x is related to the Lagrangian coordinate m as follows: x ( m, t = 0) = 12 ln (cid:18) m − sin m (cid:19) . (18) ρ / ρ x/l -4-3-2-10 v / ( T / ) / FIG. 1: The initial values of the hydrodynamic fields in theexample of attempted density blowup considered in this work.Shown are the rescaled initial density and velocity of the gas[see Eqs. (15) and (17)] versus the rescaled Eulerian coordi-nate x/l along the channel. Only the right half of the systemis shown. The values of ρ , T and l , used in our molecu-lar dynamic simulations, are presented in the beginning ofsubsection III.D. Note that the infinite Eulerian interval −∞ < x < + ∞ corresponds, in this example, to a finite interval of m : − π/ < m < π/
2, that is to a finite (rescaled) total massof the gas, equal to π . In the Lagrangian coordinates,the ideal equations (11) and (12) are ∂∂t (cid:18) ρ (cid:19) = ∂v∂m , ∂v∂t = − ∂p∂m , (19) ∂p∂t = − pρ ∂v∂m − √ p / ρ / , (20)where the dilute gas pressure p = ρT has been used in-stead of the temperature. The exact solution in the La-grangian coordinates is the following [16, 17]: ρ ( m, t ) = cos m (1 − t cos m ) , p ( m, t ) = 2 cos m, (21) v ( m, t ) = − m + 2 t sin m. (22)As x ( m, t ) can be calculated explicitly [17], x ( m, t ) = 12 ln (cid:18) m − sin m (cid:19) − tm + t sin m, (23)equations (21)-(23) describe the time history of the hy-drodynamic fields in the x -coordinate in a closed para-metric form. Let us consider some important features ofthis simple exact solution. B. Density blowup and its properties
A momentary look at the first of Eqs. (21) reveals thatthe density at the origin blows up in a finite time. Backin the dimensional variables one obtains ρ ( x = 0 , t ) = ρ (1 − t/τ ) . (24)The gas temperature becomes zero at the singularity, thevelocity gradient ∂v ( x, t ) /∂x blows up as ∼ − ( τ − t ) − ,whereas the velocity itself remains continuous and formsa cusp. This behavior at singularity is quite differentfrom that of the free flow singularity (that one observesfor the same initial velocity profile but a zero gas pres-sure). In particular, for the free flow singularity the den-sity blows up as (1 − t/τ ) − , while the velocity developsa shock discontinuity [22]. See Refs. [16, 17] for moredifferences between the two types of the singularities.It was found in Refs. [16, 17] that the finite timeblowup of the gas density and of the velocity gradient isnot a consequence of specially chosen initial conditions.Rather, it appears, in the framework of the ideal equa-tions (11) and (12) [or, equivalently, of Eqs. (19) and(20)], for quite general initial conditions. A robust lo-cal feature of this singularity is a finite, non-zero valueof the gas pressure at the time of the density blowup.The singularity is universal because the sound travel timethrough it is much shorter than the characteristic inelas-tic cooling time of the gas. As a result, the pressuregradient (almost) vanishes in the singularity region, andthe local features of this isobaric singularity become in-dependent of the flow details at large distances.As an infinite density cannot be reached in a gas withfinite-size particles, it is clear that, sufficiently close tothe attempted singularity, some of the assumptions madeon the way to the ideal equations (11)-(12) break down.However, independently of the precise way of regulariz-ing the infinite density, an attempted singularity impliesthe formation of a region with a very high density. Fur-thermore, the ideal solution should accurately describethe evolution of the system for times not too close to theattempted density blowup time τ . Before we look intowhere the ideal theory breaks down, let us consider some global characteristics of the exact solution. For these thesingularity time t = τ turns out not to be special. First,we note that the solution describes a gas with a (con-stant) finite number N of particles, given by N = L y Z ∞−∞ ρ ( x ) dx = πρ L y l = 2 √ π L y (1 − r ) σ , (25)where L y is the channel width. Note that N is indepen-dent of ρ : for larger ρ the initial density profile hasa higher peak but a smaller width, so that N remainsconstant. Another global characteristics of the solutionis the total energy of the gas E ( t ) = L y Z ∞−∞ (cid:18) ρv ρT (cid:19) dx . For the thermal part of the energy we find, after a simple algebra, E th ( t ) L y = Z ∞−∞ ρT dx = 2 Z ∞ ρT dx = 2 ρ T l Z π/ [1 − ( t/τ ) cos m ] dm = ρ T l " π − (cid:18) tτ (cid:19) + π (cid:18) tτ (cid:19) . (26)Similarly, for the kinetic energy of the macroscopic mo-tion we find E kin ( t ) L y = ρ T l " π − (cid:18) tτ (cid:19) + π (cid:18) tτ (cid:19) . (27)Summing the two and using Eq. (25), we obtain [17] E ( t ) = N T " π − π (cid:18) tτ (cid:19) + (cid:18) tτ (cid:19) . (28)As time increases, E ( t ) is monotone decreasing for t ≤ τ . We also observe that t = τ is a regular point of E ( t ), where nothing dramatic happens. Note that theparabolic-in-time law of the energy decay is quite differ-ent from Haff’s law E ( t ) = E (0) / (1 + 2 t/τ ) obtained forfreely cooling homogeneous granular gas with density ρ [23]. C. Breakdown of ideal theory
For a one-dimensional flow, the applicability of the so-lution (21)-(23) is determined by the applicability of theideal equations (11)-(12) that it solves exactly. Analysisin Ref. [17] shows that, sufficiently close to the attemptedsingularity, the ideal equations become invalid. This hap-pens because of one of two reasons (or both): either thegas becomes dense, so that criterion (3) breaks down, orthe heat diffusion becomes important, invalidating theideal equations (11)-(12). The time t br at which the idealtheory breaks down can be estimated as follows [17]:1 − t br τ ∼ max( p ρ σ , − r ) . (29)Therefore, the “bottleneck” for the validity of the equa-tions is set by the initial conditions: if the maximumis determined by the first (correspondingly, the second)term in the right hand side of Eq. (29), the ideal equa-tions become invalid because of the finite gas density(correspondingly, the finite heat diffusion). As each ofthe two terms is very small by assumption, the solutionis expected to break down only close to the attemptedsingularity [24].The one-dimensional solution may also become invalidbecause of instability with respect to small initial pertur-bations that are inevitably present in MD simulations.Numerical solutions of the hydrodynamic equations, re-ported in Refs. [16, 17], strongly suggest that the idealsolution is stable with respect to small longitudinal per-turbations. This does not exclude possible instabilitywith respect to small transverse perturbations. The onlyavailable analytic result here is the one obtained fromthe stability condition for a homogeneous cooling state,see Refs. [2, 3, 10, 11]. That stability criterion comesfrom a competition between the (destabilizing) inelasticcooling and the (stabilizing) heat diffusion and viscosityin the transverse direction. The stability criterion de-mands that L y be less than a threshold value dependingon 1 − r , ρ and σ . The stability problem for the stronglyinhomogeneous and time-dependent exact solution is ob-viously more complicated, and its complete analytic so-lution does not seem feasible. It is therefore importantthat our MD simulations, presented in the next Section,strongly suggest that, for sufficiently narrow channels,no instability in the transverse direction occurs for thetime-dependent flow we are working with.Let us summarize the main theoretical predictions. Forthe initial conditions (13) and (14), a nonlinear time-dependent flow sets in, described by the ideal exact solu-tion: Eqs. (21)-(23). This flow “attempts” to develop adensity blowup. However, close to the attempted singu-larity one (or both) of the two factors: the finite densityand the heat diffusion, invalidates the solution. The rel-ative importance of the two factors is determined, viaEq. (29), by the initial conditions. III. MOLECULAR DYNAMICS SIMULATIONSA. General
To test the theoretical predictions, we performed MDsimulations of a free cooling granular gas in a narrow two-dimensional channel. The initial conditions correspondto hydrodynamic profiles (13) and (14) and satisfy thestrong inequalities (1)-(3). According to the theory, theyare expected to generate the nonlinear time-dependentflow described by Eqs. (21)-(23). Our MD simulationcalculates the evolution of a gas of N ′ identical inelastichard disks of unit mass, with diameter σ , in a channel ofdimensions L ′ x = L x / L y . As the expected hydro-dynamic flow is symmetric with respect to x = 0, onlyone half of the system, x ∈ [0 , L x / N ′ = N/
2. Each wall of (the one half of) the channelis solid and reflects elastically the disks colliding with it.The particles move freely until a collision (“event”) oc-curs when two disks i and j find themselves at a distanceequal to σ . The collision is resolved instantaneously, leav-ing the positions of the particles unaltered and updat-ing their velocities from ( v i , v j ), before the collision, to( v ′ i , v ′ j ), after the collision. The update rule conservesthe total momentum and reduces the total kinetic energy, with a constant coefficient of normal restitution r ∈ [0 , v ′ i = v i − r g · ˆ σ ) ˆ σ , (30) v ′ j = v j + 1 + r g · ˆ σ ) ˆ σ , (31)where g = v i − v j and ˆ σ is the unit vector joining thecenters of the two disks. The hard-core interactions makepossible the following optimization of the algorithm. Itis sufficient to calculate the first collision times of allparticles and then select the absolute first one. Thesystem is freely evolved up to that time, then the col-lision is resolved, and a new list of collision times iscomputed. With standard optimization techniques of thesearch procedure it is possible to achieve fast computa-tion times [25]. Nevertheless, the time performance isproportional to the number of collisions occurred, so theratio between the physical time and the cpu-time goesdown when dense clusters emerge in the system. B. Initial conditions
The initial position and velocity of each of the N ′ disksare chosen randomly with probability distributions corre-sponding to the desired initial hydrodynamic fields. Thiswas implemented with the following procedure. For eachdisk i
1. the longitudinal position x i is chosen with probabil-ity proportional to ρ ( x,
0) from Eq. (15) for x ≥ x i -position is generated with uni-form probability on the interval [ σ/ , L ′ x − σ/ z with uniform probabilityon [0 , max { ρ ( x, } ] is compared with ρ ( x i , z < ρ ( x i ,
0) the position is accepted, other-wise the procedure is repeated from 1a,2. then the vertical position y i is chosen with uniformprobability on the interval [ σ/ , L y − σ/ x i , y i ) and all the previouslyplaced disk centers must be greater then σ : if thecondition is not satisfied, the procedure is repeatedfrom 1a,4. the velocity components v xi and v yi are chosen froma Gaussian distribution with zero mean and vari-ance equal to T ; then the longitudinal componentis shifted by an amount v ( x,
0) from Eq. (17) with x ≥ C. Lagrangian coordinate and hydrodynamic fields
We verified that, for our choice of the channel dimen-sions, the gas remained homogeneous in the y -direction.By virtue of this observation, it was sufficient to deal withone-dimensional hydrodynamic profiles depending on x .For a direct comparison with the analytical solution ofthe IGHD, the hydrodynamic profiles were obtained us-ing a uniform binning in the Lagrangian mass coordinate.Using the same notation as in Section I, we define the La-grangian mass interval for the simulated flow as [0 , π/ π/ N ′ = N/
2. Let n bin be the number of bins chosen to sam-ple the hydrodynamic profiles and N bin = N ′ /n bin be theaverage number of particles per bin. At a given time t all particles are ordered so that x i < x i +1 , i ∈ [0 , N ′ − j ∈ [1 , n bin ] has its leftmost border at x [( j − N bin ] and its rightmost border at x [ jN bin − . Thesebins are non-uniform in the x coordinate, but are uni-form in the Lagrangian mass coordinate as each containsthe same mass N bin . The position of the j -th bin is m j = π ( j − / n bin . All particles belonging to the j -th bin contribute to thevalue of the hydrodynamic fields: ρ ( m j , t ) = N bin L y L j , (32) v ( m j , t ) = P i ∈ j v xi N bin , (33) T ( m j , t ) = P i ∈ j (cid:2) ( v xi ) + ( v yi ) (cid:3) N bin − v ( m j , t ) , (34)where we have used the shorthand notation i ∈ j todenote particles in the j -th bin, and L j to denote thelength of the j -th bin. The pressure field p ( m j , t ) = ρ ( m j , t ) T ( m j , t ) is obtained straightforwardly.The hydrodynamic fields, computed for individual re-alizations, exhibit a significant noise. To get rid of thenoise, all hydrodynamic profiles presented in the nextSection were obtained after an averaging over 100 MDsimulations with different initial conditions, correspond-ing to the same initial hydrodynamic fields and obtainedwith the procedure described in subsection B. D. MD simulations versus ideal solution
The following parameters were chosen for the simula-tions: σ = 1, ρ = 10 − , T = 1, N ′ = 5 × , and L y = 125. For convenience, the coefficient of normalrestitution r was chosen so that 1 − r = p π/ × − ,i.e. r = 0 . . . . . This choice of parameters sets l ≃ . × , L ′ x = 10 l and τ ≃ . × . Theevolution of the density field, as obtained in the simula-tions and as predicted by the ideal theory, is displayed -5 -4 -3 ρ m -5 -4 -3 ρ m m t/ τ =0 t/ τ =0.2 t/ τ =0.4t/ τ =0.6 t/ τ =0.7 t/ τ =0.8 FIG. 2: (Color online) Evolution of the density ρ ( m, t ) inthe Lagrangian frame at initial times. The circles: the re-sults of MD simulations. The solid line: the prediction ofthe analytic solution, first of Eqs. (21), back in the physi-cal variables. The dashed line: a numerical solution of thenon-ideal hydrodynamic equations (8)-(10) that account forthe gradient-dependent transport. The dashed line is indis-tinguishable from the solid line. in Fig. 2 for times up to t = 0 . τ , and in Fig. 3 forlater times, up to time t = 1 . τ . The figures showthat the ideal solution is in remarkable agreement withthe MD simulations up to times t ≃ . τ . At later times,when the density peak exceeds ≃ − , the ideal solutionstarts to deviate from the MD simulation in the neigh-borhood of m = 0. The actual density peak continuesto grow, but slower than predicted by the ideal solution.At time t = τ , when the ideal solution predicts the den-sity blowup in x = 0, the actual density ρ (0 , τ ) ≃ . ρ c = 2 / ( √ σ ) is reached at t ≃ . τ , see the last frame of Figure 3. Sufficiently farfrom m = 0, the ideal solution remains very accurate.(We checked that this statement remains true even be-yond the attempted singularity time: until the end of theMD simulations.)The gas velocity profiles, shown in Figs. 4 and 5 ina linear and logarithmic scale, respectively, are very ac-curately predicted by the ideal solution, Eq. (21), untillate times. Surprisingly, the excellent agreement remainseven at times greater than 0 . τ , when the density peakalready significantly deviates from the theoretical one.To be able to see the small deviations from the theory,we had to use, in Fig. 5, a logarithmic scale.Similarly, an inspection of the pressure profiles, seeFig. 6, shows an excellent agreement with the predictionof the ideal theory, p ( m, t ) = ρ T cos( m ), see the secondof Eqs. (22). Discrepancies of about 2% appear only atlate times, when the density is already about one half ofthe theoretically predicted value. At times close to τ , thepressure field, as found in the MD simulations, developsa dip close to x = m = 0, reaching a value about 15%lower than the theoretically expected value p (0) = 10 − . -3 -2 -1 ρ m -3 -2 -1 ρ m m t/ τ =0.9 t/ τ =0.93 t/ τ =1.055t/ τ =0.96t/ τ =0.98 t/ τ =1 closepacking FIG. 3: (Color online) Evolution of the density ρ ( m, t ) in the Lagrangian frame at late times. The m -axis is zoomed in,in order to focus on the high density region. The circles: the results of MD simulations. The solid line: the prediction ofthe first of Eqs. (21), back in the physical variables. The dashed line: a numerical solution of the non-ideal hydrodynamicequations (8)-(10) that account for the gradient-dependent transport. The dash-dotted line marks the close packing density ρ c = 2 / ( √ σ ). -1.5-1-0.50 v m -0.6-0.30 v m m t/ τ =0 t/ τ =0.4 t/ τ =1t/ τ =0.8t/ τ =0.9 t/ τ =0.984 FIG. 4: (Color online) Evolution of the longitudinal velocity v ( m, t ) in the Lagrangian frame. The circles: the results ofMD simulations. The solid line: the prediction of Eq. (22),back in the physical variables. The dashed line: a numer-ical solution of the non-ideal hydrodynamic equations (8)-(10) that account for the gradient-dependent transport. Thedashed line is indistinguishable from the solid line. A direct characterization of the attempted gas densityblowup is provided by the time history of the density at x = 0. The ideal solution predicts, see Eq. (24), that1 − r ρ ρ (0 , t ) = tτ . (35)This prediction is extremely well supported by the MD -3 -2 -1 - v m -3 -2 -1 - v m m t/ τ =0 t/ τ =0.4 t/ τ =1t/ τ =0.8t/ τ =0.9 t/ τ =0.984 FIG. 5: (Color online) Evolution of the (minus) longitudinalvelocity v ( m, t ) in the Lagrangian frame, in a semi-logarithmicscale. The circles: the results of MD simulations. The solidline: the prediction of Eq. (22), back in the physical vari-ables. The dashed line: a numerical solution of the non-ideal hydrodynamic equations (8)-(10) that account for thegradient-dependent transport. simulations in Fig. 7, until t ≃ . τ . The subsequentdeviation from the theory appears as saturation of thequantity 1 − p ρ /ρ (0 , t ) at the value 1 − p ρ /ρ c cor-responding to the close packing density ρ c . The sameFig. 7 also depicts a different quantity, 1 − p T (0 , t ). Inview of the theoretical expectation p (0 , t ) = ρ T = ρ ,this quantity is also expected to grow as t/τ , and Fig. 7 × -5 × -5 × -5 p m × -5 × -5 × -5 p m m t/ τ =0 t/ τ =0.4 t/ τ =1t/ τ =0.8t/ τ =0.9 t/ τ =0.984 FIG. 6: (Color online) Evolution of the gas pressure, obtainedusing the ideal equation of state p ( m, t ) = ρ ( m, t ) T ( m, t ), inthe Lagrangian frame. The circles: the results of MD simula-tions. The solid line: the prediction of the second of Eqs. (21),back in the physical variables. The dashed line: a numeri-cal solution of the non-ideal hydrodynamic equations (8)-(10)that account for the gradient-dependent transport. τ -[ ρ / ρ ( , t )] / , - T ( , t ) / t/ τ FIG. 7: (Color online) The time history of the density maxi-mum ρ (0 , t ) and the temperature minimum T (0 , t ). The sym-bols show the results of MD simulations: black circles depictthe quantity 1 − p ρ /ρ (0 , t ), green diamonds depict the quan-tity 1 − p T (0 , t ). The solid line is the prediction of Eq. (35):an immediate consequence of Eq. (24). The dashed line is anumerical solution of the non-ideal hydrodynamic equations(8)-(10) that account for the gradient-dependent transport.The inset zooms in at later times. indeed shows this growth until t ≃ . τ . The quantity1 − p T (0 , t ) also saturates at later times, but at a valueslightly different from 1 − p ρ /ρ c , see the inset of Fig. 7.This is consistent with the pressure deviation from ρ atvery late times.We also present, in Fig. 8, the time history of the totalenergy per particle, E ( t ) N ′ = 1 N ′ N ′ X i =1 | v | , (36)as found in the MD simulations, and compare it with the t/ τ E t o t / N ’ t/ τ -3-2.5-2-1.5-1-0.50 N ’ - d E / d t FIG. 8: (Color online) Decay of the total energy per particle E ( t ) /N ′ as measured in the MD simulations (the circles) andpredicted by the exact solution (the solid line). The timederivative ( N ′ ) − dE/dt , displayed in the inset (the notationfor the circles and line is the same), shows that the energydecay remains smooth at all simulated times. theoretical prediction, Eq. (28). Here the agreement isvery good at all times, with a 3% error at late times.The (numerical) time derivative of E ( t ), depicted in theinset of Fig. 8, remains smooth also close to the singular-ity time τ . Actually, this is not surprising, as the maincontribution to the thermal energy of the gas is made bythe peripheral gas (in the Lagrangian frame), which ishotter and more dilute than the gas in the region closeto m = 0. Furthermore, the main contribution to the ki-netic energy of macroscopic motion is again made by theperipheral gas (in the Lagrangian frame) which movesfaster than the gas in the region close to m = 0. Theperipheral gas continues to follow the ideal theory at all simulation times, and this explains the remarkable suc-cess of the ideal solution in predicting the total energyhistory.To conclude this Section, the MD simulations clearlyshow, over an extended period of time, the existence ofthe “attempted” density blowup regime. The ideal gran-ular hydrodynamics (IGHD) predicts very accurately thehydrodynamic profiles, observed in the MD simulations,up to times close to the attempted singularity. The den-sity field, as measured in the MD simulations, starts todeviate from the ideal theory at time t ≃ . τ . Some-what surprisingly, the rest of the hydrodynamic fieldscontinue to show good agreement with the theory un-til even closer to the attempted singularity time τ . Inthe following Section we will see that the agreement withtheory at later times improves significantly when the non-ideal hydrodynamic equations (8)-(10), that account forthe gradient-dependent transport, are employed. IV. NON-IDEAL HYDRODYNAMICS
To investigate the mechanism of breakdown of the idealtheory, we extended the hydrodynamic theory in two ρ ( m ,t ) (a) ρ ( m ,t ) (b) FIG. 9: (Color online) The density profiles for t/τ = 0 . t/τ = 0 .
966 (b) as obtained from the MD simulations(the circles), the analytic solution (the dashed-dotted line),the numerical solution of the dilute NIGHD equations (thedashed line), and the numerical solution of finite-density hy-drodynamic equations with the gradient-dependent transportterms neglected (the solid line). separate ways. In the first one we took into accountthe gradient-dependent transport: the heat diffusion andviscosity, but continued to assume the the gas is dilute.In the second one we neglected the gradient-dependenttransport but took into account finite density corrections.The hydrodynamic equations were solved numerically inLagrangian coordinates, using an accurate variable meshand variable time step solver [26].First, we solved (the Lagrangian form of) Eqs. (8)-(10)of non-ideal granular hydrodynamics (NIGHD). Theseequations account for the viscous and heat diffusionterms but still assume a dilute gas.The numerically obtained NIGHD profiles are pre-sented, together with the MD simulations and the idealanalytical solution, in Figs. 2-7. As expected, theNIGHD profiles coincide with the MD simulations andwith the ideal theory for early and intermediate times.At late times, the gradient-dependent transport terms become significant, and the NIGHD profiles approximatethe MD simulation results much better than the ideal the-ory. As the maximum density continues to increase (andfinally approaches the close packing density ρ c ), the di-lute NIGHD description ultimately breaks down. In Fig.3 it occurs at about t/τ ≃ . p → ρT (cid:20) πρ √ g ( ρ ) (cid:21) , Λ → Λ g ( ρ ) , (37)where g ( ρ ) = 1 − πρ √ (cid:16) − πρ √ (cid:17) is the equilibrium pair correlation function of hard disksat contact. In the dilute limit ρ → g = 1and recovers the ideal equation of state p = ρT .In Fig. 9 we compare four different results for the gasdensity: the MD simulations, the ideal analytical solu-tion, the numerical solution of the first type (the NIGHD-equations), and the numerical solution of the second type.It is clearly seen that, for the choice of parameters used inour MD simulations, the numerical solution of the secondtype is not as successful as that of the first type. Thiscould be expected, as for the times when the maximumdensity is still much smaller than the close packing den-sity, the finite-density corrections (which are of the orderof ρ/ρ c ) are still negligible. In contrast, the numericalresults from the NIGHD equations agree well with theMD simulations and show that, as the attempted densityblowup is approached, the heat conduction and viscosityeffects can become important when the gas density is stillsmall. V. SUMMARY AND DISCUSSION
Our MD simulations proved the existence of an at-tempted density blowup regime as described by an exactsolution of the ideal granular hydrodynamic equations.We found the ideal solution to be in remarkable quan-titative agreement with the MD simulations over an ex-tended time interval, but not too close to the attemptedsingularity. As the attempted singularity is approached,the exact solution breaks down. A more complete hy-drodynamic theory, that accounts for the heat diffusionand viscosity, but still assumes a dilute gas, continues0to agree with the MD simulations until the gas densitybecomes a fraction of the close packing density of disks.Let us put the results of this work into a more generalcontext of clustering instability of a freely cooling granu-lar flow. As we have already noted, the local properties ofthe density blow-up, exhibited by the exact solutions ofthe IGHD equations, are the same as the local propertiesof the density blow-up exhibited by a low Mach number flow of a freely cooling granular gas with the heat dif-fusion neglected [18]. A low Mach number flow emergeswhen the pressure balance sets in on a shorter time scalethan the temperature balance. In this case any local in-elastic cooling causes a (low Mach number) gas inflowinto the colder region so as to increase the local gas den-sity there and keep the pressure gradient (almost) zero.The resulting density instability develops on the back-ground of an (almost) homogeneous gas pressure. Thisis consistent with the MD simulations presented here, seeFig. 6, where the pressure in the vicinity of the densitymaximum hardly changes up to times very close to the at-tempted singularity time. For brevity we will call the lowMach number flow instability scenario Scenario 1. Sce-nario 1 first appeared in astrophysics and plasma physicsin the context of condensation instabilities in gases andplasmas that cool by their own radiation [29].As many as four additional hydrodynamic scenariosof clustering in a freely cooling granular gas have beendiscussed in the literature. We start with the pressureinstability scenario, or Scenario 2. It was discussed, inthe context of the granular clustering, by Goldhirsch andZanetti [2], although it was also known earlier to the as-trophysics and plasma physics communities, see Ref. [29]for a review. Scenario 2 is usually presented in the fol-lowing way. Let us consider a small local increase in thegas density. This increase causes an increase in the col-lisional energy loss. As a result, the gas pressure fallsdown, a gas inflow develops, causing a further densityincrease, and the process continues. Importantly, Sce-nario 2 assumes that the inelastic cooling time is muchshorter than the sound travel time. In other words, it isthe local temperature balance that sets in rapidly here,and the resulting pressure gradient drives the flow on arelatively slow time scale.As of present, there has been no detailed nonlinearanalysis behind Scenario 2. The (well established) linearstability theory of the homogeneous cooling state [3, 5]indicates that Scenarios 1 and 2 operate in two oppo-site limits: for sufficiently short and long perturbationwavelengths, respectively. The physics behind this is thefollowing. The characteristic cooling time due to the in-elastic collisions is independent of the length scale of theinitial perturbation, whereas the sound travel time scaleis proportional to it. As a result, when all other parame-ters are fixed, Scenario 1 corresponds to an intermediatewavelength limit of the clustering instability, while Sce-nario 2 corresponds to the long wavelength limit. (In theshort wavelength limit the homogeneous cooling state ofthe gas is stable, as the clustering instability is suppressed by the heat diffusion [3, 18].)Now let us consider Scenario 3 that also assumes along-wavelength limit. As the gas temperature falls downrapidly because of the inelastic cooling, the flow is de-scribable by the zero pressure (or flow by inertia) ap-proximation [10]. Were it not regularized by close pack-ing effects, such a flow would develop a finite-time densityblowup (of a different type than the low Mach numberflow) [10, 22]. If the compressional heating interferes ear-lier than the close packing effects, the pressure becomesrelevant again, and Scenario 3 gives way to the Scenario1 [16]. In the opposite case the late time dynamics ofthe system is describable by the Burgers equation [11].Which of the two regimes is realized in a particular set-ting depends on the initial conditions.Scenarios 1-3 do not invoke the shear mode instability,and so they can operate both in one-dimensional, andmulti-dimensional settings. On the contrary, Scenarios 4and 5 do invoke the shear mode, and so they are intrin-sically multi-dimensional (and intrinsically non-linear).Scenario 4 exploits the fact that the unstable shear modemay contribute, via a non-linear coupling, to the growthof the clustering mode. Obviously, the nonlinear couplingis the only hydrodynamic mechanism of driving the clus-tering mode if the system size is larger than the criticalsize for the shear mode instability, but smaller than thecritical size for the clustering mode instability. Further-more, numerical analysis, performed in Ref. [6], indi-cated that the nonlinear coupling plays a dominant rolein the initial density growth also in the case when thesystem size is comparable to the critical system sizes forthe clustering and shear instabilities. (More precisely,the wave number of the monochromatic test perturba-tion of the transverse velocity in Ref. [6] was withinthe instability regions of both the shear, and the clus-tering modes. However, twice the wave number alreadycame out of the instability region.) What happens inmuch larger systems is presently under investigation. Itturns out that well above the clustering mode instabil-ity threshold the nonlinear coupling with the shear modedoes not dominate the density growth (though it doesmake the theory more cumbersome). In the channel ge-ometry, that we adopted in this paper and in the previousworks [10, 11, 16, 17, 18], the shear mode is suppressed,and the clustering instability develops in its pure andsimplest form.Now we proceed to Scenario 5 [2] that exploits thefact that the unstable shear mode heats the gas in someregions. Scenario 5 assumes that this heating can bebalanced by the inelastic energy loss, rendering (quitea complicated) steady state. It is furthermore assumedthat this steady state is unstable with respect to smallperturbations, and it is this instability that causes thegranular clustering. We are unaware of a quantitativetheory that would support Scenario 5, or of any quantita-tive test of Scenario 5 in MD simulations or in numericalsolutions of hydrodynamic equations.To complete the comparison of the five hydrodynamic1scenarios of clustering we note that the only scenariosthat have addressed, up to date, a strongly nonlinearstage of the clustering process quantitatively is the lowMach number flow instability scenario (Scenario 1) [18]and the zero pressure scenario (Scenario 3) [10, 11]. Itis the consideration of a strongly nonlinear stage thatenables one to identify attempted finite-time densityblowups: prototypes of the dense granular clusters.Which results of this work will withstand a general-ization to more realistic granular flow conditions: for ex-ample, rotational degrees of freedom and tangential in-elasticity of collisions? Including the rotational degreesof freedom and tangential inelasticity of collisions intoa hydrodynamic description is possible under some lim-itations [13, 20]. Solving the corresponding nonlinearhydrodynamic equations analytically will of course bebeyond our reach. It is likely that, when the gradient-dependent transport is negligible, these nonlinear equa-tions will again exhibit a finite-time density blowup. In-deed, the development of closely packed granular clus-ters in a granular flow is a robust phenomenon. There-fore, it is natural to conjecture, based on results of thiswork, that more realistic granular clusters (those emerg- ing when the rotational degrees of freedom are taken intoaccount) will still be describable as regularized attempteddensity blowups.In summary, the results of the present work gives sup-port to the notion of a granular cluster as a regularizeddensity blowup of ideal granular hydrodynamic equa-tions, put forward in Refs. [10, 11, 15, 16, 17, 18].In more general terms, they present additional evidencethat granular hydrodynamics is a powerful and accuratequantitative theory of granular flows, especially once itis employed within its limits of applicability but, luckily,sometimes even beyond them.
Acknowledgments
Our work was supported by the Italian MIUR (PRINgrant n. 2005027808 003), by the Israel Science Founda-tion (grant No. 107/05) and by the German-Israel Foun-dation for Scientific Research and Development (GrantI-795-166.10/2003). [1] M. A. Hopkins and M. Y. Louge, Phys. Fluids A , 47(1991).[2] I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. , 1619(1993); I. Goldhirsch, M.-L. Tan, and G. Zanetti, J. Sci.Comp. , 1 (1993).[3] S. McNamara, Phys. Fluids A , 3056 (1993).[4] S. McNamara and W. R. Young, Phys. Rev. E , 5089(1996).[5] R. Brito and M. H. Ernst, Europhys. Lett. , 497 (1998).[6] J. J. Brey, M. J. Ruiz-Montero, and D. Cubero, Phys.Rev. E , 3150 (1999).[7] S. Luding and H. J. Herrmann, Chaos , 673 (1999).[8] T.P.C. van Noije and M.H. Ernst, Phys. Rev. E , 1765(2000).[9] X. B. Nie, E. Ben-Naim, and S. Y. Chen, Phys. Rev.Lett. , 204301 (2002).[10] E. Efrati, E. Livne, and B. Meerson, Phys. Rev. Lett. ,088001 (2005).[11] B. Meerson and A. Puglisi, Europhys. Lett. , 478(2005).[12] V. Garz´o, Phys. Rev. E , 021106 (2005).[13] I. Goldhirsch, Annu. Rev. Fluid Mech. , 267 (2003).[14] Y. Bromberg, E. Livne, and B. Meerson, in GranularGas Dynamics , edited by T. P¨oschel and N.V. Brilliantov(Springer, Berlin, 2003), p. 251; cond-mat/0305557.[15] D. Volfson, B. Meerson, and L. S. Tsimring, Phys. Rev.E , 061305 (2006).[16] I. Fouxon, B. Meerson, M. Assaf, and E. Livne, Phys. Rev. E , 050301(R) (2007).[17] I. Fouxon, B. Meerson, M. Assaf, and E. Livne, Phys.Fluids , 093303 (2007).[18] B. Meerson, I. Fouxon, and A. Vilenkin, arXiv:0708.3085[cond-mat.soft].[19] L. P. Kadanoff, Physics Today , 11 (1997).[20] N.V. Brilliantov and T. P¨oschel, Kinetic Theory of Gran-ular Gases (Oxford University Press, Oxford, 2004).[21] Ya. B. Zel’dovich and Yu. P. Raizer,
Physics of ShockWaves and High Temperature Hydrodynamic Phenom-ena , vol. 1 (Academic Press, New York, 1966).[22] G. B. Whitham,
Linear and Nonlinear Waves (Wiley,New York, 1974), Chapter 2.[23] P.K. Haff, J. Fluid Mech. , 401 (1983).[24] The assumption of a small Knudsen number also breaksdown close to the singularity. However, the finite heat dif-fusion invalidates the solution much earlier, see Ref. [17].[25] T. P¨oschel and T. Schwager,
Computational GranularDynamics (Springer, Berlin, 2005).[26] J.G. Blom and P.A. Zegeling, ACM Trans. Math. Soft-ware , 194 (1994).[27] N. F. Carnahan and K. E. Starling, J. Chem. Phys. ,635 (1969).[28] J.T. Jenkins and M.W. Richman, Phys. Fluids , 3485(1985); Arch. Rat. Mech. Anal. , 355 (1985).[29] B. Meerson, Rev. Mod. Phys.68