Table Of ContentAttempted density blowup in a freely cooling dilute granular gas: hydrodynamics
versus molecular dynamics
Andrea Puglisi1, Michael Assaf2, Itzhak Fouxon2, and Baruch Meerson2
1Dipartimento di Fisica, Universit`a La Sapienza, p.le Aldo Moro 2, 00185 Roma, Italy
2Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel
It has been recently shown (Fouxon et al. 2007) that, in the framework of ideal granular hydro-
dynamics(IGHD),an initially smooth hydrodynamicflowofagranulargascanproduceaninfinite
gas density in a finite time. Exact solutions that exhibit this property have been derived. Close to
thesingularity,thegranulargaspressureisfiniteandalmostconstant. Thisworkreportsmolecular
8 dynamics (MD) simulations of a freely cooling gas of nearly elastically colliding hard disks, aimed
0 at identifying the “attempted” density blowup regime. The initial conditions of the simulated flow
0 mimicthoseofoneparticularsolution oftheIGHDequationsthatexhibitsthedensityblowup. We
2 measure the hydrodynamic fields in the MD simulations and compare them with predictions from
n the ideal theory. We find a remarkable quantitative agreement between the two over an extended
a time interval, proving theexistence of the attempted blowup regime. As theattempted singularity
J isapproached,thehydrodynamicfields,asobservedintheMDsimulations,deviatefromthepredic-
5 tionsof theideal solution. Toinvestigate themechanism of breakdownof theideal theory nearthe
1 singularity,weextendthehydrodynamictheorybyaccountingseparatelyforthegradient-dependent
transport and for finite density corrections.
]
t PACSnumbers: 45.70.Qj
f
o
s
t. I. INTRODUCTION diffusion and viscosity. Fouxon et al. derived a class of
a exact solutions to the ideal equations, for which an ini-
m
tially smoothflow develops a finite-time density blowup.
Spontaneous clustering of particles in granular gases
- Close to the blowup time τ, the maximum gas density
d has attracted much recent interest [1, 2, 3, 4, 5, 6, 7, 8,
exhibits a power law behavior (τ t) 2. The velocity
on 9cl,u1s0te,r1in1,g1i2n]s.taAbsiloittyheorfpaaftrteeerlny-fcooromliinnggginrastnaublailritgieass,hthaes gradient blows up as ∼ −(τ −∼t)−1−, wh−ereas the veloc-
c ity itself remains continuous and forms a cusp, rather
served as a sensitive probe of theoretical modeling and,
[ than a shock discontinuity, at the singularity. The gas
first of all, of the Navier-Stokes granular hydrodynam-
temperaturevanishesatthesingularity,but the pressure
2 ics (GHD). Although the formal criteria of its validity
v remains finite and almost constant. Extensive numeri-
may be quite restrictive, see below, GHD has a great
3 cal simulations with the ideal hydrodynamic equations
power, sometimes far beyond its formal validity limits
3 showed that the singularity, exhibited by the exact solu-
[13],inpredictingahostofcollectivephenomenaingran-
4 tions, is universal, as it develops for quite general initial
3 ular flows, such as shocks, vortices and clusters. In the
conditions [16, 17]. The reason behind this universality
. recent years, GHD has been applied to a variety of non-
9 is in the fact that the sound travel time through the re-
stationary dilute granular flows [10, 14, 15, 16, 17, 18].
0 gion of the developing singularity is much shorter than
Non-stationary flows are both appealing and challeng-
7
thecharacteristicinelasticcoolingtimeofthegasinthat
0 ing for continuum modeling of granulardynamics. As in
region. As a result, the pressure gradient (almost) van-
: otherareasofcontinuummodeling,thisisespeciallytrue
v ishes in the singularity region, and the local features of
when a non-stationary flow develops a finite-time singu-
i this isobaric singularity become essentially independent
X laritiy [19]. Examples are provided by the finite-time
ofhowtheflowwasinitiatedandhowitbehavesatlarge
r blowup of the gas density: at zero gravity [10, 16, 17]
a distances from the singularity. This singularity is of the
(as described by the ideal GHD discussed below, IGHD
same type as the one that develops, in the framework of
from now on), and at finite gravity [15] as described by
the IGHD equations,ina generallow Mach number flow
the more complete, non-ideal GHD. Of course, a den-
of a freely cooling granular gas [18].
sity blowupina gaswith finite-size particlescanonly be
anintermediateasymptotics,astheblowupisultimately Here we perform molecular dynamics (MD) simula-
arrested: either by close-packing effects [11], or by the tions of a freely cooling gas of nearly elastically colliding
gradient-dependent transport [18]. Still, the attempted harddisks,aimedatidentifyingthe“attempted”density
blowup regimes, signaling the development of high den- blowup regime, predicted by the ideal analytical solu-
sity regions in the gas, are fascinating and worth a de- tions [16,17]. We simulate afreely evolvingdilute gasof
tailed investigation. One such regime has been recently nearlyelasticallycollidingharddisksinanarrowchannel
addressed by Fouxon et al. [16, 17]. They dealt with with perfectly elastic side walls. In this geometry both
a macroscopically one-dimensional, freely cooling, dilute the clustering mode inthe transversedirections,and the
granular flow in the framework of ideal GHD that ne- shear mode are suppressed (see Refs. [2, 3, 10, 11] for
glects the gradient-dependenttransporteffects: the heat detailed criteria). As a result, the coarse-grained,or hy-
2
drodynamic flow depends only on the longitudinal coor- this Section that the local gas density ρ is much smaller
dinate along the channel (and time), as it was assumed than the close-packing density of disks ρ =2/(√3σ2):
c
in Refs. [16, 17]. We choose the initial conditions of the
MD simulations so that the coarse-grained density, ve-
locity and temperature fields are those producing one of ρσ2 1. (3)
≪
the exact blowup solutions of the ideal GHD equations.
Then we follow the time history of the hydrodynamic
The strong inequalities (2) and (3) need to be verified a
fields in the MD simulations and compare it with that
posteriori, once the hydrodynamic problem in question
predicted by the idealexact solutionandwith numerical
is solved. The strong inequalities (1)-(3) enable one to
solutions of non-ideal GHD equations.
employ the well-established equations of Navier-Stokes
The remainder of this paper is organized as follows.
granular hydrodynamics (see, e.g. [13, 20]) that deal
In Section II we briefly summarize the ideal GHD anal-
withthreecoarse-grainedfields: themassdensityρ(x,t),
ysis [16, 17] of the density blowup: we present the ideal
the meanflowvelocityv(x,t) andthe granulartempera-
GHD equations, one of their exact solutions, its main
ture T(x,t). In a sufficiently narrow channel these fields
features and expected limits of its validity. In Section
depend only on the longitudinal coordinate x, and the
III wedescribe ourMD simulationsandcomparethe hy-
hydrodynamic equations become
drodynamic quantities, computed from the simulations,
with the exact solution of the ideal GHD equations. We
findthattheexactsolutionis inremarkablequantitative ∂ρ ∂(ρv)
agreement with the MD simulations over an extended + =0, (4)
∂t ∂x
timeinterval,provingtheexistenceoftheattemptedden-
∂v ∂v ∂(ρT) ∂ ∂v
sity blowup regime. As the attempted singularity is ap- ρ +v = +ν √T , (5)
0
∂t ∂x − ∂x ∂x ∂x
proached, the hydrodynamic fields, as observed in the (cid:18) (cid:19) (cid:18) (cid:19)
MD simulations, deviate from the predictions of the ex- ∂T +v∂T = T∂v ΛρT3/2+ κ0 ∂ √T∂T
actsolution. Toinvestigatethemechanismofbreakdown ∂t ∂x − ∂x − ρ ∂x ∂x
(cid:18) (cid:19)
of the ideal solution, we extend the hydrodynamic the- ν √T ∂v 2
ory, in Section IV, in two separate ways. In the first one + 0 , (6)
ρ ∂x
we take into account the gradient-dependent transport: (cid:18) (cid:19)
the heat diffusion and viscosity, but continue to assume
that the gas is dilute. In the second one we neglect the
where Λ = √π(1 r2)σ, ν = (2√πσ) 1 and κ =
gradient-dependent transport but take into account, in − 0 − 0
2/(√πσ), see e.g. Refs. [6, 20]. Equations (4)-(6) differ
a semi-phenomenological way, finite density corrections.
from the hydrodynamic equations for a gas of elastically
Section V summarizes our results and puts them into a
colliding disks only by the presence in Eq. (6) of the in-
more general context of hydrodynamic scenarios of clus-
elastic loss rate term ΛρT3/2, that describes the pro-
tering in freely evolving granular gases. −
portionality of the energy loss per particle to the num-
ber of particle collisions per unit time (proportional to
ρT1/2) and to the energy loss per collision (proportional
II. HYDRODYNAMIC THEORY AND DENSITY
toT). Thisadditionalterm,however,bringsawholenew
BLOWUP
physics (and mathematics) into the problem.
A. Ideal granular hydrodynamics and exact Let us rewrite Eqs. (4)-(6) in dimensionless variables.
solution We will measure the gas density, temperature and the
velocityintheunitsofρ ,T /2and T /2,respectively,
0 0 0
We consider a two-dimensional granular gas of identi- where ρ0 and T0 are some characteristic values of the
p
calhardandsmoothdisks with diameter σ andmassset initial density and temperature. The time and distance
to unity, and adopt a simple model where the inelastic will be measured in the units of
particle collisions are characterized by a constant coeffi-
cient of normal restitution r. Throughout this paper we
4 T
will only deal with nearly elastic collisions, τ = and l =τ 0 , (7)
Λρ √T 2
0 0 r
1 r 1, (1)
− ≪
respectively. As one can see from Eq. (6), τ is the char-
and assume a very small Knudsen number:
acteristic cooling time of the gas due to the collisional
l /L 1. (2) energyloss, while l is the characteristicdistance a sound
free
≪ wave travels during time τ. The numerical factors in
Herel isthe meanfreepathofthe particles,andLis Eqs.(7) are chosenfor convenience. We will keepidenti-
free
thecharacteristiclengthscaleofthehydrodynamicfields calnotationforthe rescaledandphysicalquantities,and
thatmaydependontime. Inaddition,wewillassumein take care that no confusion arises. Using the rescaled
3
1
quantities, we rewrite Eqs. (4)-(6) as
0.8
∂ρ ∂(ρv)
∂t + ∂x =0, (8) ρ00.6
∂v ∂v ∂(ρT) 1 r2 ∂ ∂v ρ/ 0.4
ρ +ρv = + − √T ,(9)
∂t ∂x − ∂x 4√2 ∂x ∂x 0.2
(cid:18) (cid:19)
∂T ∂T ∂v 00
+v = T 2√2ρT3/2
∂t ∂x − ∂x−
2 -1
+1−r2 ∂ √T∂T + (1−r2)√T ∂v 2, (10) 1/2)
√2ρ ∂x(cid:18) ∂x(cid:19) 4√2ρ (cid:18)∂x(cid:19) (T/0-2
Let us assume that the characteristic magnitudes of the v/ -3
rescaled hydrodynamic fields, and of their spatial and
-4
temporalderivatives,areoforderunity (this assumption 0 1 2 3 4 5 6 7 8 9 10
x/l
needs to be checked a posteriori). Then we can neglect
the viscous and thermal conduction terms, as they scale
FIG. 1: The initial values of the hydrodynamic fields in the
as 1 r2 1, and arrive at the IGHD equations:
− ≪ exampleofattempteddensityblowupconsideredinthiswork.
∂ρ ∂(ρv) ∂v ∂v ∂(ρT) Shown are therescaled initial density and velocity of thegas
+ =0, ρ +v = ,(11) [see Eqs. (15) and (17)] versus the rescaled Eulerian coordi-
∂t ∂x ∂t ∂x − ∂x
(cid:18) (cid:19) natex/l along thechannel. Onlytherighthalfofthesystem
∂T +v∂T = T∂v 2√2ρT3/2. (12) is shown. The values of ρ0, T0 and l, used in our molecu-
∂t ∂x − ∂x − lar dynamic simulations, are presented in the beginning of
subsection III.D.
TheseequationswereinvestigatedinRefs.[16,17],where
a family of exact analytic solutions was derived. Here
we will consider a representative and simple particular Note that the infinite Eulerian interval < x < +
solutionthatevolvesfromthefollowinginitialconditions: corresponds, in this example, to a finit−e∞interval of m∞:
ρ(x,t=0)=cosh−1x, T(x,t=0)=2, (13) −π/2<m<π/2,thatistoafinite(rescaled)totalmass
of the gas, equal to π. In the Lagrangian coordinates,
v(x,t=0)= 2arcsin(tanhx). (14)
the ideal equations (11) and (12) are
−
To remind the reader, we are using rescaled variables
∂ 1 ∂v ∂v ∂p
here. Back in the physical variables, the initial profiles = , = , (19)
are ∂t ρ ∂m ∂t −∂m
(cid:18) (cid:19)
ρ ∂p ∂v
ρ(x,0)= 0 , (15) = 2pρ 2√2p3/2ρ1/2, (20)
cosh(x/l) ∂t − ∂m −
T(x,0)=T0, (16) where the dilute gas pressure p = ρT has been used in-
x stead of the temperature. The exact solution in the La-
v(x,0)= 2T arcsin tanh , (17)
− 0 l grangiancoordinates is the following [16, 17]:
p h (cid:16) (cid:17)i
wherelisdefinedinEq.(7). Thatis,att=0thedensity cosm
profile has a maximum ρ and width l, the temperature ρ(m,t)= , p(m,t)=2cosm, (21)
0 (1 tcosm)2
T is uniform, and there is an inflow of the gas towards −
0 v(m,t)= 2m+2tsinm. (22)
the origin with v(x ,t = 0) = π T0/2. The −
→ ±∞ ∓
initial scale of variation of the fields, l, is by a factor
p As x(m,t) can be calculated explicitly [17],
1/(1 r2) greater than the mean free path of the gas,
−
justifying the use of the hydrodynamic description. Fur-
1 1+sinm
thermore,both the magnitudes, and the spatial scalesof x(m,t)= ln 2tm+t2sinm,(23)
2 1 sinm −
the rescaled fields are of order unity which justifies, at (cid:18) − (cid:19)
least for finite times, the use of the ideal equations (11)
equations (21)-(23) describe the time history of the hy-
and (12). Figure 1 shows the initial density and velocity
drodynamic fields in the x-coordinate in a closed para-
fields of the flow in the Eulerian coordinate.
metric form. Let us consider some important features of
Now we go overfrom the Euleriancoordinatex to the
this simple exact solution.
Lagrangianmass coordinate m= xρ(x,0)dx, see e.g.
0 ′ ′
[21]. For the initial density profile (13), the Eulerian
R
coordinate x is related to the Lagrangian coordinate m B. Density blowup and its properties
as follows:
1 1+sinm AmomentarylookatthefirstofEqs.(21)revealsthat
x(m,t=0)= ln . (18)
2 1 sinm the density at the origin blows up in a finite time. Back
(cid:18) − (cid:19)
4
in the dimensional variables one obtains algebra,
ρ
ρ(x=0,t)= 0 . (24) Eth(t) ∞ ∞
(1 t/τ)2 = ρTdx = 2 ρTdx
L
− y Z−∞ Z0
π/2
Thegastemperaturebecomeszeroatthesingularity,the = 2ρ T l [1 (t/τ)cosm]2dm
velocity gradient ∂v(x,t)/∂x blows up as (τ t) 1, 0 0 −
∼ − − − Z0
whereas the velocity itself remains continuous and forms 2
t π t
a cusp. This behavior at singularity is quite different = ρ T l π 4 + . (26)
0 0
from that of the free flow singularity (that one observes " − (cid:18)τ(cid:19) 2 (cid:18)τ(cid:19) #
for the same initial velocity profile but a zero gas pres-
Similarly, for the kinetic energy of the macroscopic mo-
sure). Inparticular,forthe freeflow singularitythe den-
sity blows up as (1 t/τ) 1, while the velocity develops tion we find
−
−
a shock discontinuity [22]. See Refs. [16, 17] for more
E (t) π3 t π t 2
differences between the two types of the singularities. kin =ρ T l 4 + . (27)
0 0
It was found in Refs. [16, 17] that the finite time Ly "12 − (cid:18)τ(cid:19) 2 (cid:18)τ(cid:19) #
blowup of the gas density and of the velocity gradient is
not a consequence of specially chosen initial conditions. Summing the two and using Eq. (25), we obtain [17]
Rather, it appears, in the framework of the ideal equa-
tions (11) and (12) [or, equivalently, of Eqs. (19) and π2 8 t t 2
E(t)=NT 1+ + . (28)
(20)], for quite general initial conditions. A robust lo- 0" 12 − π (cid:18)τ(cid:19) (cid:18)τ(cid:19) #
cal feature of this singularity is a finite, non-zero value
of the gas pressure at the time of the density blowup. As time increases, E(t) is monotone decreasing for t
Thesingularityisuniversalbecausethesoundtraveltime τ. We also observe that t = τ is a regular point o≤f
throughit is muchshorterthanthe characteristicinelas- E(t), where nothing dramatic happens. Note that the
tic cooling time of the gas. As a result, the pressure parabolic-in-time law of the energy decay is quite differ-
gradient (almost) vanishes in the singularity region, and entfromHaff’slawE(t)=E(0)/(1+2t/τ)2 obtainedfor
the local features of this isobaric singularity become in- freely cooling homogeneous granular gas with density ρ
0
dependent of the flow details at large distances. [23].
As an infinite density cannot be reached in a gas with
finite-size particles, it is clear that, sufficiently close to
theattemptedsingularity,someoftheassumptionsmade C. Breakdown of ideal theory
on the way to the ideal equations (11)-(12) break down.
However, independently of the precise way of regulariz- Foraone-dimensionalflow,the applicabilityofthe so-
ing the infinite density, an attempted singularity implies lution (21)-(23) is determined by the applicability of the
the formation of a region with a very high density. Fur- ideal equations (11)-(12) that it solves exactly. Analysis
thermore, the ideal solution should accurately describe inRef.[17]showsthat,sufficientlyclosetotheattempted
the evolutionofthe systemfortimes nottoocloseto the singularity,theidealequationsbecomeinvalid. Thishap-
attempted density blowup time τ. Before we look into pens because of one of two reasons (or both): either the
wheretheidealtheorybreaksdown,letusconsidersome gas becomes dense, so that criterion (3) breaks down, or
global characteristicsof the exactsolution. For these the the heat diffusion becomes important, invalidating the
singularity time t=τ turns out not to be special. First, idealequations(11)-(12). Thetimet atwhichtheideal
br
we note that the solution describes a gas with a (con- theory breaks down can be estimated as follows [17]:
stant) finite number N of particles, given by
t
1 br max( ρ σ2,1 r). (29)
∞ 2√2πLy − τ ∼ 0 −
N =L ρ(x)dx=πρ L l= , (25)
y 0 y (1 r2)σ p
Therefore, the “bottleneck” for the validity of the equa-
Z−∞ −
tions is set by the initial conditions: if the maximum
where Ly is the channel width. Note that N is indepen- is determined by the first (correspondingly, the second)
dent of ρ : for larger ρ the initial density profile has
0 0 term in the right hand side of Eq. (29), the ideal equa-
a higher peak but a smaller width, so that N remains
tions become invalid because of the finite gas density
constant. Another global characteristics of the solution
(correspondingly, the finite heat diffusion). As each of
is the total energy of the gas
the two terms is very small by assumption, the solution
is expected to break down only close to the attempted
∞ ρv2 singularity [24].
E(t)=L +ρT dx.
y
2 Theone-dimensionalsolutionmayalsobecomeinvalid
Z−∞(cid:18) (cid:19) becauseofinstabilitywithrespecttosmallinitialpertur-
Forthethermalpartoftheenergywefind,afterasimple bations that are inevitably present in MD simulations.
5
Numerical solutions of the hydrodynamic equations, re- withaconstantcoefficientofnormalrestitutionr [0,1]:
∈
ported in Refs. [16, 17], strongly suggest that the ideal
solution is stable with respect to small longitudinal per- 1+r
turbations. This does not exclude possible instability vi′ =vi− 2 (g·σˆ)σˆ, (30)
withrespecttosmalltransverse perturbations. Theonly 1+r
available analytic result here is the one obtained from vj′ =vj + 2 (g·σˆ)σˆ, (31)
the stability condition for a homogeneous cooling state,
see Refs. [2, 3, 10, 11]. That stability criterion comes where g = vi vj and σˆ is the unit vector joining the
−
from a competition between the (destabilizing) inelastic centersofthetwodisks. Thehard-coreinteractionsmake
cooling and the (stabilizing) heat diffusion and viscosity possible the following optimization of the algorithm. It
in the transverse direction. The stability criterion de- is sufficient to calculate the first collision times of all
mands that L be less than a threshold value depending particles and then select the absolute first one. The
y
on1 r,ρ andσ. Thestabilityproblemforthestrongly system is freely evolved up to that time, then the col-
0
−
inhomogeneousandtime-dependent exactsolutionisob- lision is resolved, and a new list of collision times is
viously more complicated, and its complete analytic so- computed. Withstandardoptimizationtechniquesofthe
lution does not seem feasible. It is therefore important search procedure it is possible to achieve fast computa-
that our MD simulations, presented in the next Section, tion times [25]. Nevertheless, the time performance is
strongly suggest that, for sufficiently narrow channels, proportionalto the number of collisions occurred, so the
no instability in the transverse direction occurs for the ratio between the physical time and the cpu-time goes
time-dependent flow we are working with. down when dense clusters emerge in the system.
Letussummarizethemaintheoreticalpredictions. For
the initial conditions (13) and (14), a nonlinear time-
dependentflowsetsin,describedbythe idealexactsolu- B. Initial conditions
tion: Eqs. (21)-(23). This flow “attempts” to develop a
density blowup. However, close to the attempted singu-
TheinitialpositionandvelocityofeachoftheN disks
′
larity one (or both) of the two factors: the finite density
arechosenrandomlywithprobabilitydistributionscorre-
and the heat diffusion, invalidates the solution. The rel-
spondingtothedesiredinitialhydrodynamicfields. This
ative importance of the two factors is determined, via
wasimplementedwith the followingprocedure. Foreach
Eq. (29), by the initial conditions.
disk i
1. thelongitudinalpositionx ischosenwithprobabil-
i
ity proportionalto ρ(x,0) from Eq. (15) for x 0,
III. MOLECULAR DYNAMICS SIMULATIONS using an acceptance/rejection method: ≥
A. General (a) a random x -position is generated with uni-
i
form probability on the interval [σ/2,L
′x −
σ/2],
To test the theoretical predictions, we performed MD
simulationsofafreecoolinggranulargasinanarrowtwo- (b) a random number z with uniform probability
dimensional channel. The initial conditions correspond on[0,max ρ(x,0) ]iscomparedwithρ(x ,0),
i
{ }
to hydrodynamic profiles (13) and (14) and satisfy the
(c) if z < ρ(x ,0) the position is accepted, other-
stronginequalities(1)-(3). Accordingtothetheory,they i
wise the procedure is repeated from 1a,
are expected to generate the nonlinear time-dependent
flow described by Eqs. (21)-(23). Our MD simulation
calculates the evolution of a gas of N′ identical inelastic 2. thenthe verticalpositionyi ischosenwithuniform
harddisks ofunit mass,with diameterσ, ina channelof probability on the interval [σ/2,Ly σ/2],
−
dimensions L = L /2 and L . As the expected hydro-
′x x y
dynamic flow is symmetric with respect to x = 0, only 3. a non-overlapcheck is performed: the distance be-
one half of the system, x [0,Lx/2], is simulated, so tweenthe diskcenter(xi,yi)andallthepreviously
N = N/2. Each wall of (∈the one half of) the channel placed disk centers must be greater then σ: if the
′
is solid and reflects elastically the disks colliding with it. conditionisnotsatisfied,theprocedureisrepeated
The particles move freely until a collision (“event”) oc- from 1a,
curswhentwodisksiandj findthemselvesatadistance
equaltoσ. Thecollisionisresolvedinstantaneously,leav- 4. the velocitycomponentsvx andvy arechosenfrom
i i
ing the positions of the particles unaltered and updat- a Gaussian distribution with zero mean and vari-
ing their velocities from (vi,vj), before the collision, to ance equal to T0; then the longitudinal component
(v ,v ), after the collision. The update rule conserves is shifted by an amount v(x,0) from Eq. (17) with
i′ j′
thetotalmomentumandreducesthetotalkineticenergy, x 0.
≥
6
C. Lagrangian coordinate and hydrodynamic fields
10-3 t/τ=0 t/τ=0.2 t/τ=0.4
We verified that, for our choice of the channel dimen-
sions, the gas remained homogeneous in the y-direction. ρ 10-4
Byvirtueofthisobservation,itwassufficienttodealwith
one-dimensional hydrodynamic profiles depending on x. 10-5
For a direct comparison with the analytical solution of
the IGHD, the hydrodynamic profiles were obtained us- 10-3 t/τ=0.6 t/τ=0.7 t/τ=0.8
ingauniformbinningintheLagrangianmasscoordinate.
UsingthesamenotationasinSectionI,wedefinetheLa- ρ 10-4
grangianmass interval for the simulated flow as [0,π/2],
whereπ/2correspondsto(onehalfof)thetotalgasmass,
-5
N =N/2. Letn bethenumberofbinschosentosam- 10
′ bin
plethehydrodynamicprofilesandN =N /n bethe 0 0.5 1 1.50 0.5 1 1.50 0.5 1 1.5
bin ′ bin m m m
average number of particles per bin. At a given time t
allparticlesareorderedsothatx <x ,i [0,N 1].
i i+1 ′ FIG. 2: (Color online) Evolution of the density ρ(m,t) in
∈ −
Then each bin j [1,nbin] has its leftmost border at the Lagrangian frame at initial times. The circles: the re-
∈
x[(j 1)Nbin] anditsrightmostborderatx[jNbin 1]. These sults of MD simulations. The solid line: the prediction of
bins−are non-uniform in the x coordinate, bu−t are uni- the analytic solution, first of Eqs. (21), back in the physi-
forminthe Lagrangianmasscoordinateaseachcontains cal variables. The dashed line: a numerical solution of the
the same mass N . The position of the j-th bin is non-ideal hydrodynamic equations (8)-(10) that account for
bin
the gradient-dependent transport. The dashed line is indis-
π(j 1/2) tinguishable from thesolid line.
m = − .
j
2n
bin
All particles belonging to the j-th bin contribute to the
in Fig. 2 for times up to t = 0.8τ, and in Fig. 3 for
value of the hydrodynamic fields:
later times, up to time t = 1.055τ. The figures show
N that the ideal solution is in remarkable agreement with
bin
ρ(mj,t)= , (32) theMDsimulationsuptotimest 0.9τ. Atlatertimes,
LyLj whenthedensitypeakexceeds 1≃0 2,theidealsolution
vx ≃ −
i j i starts to deviate from the MD simulation in the neigh-
v(mj,t)= ∈ , (33)
N borhood of m = 0. The actual density peak continues
P bin
(vx)2+(vy)2 to grow,but slower than predicted by the ideal solution.
T(mj,t)= i∈j i i v2(mj,t),(34) At time t =τ, when the ideal solution predicts the den-
2N −
P (cid:2) bin (cid:3) sity blowup in x = 0, the actual density ρ(0,τ) 0.2.
≃
where we have used the shorthand notation i j to The close-packing density ρ = 2/(√3σ2) is reached at
c
∈
denote particles in the j-th bin, and L to denote the t 1.05τ, see the last frame of Figure 3. Sufficiently far
j ≃
length of the j-th bin. The pressure field p(m ,t) = from m = 0, the ideal solution remains very accurate.
j
ρ(m ,t)T(m ,t) is obtained straightforwardly. (We checked that this statement remains true even be-
j j
The hydrodynamic fields, computed for individual re- yondtheattemptedsingularitytime: untiltheendofthe
alizations, exhibit a significant noise. To get rid of the MD simulations.)
noise, all hydrodynamic profiles presented in the next The gas velocity profiles, shown in Figs. 4 and 5 in
Section were obtained after an averaging over 100 MD a linear and logarithmic scale, respectively, are very ac-
simulations with different initial conditions, correspond- curately predicted by the ideal solution, Eq. (21), until
ing to the same initialhydrodynamicfields andobtained latetimes. Surprisingly,theexcellentagreementremains
with the procedure described in subsection B. even at times greater than 0.9τ, when the density peak
already significantly deviates from the theoretical one.
To be able to see the small deviations from the theory,
D. MD simulations versus ideal solution we had to use, in Fig. 5, a logarithmic scale.
Similarly, an inspection of the pressure profiles, see
The following parameters were chosen for the simula- Fig. 6, shows an excellent agreementwith the prediction
tions: σ = 1, ρ0 = 10−4, T0 = 1, N′ = 5×104, and oftheidealtheory,p(m,t)=ρ0T0cos(m),seethe second
Ly = 125. For convenience, the coefficient of normal of Eqs. (22). Discrepancies of about 2% appear only at
restitution r was chosen so that 1 r2 = π/2 10 2, late times, when the density is already about one half of
−
− ×
i.e. r = 0.99371367.... This choice of parameters sets thetheoreticallypredictedvalue. Attimesclosetoτ,the
p
l 1.27324 106, L = 10l and τ 1.8006 106. The pressure field, as found in the MD simulations, develops
≃ × ′x ≃ ×
evolution of the density field, as obtained in the simula- a dip close to x = m = 0, reaching a value about 15%
tions and as predicted by the ideal theory, is displayed lower than the theoretically expected value p(0)=10 4.
−
7
-1
10
t/τ=0.9 t/τ=0.93 t/τ=0.96
-2
10
ρ
-3
10
0
10
close
-1
10 packing
ρ
-2
10
-3
10 t/τ=0.98 t/τ=1 t/τ=1.055
0 0.4 0.8 0 0.4 0.8 0 0.4 0.8
m m m
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 of
the first of Eqs. (21), back in the physical variables. The dashed line: a numerical solution of the non-ideal hydrodynamic
equations (8)-(10) that account for the gradient-dependent transport. The dash-dotted line marks the close packing density
ρc =2/(√3σ2).
0
0
10
-0.5
-1
10
v -1 -v
-2
-1.5 t/τ=0 t/τ=0.4 t/τ=0.8 10 t/τ=0 t/τ=0.4 t/τ=0.8
-3
100
10
0
-1
10
v-0.3 -v10-2
-0.6 t/τ=0.9 t/τ=0.984 t/τ=1 10-3 t/τ=0.9 t/τ=0.984 t/τ=1
0 0.4 0.8 1.2 0 0.4 0.8 1.2 0 0.4 0.8 1.2 0 0.4 0.8 1.2 0 0.4 0.8 1.2 0 0.4 0.8 1.2
m m m m m m
FIG. 4: (Color online) Evolution of the longitudinal velocity FIG. 5: (Color online) Evolution of the (minus) longitudinal
v(m,t) in the Lagrangian frame. The circles: the results of velocityv(m,t)intheLagrangianframe,inasemi-logarithmic
MD simulations. The solid line: the prediction of Eq. (22), scale. The circles: the results of MD simulations. The solid
back in the physical variables. The dashed line: a numer- line: the prediction of Eq. (22), back in the physical vari-
ical solution of the non-ideal hydrodynamic equations (8)- ables. The dashed line: a numerical solution of the non-
(10) that account for the gradient-dependent transport. The ideal hydrodynamic equations (8)-(10) that account for the
dashed line is indistinguishable from thesolid line. gradient-dependenttransport.
A directcharacterizationof the attempted gas density
simulations in Fig. 7, until t 0.9τ. The subsequent
blowup is provided by the time history of the density at ≃
deviation from the theory appears as saturation of the
x=0. The ideal solution predicts, see Eq. (24), that
quantity 1 ρ /ρ(0,t) at the value 1 ρ /ρ cor-
0 0 c
− −
1 ρ0 = t. (35) rFeisgp.o7ndalisnogdteoppitchtes aclodsieffepreancktiqnugadnetintsyi,ty1 ρcp. TT(h0e,ts)a.mIne
− ρ(0,t) τ −
r view of the theoretical expectation p(0,t) = ρ T = ρ ,
0 0 0
p
This prediction is extremely well supported by the MD this quantity is also expected to grow as t/τ, and Fig. 7
8
2
0
9×10-5 -0.5
dt -1
p6×10-5 1.5 dE/-1.5
-1N’ -2
3×10-05 t/τ=0 t/τ=0.4 t/τ=0.8 E/N’tot 1 -2-.530 0.t5/τ 1
9×10-5
0.5
p6×10-5
3×10-5 00 0.5 1
t/τ=0.9 t/τ=0.984 t/τ=1 t/τ
0
0 0.4 0.8 1.2 0 0.4 0.8 1.2 0 0.4 0.8 1.2
m m m FIG. 8: (Color online) Decay of the total energy per particle
′
E(t)/N asmeasuredintheMDsimulations(thecircles)and
FIG.6: (Coloronline)Evolutionofthegaspressure,obtained predicted by the exact solution (the solid line). The time
using the ideal equation of state p(m,t) =ρ(m,t)T(m,t), in derivative (N′)−1dE/dt, displayed in the inset (the notation
theLagrangian frame. Thecircles: theresultsofMDsimula- for the circles and line is the same), shows that the energy
tions. Thesolidline: thepredictionofthesecondofEqs.(21), decay remains smooth at all simulated times.
back in the physical variables. The dashed line: a numeri-
calsolutionofthenon-idealhydrodynamicequations(8)-(10)
that account for thegradient-dependenttransport. theoretical prediction, Eq. (28). Here the agreement is
very good at all times, with a 3% error at late times.
The (numerical) time derivative of E(t), depicted in the
1
1/20,t)0.8 iintysettimofeFτig..A8,ctruemalalyi,nsthsimsoisotnhoatlssuorcplroisseintgo,tahsetshinegmulaairn-
T( contributionto the thermalenergyof the gasis made by
1/20,t)], 1-00..46 1 ththooetmtper=erai0pn.hdFemruaroltrhgeeardsmil(uoinrtee,tththheaenLmatgahrienangcgaoisnatnirnifbrtuahtmeioern)e,gtiowonhthiccehloksiies-
ρ( 0.95 netic energyofmacroscopicmotionisagainmadeby the
ρ1-[/00.2 0.90.9 0.95 1 1.05 1.1 pfaesrtieprhethraalngtahse(ginasthine tLhaegrraegnigoinanclforsaemteo) mwh=ich0.moTvhees
0 t/τ peripheral gas continues to follow the ideal theory at all
0 0.5 1
t/τ simulation times, and this explains the remarkable suc-
cess of the ideal solution in predicting the total energy
FIG. 7: (Color online) The time history of thedensity maxi- history.
mumρ(0,t)andthetemperatureminimumT(0,t). Thesym- To conclude this Section, the MD simulations clearly
bols show the results of MD simulations: black circles depict
show, over an extended period of time, the existence of
thequantity1−pρ0/ρ(0,t),greendiamondsdepictthequan- the “attempted”density blowupregime. The idealgran-
tity1−pT(0,t). Thesolid lineisthepredictionofEq.(35): ularhydrodynamics(IGHD)predictsveryaccuratelythe
an immediate consequence of Eq. (24). The dashed line is a
hydrodynamic profiles, observed in the MD simulations,
numerical solution of the non-ideal hydrodynamic equations
up to times close to the attempted singularity. The den-
(8)-(10) that account for the gradient-dependent transport.
sity field, as measured in the MD simulations, starts to
The inset zooms in at later times.
deviate from the ideal theory at time t 0.9τ. Some-
≃
what surprisingly, the rest of the hydrodynamic fields
indeed shows this growth until t 0.9τ. The quantity continue to show good agreement with the theory un-
≃
1 T(0,t)also saturates at later times, but at a value til even closer to the attempted singularity time τ. In
−
slightlydifferentfrom1 ρ /ρ ,seetheinsetofFig.7. thefollowingSectionwewillseethattheagreementwith
Thispis consistent with t−he pr0essucre deviationfrom ρ at theoryatlatertimesimprovessignificantlywhenthenon-
0
very late times. p ideal hydrodynamic equations (8)-(10), that account for
Wealsopresent,inFig.8,thetimehistoryofthetotal the gradient-dependent transport, are employed.
energy per particle,
E(t) 1 N′ v 2 IV. NON-IDEAL HYDRODYNAMICS
= | | , (36)
N N 2
′ ′ i=1
X Toinvestigatethemechanismofbreakdownoftheideal
asfoundintheMDsimulations,andcompareitwiththe theory, we extended the hydrodynamic theory in two
9
becomesignificant,andtheNIGHDprofilesapproximate
0.04
(a) theMDsimulationresultsmuchbetterthantheidealthe-
ory. As the maximumdensity continues to increase(and
finally approaches the close packing density ρ ), the di-
c
0.03
lute NIGHD descriptionultimately breaksdown. InFig.
3 it occurs at about t/τ 0.98.
m,t) In the second type of h≃ydrodynamic computations we
ρ( 0.02 discarded the viscous and heat diffusion terms but took
into account (moderate) finite density effects. This was
done by adopting, in Eqs. (19) and (20), instead of the
0.01 ideal equation of state and ideal energy loss rate, the
Carnahan-Starling equation of state [27] and a modifi-
cation of the energy loss rate, derived by Jenkins and
0 0.1 0.2 0.3 0.4
m Richman [28] in the spirit of Enskog theory:
πρ
p ρT 1+ g(ρ) ,
0.1 (b) → √3
(cid:20) (cid:21)
Λ Λg(ρ), (37)
→
0.08
where
0.06
m,t) g(ρ)= 1− 372π√ρ3
ρ( 2
0.04 1 πρ
− 2√3
(cid:16) (cid:17)
0.02 is the equilibrium pair correlationfunction of hard disks
at contact. In the dilute limit ρ 0 one obtains g = 1
→
and recovers the ideal equation of state p=ρT.
0
0 0.1 0.2 0.3 0.4 In Fig. 9 we compare four different results for the gas
m
density: the MD simulations, the ideal analytical solu-
tion,thenumericalsolutionofthefirsttype(theNIGHD-
equations),andthenumericalsolutionofthesecondtype.
FIG. 9: (Color online) The density profiles for t/τ = 0.944 Itisclearlyseenthat,forthechoiceofparametersusedin
(a)andt/τ =0.966(b)asobtainedfromtheMDsimulations ourMDsimulations,thenumericalsolutionofthesecond
(the circles), the analytic solution (the dashed-dotted line),
type is not as successful as that of the first type. This
the numerical solution of the dilute NIGHD equations (the
could be expected, as for the times when the maximum
dashed line), and the numerical solution of finite-density hy-
density is still much smaller than the close packing den-
drodynamicequationswiththegradient-dependenttransport
sity,the finite-densitycorrections(whichareofthe order
terms neglected (thesolid line).
of ρ/ρ ) are still negligible. In contrast, the numerical
c
results from the NIGHD equations agree well with the
MDsimulationsandshowthat,astheattempteddensity
separate ways. In the first one we took into account
blowupis approached,the heat conductionandviscosity
thegradient-dependenttransport: theheatdiffusionand
effectscanbecomeimportantwhenthegasdensityisstill
viscosity, but continued to assume the the gas is dilute.
small.
In the second one we neglected the gradient-dependent
transportbuttookintoaccountfinitedensitycorrections.
The hydrodynamic equations were solved numerically in
Lagrangiancoordinates,usinganaccuratevariablemesh V. SUMMARY AND DISCUSSION
and variable time step solver [26].
First,wesolved(theLagrangianformof)Eqs.(8)-(10) Our MD simulations proved the existence of an at-
of non-ideal granular hydrodynamics (NIGHD). These tempted density blowupregimeas describedby anexact
equations account for the viscous and heat diffusion solution of the ideal granular hydrodynamic equations.
terms but still assume a dilute gas. We found the ideal solution to be in remarkable quan-
The numerically obtained NIGHD profiles are pre- titative agreement with the MD simulations over an ex-
sented, together with the MD simulations and the ideal tended time interval, but not too close to the attempted
analytical solution, in Figs. 2-7. As expected, the singularity. As the attempted singularity is approached,
NIGHD profiles coincide with the MD simulations and the exact solution breaks down. A more complete hy-
with the ideal theory for early and intermediate times. drodynamic theory, that accounts for the heat diffusion
At late times, the gradient-dependent transport terms and viscosity, but still assumes a dilute gas, continues
10
to agree with the MD simulations until the gas density by the heat diffusion [3, 18].)
becomes a fraction of the close packing density of disks. Now let us consider Scenario 3 that also assumes a
Let us put the results of this work into a more general long-wavelengthlimit. Asthegastemperaturefallsdown
contextofclusteringinstabilityofafreelycoolinggranu- rapidly because of the inelastic cooling, the flow is de-
larflow. Aswehavealreadynoted,thelocalpropertiesof scribable by the zero pressure (or flow by inertia) ap-
the density blow-up, exhibited by the exact solutions of proximation [10]. Were it not regularized by close pack-
theIGHDequations,arethesameasthelocalproperties ingeffects,suchaflowwoulddevelopafinite-timedensity
of the density blow-up exhibited by a low Mach number blowup (of a different type than the low Mach number
flow of a freely cooling granular gas with the heat dif- flow)[10,22]. Ifthecompressionalheatinginterferesear-
fusion neglected [18]. A low Mach number flow emerges lier than the close packing effects, the pressure becomes
whenthe pressurebalancesetsinonashortertime scale relevant again, and Scenario 3 gives way to the Scenario
than the temperature balance. In this case any local in- 1 [16]. In the opposite case the late time dynamics of
elastic cooling causes a (low Mach number) gas inflow the system is describable by the Burgers equation [11].
into the colder regionso as to increasethe localgas den- Which of the two regimes is realized in a particular set-
sity there and keep the pressure gradient (almost) zero. ting depends on the initial conditions.
The resulting density instability develops on the back- Scenarios1-3donotinvoketheshearmodeinstability,
ground of an (almost) homogeneous gas pressure. This and so they can operate both in one-dimensional, and
isconsistentwiththeMDsimulationspresentedhere,see multi-dimensionalsettings. On the contrary,Scenarios 4
Fig. 6, where the pressure in the vicinity of the density and 5 do invoke the shear mode, and so they are intrin-
maximumhardlychangesuptotimesveryclosetotheat- sically multi-dimensional (and intrinsically non-linear).
temptedsingularitytime. Forbrevitywewillcallthelow Scenario4exploitsthefactthattheunstableshearmode
Mach number flow instability scenario Scenario 1. Sce- may contribute, via a non-linear coupling, to the growth
nario1firstappearedinastrophysicsandplasmaphysics oftheclusteringmode. Obviously,thenonlinearcoupling
in the context of condensation instabilities in gases and istheonly hydrodynamicmechanismofdrivingthe clus-
plasmas that cool by their own radiation [29]. tering mode if the system size is larger than the critical
As many as four additional hydrodynamic scenarios size for the shear mode instability, but smaller than the
of clustering in a freely cooling granular gas have been critical size for the clustering mode instability. Further-
discussed in the literature. We start with the pressure more, numerical analysis, performed in Ref. [6], indi-
instability scenario, or Scenario 2. It was discussed, in cated that the nonlinear coupling plays a dominant role
thecontextofthegranularclustering,byGoldhirschand in the initial density growth also in the case when the
Zanetti [2], although it was also known earlier to the as- system size is comparable to the critical system sizes for
trophysicsandplasmaphysicscommunities,seeRef. [29] the clustering and shear instabilities. (More precisely,
for a review. Scenario 2 is usually presented in the fol- the wave number of the monochromatic test perturba-
lowing way. Let us consider a small local increase in the tion of the transverse velocity in Ref. [6] was within
gas density. This increase causes an increase in the col- the instability regions of both the shear, and the clus-
lisional energy loss. As a result, the gas pressure falls tering modes. However, twice the wave number already
down, a gas inflow develops, causing a further density came out of the instability region.) What happens in
increase, and the process continues. Importantly, Sce- much larger systems is presently under investigation. It
nario 2 assumes that the inelastic cooling time is much turns out that well above the clustering mode instabil-
shorter than the sound travel time. In other words, it is itythresholdthenonlinearcouplingwiththeshearmode
the local temperature balance that sets in rapidly here, does not dominate the density growth (though it does
and the resulting pressure gradient drives the flow on a make the theory more cumbersome). In the channel ge-
relatively slow time scale. ometry,thatweadoptedinthispaperandintheprevious
works [10, 11, 16, 17, 18], the shear mode is suppressed,
As of present, there has been no detailed nonlinear
and the clustering instability develops in its pure and
analysis behind Scenario 2. The (well established) linear
simplest form.
stability theory of the homogeneous cooling state [3, 5]
indicates that Scenarios 1 and 2 operate in two oppo- Now we proceed to Scenario 5 [2] that exploits the
site limits: for sufficiently short and long perturbation fact that the unstable shear mode heats the gas in some
wavelengths,respectively. The physics behind this is the regions. Scenario 5 assumes that this heating can be
following. The characteristic cooling time due to the in- balanced by the inelastic energy loss, rendering (quite
elastic collisions is independent of the lengthscale of the a complicated) steady state. It is furthermore assumed
initial perturbation, whereas the sound travel time scale that this steady state is unstable with respect to small
is proportionaltoit. As aresult, whenallotherparame- perturbations, and it is this instability that causes the
ters are fixed, Scenario 1 corresponds to an intermediate granular clustering. We are unaware of a quantitative
wavelength limit of the clustering instability, while Sce- theorythatwouldsupportScenario5,orofanyquantita-
nario2correspondstothe longwavelengthlimit. (Inthe tivetestofScenario5inMD simulationsorinnumerical
short wavelength limit the homogeneous cooling state of solutions of hydrodynamic equations.
thegasisstable,astheclusteringinstabilityissuppressed To complete the comparisonof the five hydrodynamic
