A nonlinear theory of non-stationary low Mach number channel flows of freely cooling nearly elastic granular gases

0.63 MB·English

by  Baruch Meerson

A nonlinear theory of non-stationary low Mach number channel flows of freely cooling
nearly elastic granular gases
Baruch Meerson, Itzhak Fouxon, and Arkady Vilenkin
Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel
(Dated: February 1, 2008)
We employ hydrodynamicequations to investigate non-stationary channel flows of freely cooling
dilutegasesofhardandsmooth sphereswithnearlyelasticparticlecollisions. Thisworkfocuseson
theregimewherethesoundtraveltimethroughthechannelismuchshorterthanthecharacteristic
cooling time of the gas. As a result, the gas pressure rapidly becomes almost homogeneous, while
8 the typical Mach number of the flow drops well below unity. Eliminating the acoustic modes and
0 employingLagrangiancoordinates,wereducethehydrodynamicequationstoasinglenonlinearand
0 nonlocalequationofareaction-diffusiontype. Thisequationdescribesabroadclassofchannelflows
2 and,inparticular,can follow thedevelopmentoftheclusteringinstability from aweaklyperturbed
homogeneous cooling state to strongly nonlinear states. If the heat diffusion is neglected, the
n
reduced equation becomes exactly soluble, and the solution develops a finite-time density blowup.
a
The blowup has the same local features at singularity as those exhibited by the recently found
J
family of exact solutions of the full set of ideal hydrodynamicequations (Fouxon et al. 2007). The
5
heat diffusion, however, always becomes important near the attempted singularity. It arrests the
1
density blowup and brings about novel inhomogeneous cooling states (ICSs) of the gas, where the

pressurecontinuestodecaywithtime,whilethedensityprofilebecomestime-independent. TheICSs
]
t representexactsolutionsofthefullsetofgranularhydrodynamicequations. Boththedensityprofile
f
o ofanICS,andthecharacteristicrelaxationtimetowardsitaredeterminedbyasingledimensionless
s parameter thatdescribestherelativeroleoftheinelasticenergylossandheatdiffusion. At 1
L L≫
. the intermediate cooling dynamics proceeds as a competition between “holes”: low-density regions
t
a of the gas. This competition resembles Ostwald ripening (only one hole survives at the end), and
m wereportaparticularregimewherethe“holeripening”statisticsexhibitsasimpledynamicscaling
- behavior.
d
n PACSnumbers: 45.70.Qj,47.20.Ky
o
c
[ I. INTRODUCTION Non-stationary flows provide sharp tests to continuum

models of granular flows, especially when the time-
2
dependent solutions of the continuum equations tend to
v Clusteringofmatterisaspectacularexampleofstruc-
5 ture formation in nature. A fascinating example of clus- develop finite-time singularities. Examples are provided
8 by the recently predicted finite-time blowup of the gas
teringisprovidedbygranulargases: gasesofmacroscopic
0 density in freely cooling granular gases: at zero gravity
particles that lose kinetic energy in collisions. Granular
3 [12,17,18](asdescribedbyidealgranularhydrodynamic
gasisalow-densitylimitofgranularflows[1,2]. Thesim-
.
8 equations), and at finite gravity (even in the framework
plest version of the granular gas model assumes a dilute
0 of non-ideal granular hydrodynamic equations) [16].
assembly of identical smooth hard spheres (with diame-
7 We will assume in this paper that particle collisions
ter σ andunit mass)who loseenergyatbinary collisions
0 are almost elastic, the local gas density (that we denote
: in such a way that the normal component of the rela-
v byρ)ismuch smallerthantheclose-packingdensity,and
tive velocity of the colliding particles gets reduced by
i the Knudsen number is very small:
X a constant factor 0 r < 1 (the coefficient of normal
≤
r restitution) upon each collision. Granular gases exhibit 1 r 1, ρσd 1, and l /L 1. (1)
a variouspatternforminginstabilities,includingtheshear- − ≪ ≪ free ≪
ing/clusteringinstabilityofafreelycoolinghomogeneous Here d > 1 is the dimension of space, l is the mean
free
inelastic gas [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. This free path of the particles, and L is the characteristic
instability causes the generation of a macroscopic flow, length scale of the hydrodynamic fields. Under these as-
both solenoidal and potential, and formation of dense sumptions (the secondandthird onesneedto be verified
clusters of particles. a posteriori) the Navier-Stokes hydrodynamics provides
A naturaltheoretical descriptionof macroscopicgran- a quantitatively accurate leading-order theory [1, 2]. It
ular flows is provided by the Navier-Stokes granular hy- wasshown[4,5], byusinghydrodynamicequationsthat,
drodynamics [1, 2]. Although the criteria of its validity for sufficiently large systems, the homogeneous cooling
are quite restrictive, see below, granular hydrodynamics state (HCS) of the granular gas becomes unstable with
hasagreatpredictivepower,sometimesgoingfarbeyond respecttosmallperturbations. Therearetwolinearlyun-
the formal limits of applicability [2]. Recently, granu- stable modes. The shear mode correspondsto the devel-
lar hydrodynamics has been applied to a variety of non- opment of a macroscopic solenoidal flow, while the clus-
stationary flows of granular gases [12, 15, 16, 17, 18]. tering mode correspondsto the developmentof a macro-
2
scopicpotentialflowthatbringsaboutformationofclus- free-flowsingularity,compressionalheating starts to act.
ters of particles. As a result, the gas pressure again becomes important
A consistent nonlinear hydrodynamic theory of the and changes the local blowup properties.
clustering instability has not been available for quite
a long time. Solving strongly nonlinear hydrodynamic
A crucial feature of the finite-time singularity of the
equations is hard (even numerically), and one looks
ideal hydrodynamic equations is that it obeys an iso-
for additional simplifications. Following Refs. [12, 13,
baricscenario: the (finite) gaspressurebecomesuniform
17, 18], we will assume throughout this paper that the
in space in a close vicinity of the developing singularity
macroscopicflow(butnotmicroscopicmotionofthepar-
[18]. Thishintsatthepossibilityofanadditionalsimplifi-
ticles!) is one-dimensional (1d). This assumption is nat-
cationoftheproblem. Indeed,an(almost)homogeneous
ural in the geometry of a narrow channel with perfectly
pressure in a gas implies a low Mach number flow, when
elastic side walls that we adopt here. In a narrow chan-
the inertial terms in the momentum equation are small
nelboththeclusteringmodeinthetransversedirections,
comparedtothepressuregradientterm. Thisregimeap-
and the shear mode are suppressed (see Refs. 12, 13 for
pears when the sound travel time through the system is
detail). As a result, the macroscopic flow can depend
very short compared with other time scales of the prob-
only on the coordinate along the channel and time, and
lem, and one is interested in the dynamics of the system
we can focus on the development of the pure clustering
at the long time scales [21, 22, 23, 24, 25, 26, 27, 28].
mode as it enters a strongly nonlinear regime.
In particular, this regime appears naturally in the linear
Efrati et al. [12] investigated numerically the long-
theory of the clustering instability of the HCS for inter-
wavelength limit of such a quasi-1d clustering instabil-
mediatewavelengthsoftheperturbations,seebelow. Itis
ity. In this limit the inelastic energy loss of the gas is
this (almost) spatially independent pressure regime that
thefastestprocess,sothegaspressurerapidlydropstoa
we will be considering in the present work.
verysmallvalue. Thefurtherdynamicsbecomes(almost)
purely inertial which (almost) brings about a finite-time
blow-up of the velocity gradient and, therefore, of the The remainderof the paper is organizedas follows. In
density [19]. The signatures of this finite-time singular- SectionIIwestartwithafullsetofequationsofgranular
ity wereindeedobservedinthe numericalsolutionofthe hydrodynamicforadilutegranularflowinachanneland
hydrodynamic equations [12] until the growing gas den- reduce them, for sufficiently short channels, to the low
sity became so high that the numerical scheme lost ac- Mach number flow equations. In Section III we employ
curacy. The numerical results of Ref. [12] were tested Lagrangiancoordinateswhichenableustoexactlyreduce
in molecular dynamics (MD) simulations [13]. The MD the low Mach number flow equations to a single nonlin-
simulationssupportedthe free-flowblow-upscenarioun- ear and nonlocal equation, of a reaction-diffusion type,
til the time when the gas density approachedthe hexag- for the square root of the inverse gas density. The new
onal close-packing value, and the further density growth equationis tested in Section IV on two simple problems:
stopped. the HCS and the linear theory of clustering instability
Recently, Fouxon et al. [17, 18] analyzed, analytically in short channels. In Section V we show that, when the
and numerically, the one dimensional flow in the frame- heat diffusion is neglected, the new equation becomes
work of equations of ideal hydrodynamics (that is, ne- exactly soluble, and the solution develops a finite-time
glecting the heat diffusion and viscosity effects). They densityblowupwiththesameuniversalfeaturesatsingu-
derived a family of exact solutions to these equations, larity as those exhibited by the family of exact solutions
with and without shocks, for which an initially smooth of the full set of ideal granular hydrodynamic equations
flow develops a finite-time density blowup. Close to the [17,18]. SectionVI presentsananalyticalandnumerical
blow-up time t , the maximum density exhibits a power analysisthatshowsthattheheatdiffusionterm,nomat-
c
law behavior (t t) 2. The velocity gradient blows ter how small in the beginning, becomes important near
c −
∼ −
up as (t t) 1, whereas the velocity itself remains the attempted density blowup. As a result, the density
c −
∼ − −
continuousanddevelops a cusp, ratherthana shock dis- blowup is arrested, and a novel, inhomogeneous cooling
continuity, at the singularity. The gas temperature van- state (ICS) of the gas emerges, with a time-independent
ishes at the singularity, but the pressure remains finite. inhomogeneous density profile. Importantly, the ICSs
Extensive numericalsimulations with the ideal hydrody- represent exact solutions of hydrodynamic equations. A
namicequationsshowedthatthesingularityexhibitedby limitingformofthenovelcoolingstateiswhatwecallthe
the exactsolutionsis universal,as itdevelopsforgeneric “hole”, and we investigate its properties and the relax-
initial conditions. Very recently, the existence of the at- ationdynamics towardsit. For sufficiently long channels
temptedblowupregimehasbeenprovedinmoleculardy- (other parameters being fixed) the cooling dynamics of
namic simulations of a gas of nearly elastically colliding the systemtakesthe formofacompetitionbetween,and
harddisksinachannelgeometry[20]. TheresultsofRefs. “ripening” of, holes. Therefore, in Section VII we inves-
[17, 18] also imply that, for long wavelength initial con- tigatethedynamicsandstatisticsofthiscompetition. In
ditions, the free flow regime may not hold all the way to SectionVIIIwesummarizeourresultsandputtheminto
the density blowup [17, 18]. Very close to the attempted a perspective.
3
II. GRANULAR HYDRODYNAMICS AND A The second characteristic length scale is the heat dif-
LOW MACH NUMBER FLOW fusion length
For flows depending on a single spatial coordinate x κ p1/2t 1/2 κ1/2
l 0 0 c 0
and time t the granular hydrodynamic equations can be d ∼ ρ3/2 ! ∼ Λ1/2ρ0
written as follows: 0
which, up to a numerical pre-factor, coincides with the
∂ρ ∂(ρv)
+ =0, (2) critical length
∂t ∂x
∂v ∂v ∂(ρT) ∂ ∂v
ρ +v = +ν √T , (3) 2κ
∂t ∂x − ∂x 0∂x ∂x l = 0 , (7)
∂(cid:18)T ∂T (cid:19) ∂v (cid:18) (cid:19) cr sΛρ20
+v = (γ 1)T ΛρT3/2
∂t ∂x − − ∂x − predictedbythelineartheoryoftheclusteringinstability.
+κρ0∂∂x √T∂∂Tx + ν0(γ−ρ1)√T ∂∂xv 2 . (4) wTeheharavteioallrse/alddyisasosfuomrdedera(κst0rΛo)n−g1i/n2e∼qu(a1li−tyr12)−r12/2. A1s,
(cid:18) (cid:19) (cid:18) (cid:19) this ratio is very large: l /l 1. Throughou−t the≪rest
s d
≫
Here γ is the adiabatic index of the gas (γ = 2 and 5/3 of the paper we will also assume that the channel length
for d = 2 and d = 3, respectively), Λ = 2π(d 1)/2(1 Lismuchshorterthanthesoundtraveldistancel . This
r2)σd−1/[dΓ(d/2)] (see e.g. [8]), Γ(...) is the−gamma−- hierarchy of length scales brings about a reducedsset of
function, and d 2 is the dimension of space, so that equations,inmuchthe samewayasinhydrodynamicsof
≥
d = 2 corresponds to disks, and d = 3 to hard spheres. optically thin gases and plasmas that cool by their own
Furthermore, ν = (2σ√π) 1 and κ = 4ν in 2D, and radiation [21, 23, 24, 25, 28]. Note that the length scale
0 − 0 0
ν0 = 5(3σ2√π)−1 and κ0 = 15ν0/8 in 3D [1]. Equa- separation L ls is equivalent to a time scale separa-
≪
tions(2)-(4)differfromthehydrodynamicequationsfora tion: the soundtraveltime throughthe channel,L/c ,is
s
dilutegasofelastically collidingspheresonlybythepres- much shorter than the characteristic cooling time t . As
c
ence of the inelastic loss term ΛρT3/2 which is propor- aresult,soundwavesrapidlymakethepressure(almost)
tionaltotheaverageenergylos−spercollision, (1 r2)T, homogeneous throughout the channel. The subsequent
∼ −
and to the collision rate, ρT1/2. slowerevolutionofthegasproceedsonthebackgroundof
∼
It will be convenient for our purposes to rewrite analmosthomogeneous(but in generaltime-dependent)
Eqs.(2)-(4) interms ofthe pressurep=ρT,ratherthan gaspressure,while typicalMachnumbersof the flow are
the temperature. The energy equation (4) becomes much less than unity. In a more formal language, this
reductionof the hydrodynamic equations correspondsto
∂p ∂p ∂v
+v = γp Λρ1/2p3/2 elimination of acoustic modes.
∂t ∂x − ∂x − Before we perform the reduction procedure, let us in-
∂ p ∂ p p ∂v 2 troduce rescaledvariables. We will measure the distance
+κ +ν (γ 1) .(5)
0∂x ρ ∂x ρ 0 − ρ ∂x along the channel in the units of lcr, rescale time by tc,
(cid:20)r (cid:18) (cid:19)(cid:21) r (cid:18) (cid:19) andmeasurethegasdensity,pressureandvelocityinthe
A set of hydrodynamic equations can be simplified if units ofρ , p andl /t , respectively. Keeping the orig-
0 0 cr c
there is a time scale separationor, equivalently, a length inal notation for the rescaled variables, we observe that
scale separation, in the problem. For a freely cooling Eq. (2) does not change, while Eqs. (3) and (5) become
granulargas,abasictime scaleis thecharacteristiccool-
ing time ∂v ∂v ∂p ∂ p∂v
ε ρ +v = +ε , (8)
1 2
∂t ∂x −∂x ∂x ρ∂x
2 (cid:18) (cid:19) (cid:18)r (cid:19)
tc = Λρ1/2p1/2 , (6) ∂p +v∂p = γp∂v 2ρ1/2p3/2
0 0 ∂t ∂x − ∂x −
where ρ is the average gas density (the total gas mass ∂ p ∂ p p ∂v 2
0 + +ε (γ 1) ,(9)
divided by the volume of the channel), and p0 is a char- ∂x ρ ∂x ρ 2 − ρ ∂x
(cid:20)r (cid:18) (cid:19)(cid:21) r (cid:18) (cid:19)
acteristic value of the initial pressure. There are two
characteristic length scales related to t . The first is the where ε = κ Λ/2, and ε = ν Λ/2, and ε ε
sound travel distance c 1 r2 1 1. W0 e will limit2 ourse0lves to the 1ze∼roth2o∼r-
− ≪
der approximation with respect to this small parameter
γ√2
and send ε and ε to zero. The continuity equation (2)
l = c t , 1 2
s s c
Λρ0 ∼ does not change. The momentum equation (8) becomes
∂p/∂x = 0, therefore p = p(t) is independent of x. The
which is of the order of the distance a sound wave with
energy equation becomes
speed c = (γp /ρ )1/2 travels during the time t . The
s 0 0 c
qinuaRnetfist.y[1ls7,is18th].e same as the length scale l introduced p˙(t) = γp(t)∂v 2ρ1/2p(t)3/2
− ∂x −
4
+ p(t)3/2 ∂ 1 ∂ 1 . (10) boundary conditions (BCs) one can always achieve this
∂x √ρ ∂x ρ byexploitingtheGalilianinvarianceofthehydrodynamic
(cid:20) (cid:18) (cid:19)(cid:21)
equations to get rid of the center-of-mass motion. This
The rescaled length of the channel is sets v(x = 0,t) = 0, where x = 0 is the center-of-mass
coordinate. For the no-flux BC (impenetrable walls), a
= L = (√π/2)(1−r2)1/2ρ0σL in 2d, (11) natural choice of x = 0 is at one of the walls, where the
L lcr (cid:26) 16π/75(1−r2)1/2ρ0σ2L in 3d. gas velocity is again zero. Then a convenient choice of
the Lagrangianmass coordinate is
p
Note that, in the rescaled variables, the rescaled length
x
of the channel coincides with the rescaled total mass
L m(x,t)= ρ(x′,t)dx′, (13)
of the gas, Lρ(x,t)dx.
0 Z0
Togetanexplicitexpressionforp˙weintegrateEq.(10)
R which is simply the (rescaled) mass content between the
overthewholechannel. Assuming eitherperiodic,orno-
Eulerian points 0 and x. The inverse transformation is
fluxboundaryconditions(BCs)atthechannelendsx=0
and x= , we obtain m dm
L ′
x(m,t)= . (14)
ρ(m,t)
p˙(t) = 2 ρ1/2(x,t) , (12) Z0 ′
p(t)3/2 − x In the Lagrangiancoordinates Eqs. (2) and (10) become
D E
where we have introduced the spatial average ∂ 1 ∂v
= , (15)
∂t ρ ∂m
... = 1 L(...)dx. (cid:18) (cid:19) ∂v
h ix p˙ = γpρ 2p3/2ρ1/2
LZ0 − ∂m −
∂ ∂ 1
ForthelowMachnumberflow,Eq.(12)describes,inthe + p3/2ρ √ρ . (16)
leading order, the global energy balance of the gas, see ∂m ∂mρ
(cid:18) (cid:19)
Section VI C below. Equations (2), (10) and (12) for
As the total rescaled mass of the gas is equal to the
ρ(x,t), v(x,t) and p(t) make a complete set of reduced
rescaled channel length , we define the spatial average
but fully nonlinear equations for the low Mach number L
in the Lagrangiancoordinate as
flow of a freely cooling granular gas in a channel geom-
etry. As is usually the case for low Mach number flows, 1 L
theviscoustermsdroppedfromthereducedformulation, ... = (...)dm,
h i
while the heat diffusion term remains. LZ0
Therescaledlength/massofthesystem ,seeEq.(11), and rewrite Eq. (12) as
L
is determined by the relative role of the inelastic energy
lossandheatdiffusion. Aswewillseeshortly, controls p˙(t) 1
= 2 . (17)
the main properties of the cooling dynamics.LFor com- p(t)3/2 − ρ1/2(m,t)
(cid:28) (cid:29)
parison,thecharacteristicinitialpressurep onlysetsthe
0 It is convenient to introduce a new rescaled variable
time scaleforthe dynamics. Tofacilitatefuture compar-
w(m,t)=ρ 1/2(m,t) and a new rescaled time
isons of the theory with MD simulations, we rewrite the −
parameter in a slightly different form: 1 t
L τ = p1/2(t)dt . (18)
′ ′
γ
√π(1−r2)Npσ in 2d, Z0
= 2Ly Then, by eliminating ∂ v and p˙ from Eqs. (15)-(17), we
L  √16π(1−r2)Npσ2 in 3d. canreducetheseequatimonstoasingleintegro-differential
 √75LyLz
equation of a reaction-diffusion type:
Here N is the total number of particles in the channel,
p ∂w ∂2w
and Ly and Lz are the transverse channel dimensions. w = w+w2 w + . (19)
∂τ − h i ∂m2
Equation(19) describes a broadclass of slow 1d flows in
III. LAGRANGIAN DESCRIPTION AND
freelycoolingnearlyelasticgranulargases. Inparticular,
NONLOCAL REACTION-DIFFUSION
this equation encodes the development of the clustering
EQUATION
instability: from a weakly perturbed HCS (after a brief
acoustic transient) all the way to the strongly nonlinear
Remarkably, it is possible to bring the three equa-
stage. Indeed,letusrewriteEq.(17)intermsofthenew
tions(2),(10)and(12)toasinglenonlocalequationofa
variable w(m,τ) and new time τ:
reaction-diffusiontype. LetusfirstintroduceLagrangian
mass coordinates [29]. It is convenient to choose a ref- 1 dp
= 2γ w(m,τ) . (20)
erence frame so that v(x = 0,t) = 0. For the periodic p(τ)dτ − h i
5
Once Eq. (19) for w(m,τ) is solved, we can calculate perturbationδρ(m,τ)= 2δw(m,τ).]Onecanrepresent
−
the pressure p(τ) from Eq. (20) and then return to the δw(m,t) as a linear superposition of sines and cosines of
(rescaled) physical time t using Eq. (18): km with different (rescaled) wave numbers k. This fact,
in conjunction with the BCs at the ends of the channel,
τ dτ′ guarantees that δw(m,τ) =0. Then Eq. (19) yields
t=γ . (21) h i
p1/2(τ )
Z0 ′
∂ ∂2
Furthermore, using Eq. (15) and the condition v(m = δw(m,τ)=δw(m,τ)+ δw(m,τ). (26)
∂τ ∂m2
0,t) = 0, we can find the gas velocity: v(m,t) =
0m∂tw2(m′,t)dm′. Finally, we can return to the For a single mode perturbation with wave number k we
Eulerian coordinate by using Eq. (14): x(m,t) = obtain
Rmw2(m,t)dm.
0 ′ ′
Notably, equation (19) is parameter-free: the only pa- δw(m,τ)=δw(m,0)eΓˆkτ (27)
R
rameter entering the problem [except possible parame-
ters introduced by the initial condition w(m,0)] is the with the growth/damping rate
rescaledsystemlength/mass . Conservationofthetotal
L
massofthegasinthechannelappearsintheLagrangian Γˆ =1 k2. (28)
k
formulation as the conservation law −
For k < k = 1 (correspondingly, k > k = 1) Eqs. (27)
hw2(m,τ)i=1, (22) and (28) ∗describe an exponential grow∗th (correspond-
ingly,decay)ofa smallsingle-modeperturbationintime
easily verifiable from Eq. (19).
τ. Recalling that we rescaled the coordinate to the crit-
ical length l , provided by the complete (unreduced)
cr
linear theory, we immediately notice that Eq. (28) cor-
IV. SIMPLE TESTS: HCS AND LINEAR
rectly predicts the instability threshold. To go back to
THEORY OF CLUSTERING INSTABILITY
the physical time t we substitute, in the leading order,
Haff’slaw(25)intoEq.(18)andobtain,afterelementary
As a first test of Eqs. (19) and (20), let us consider
integration,
a HCS. Here at t = 0 we have (in the physical units)
ρ(m,t = 0) = ρ = const, T(m,t = 0) = T = const
0 0 1
and v(m,t = 0) = 0 and, therefore, w(m,t = 0) = 1 τ = ln(1+t) . (29)
γ
and p(t = 0) = p = ρ T . As the gas density remains
0 0 0
constant in space at t>0, we can rewrite Eq. (19) as
Plugging it into Eq. (27), we obtain an algebraic growth
of the small perturbations in the physical time:
dw(τ)
= 1+w2(τ). (23)
dτ −
δw(m,t)=δw(m,0) (1+t)Γˆk/γ . (30)
The solution of this equation with the initial condition
w(0)=1isofcoursew(τ)=1: the gasremainsspatially ThegrowthexponentΓ=Γˆ /γ,withΓˆ fromEq.(28)
k k
homogeneous. Now we use Eq. (17) and obtain
coincides with that obtained from the complete linear
stability analysis [5], if we assume there kl 1 (in the
p˙(t) s ≫
= 1, (24) physical units) and consider the clustering mode, rather
p(t)3/2 −
than the two decaying acoustic modes. Figure 1 shows
this comparison in a graphic form. At kl < 1 the iso-
which yields, in the rescaled variables, Haff’s law for the s
baric growth rate underestimates the true g∼rowth rate,
gas pressure:
but in the region of kl 1 excellent agreement is ob-
s
≫
1 served. The comparison with the complete linear stabil-
p(t)= . (25)
(1+t)2 ity analysis is instructive for two more reasons. First, as
was observed by McNamara [5], for kl 1 the pres-
s
≫
The next test of Eq. (19) is the linear stability anal- sure perturbations of the clustering mode vanish in the
ysis of a HCS. While the reduced Eq. (19) is not sup- leading order in 1/(kls). That is, the linear density and
posed to capture the evolution of small perturbations temperature perturbations grow on the background of
with an arbitrary polarization, it must reproduce cor- an (almost) constant pressure. Second, the viscosity ef-
rectly the evolution of the clustering mode in the limit fectsdonotaffectthegrowthexponentinthisregime[5].
when the perturbation wavelengths are small compared As our reduced formalism shows, the last two properties
with the sound travel distance l . Let us show it to be persist, for the low Mach number flow, in the nonlinear
s
indeed the case. We look for the solution of rescaled regime as well.
Eq. (19) in the form w(m,τ) = 1 + δw(m,τ), where Having successfully tested our reduced model in these
δw(m,τ) 1. [Correspondingly, the rescaled density two simple cases, we now consider nonlinear evolution.
| | ≪
6
1 a behaviorofthe(rescaled)gasdensitynearthesingularity
0.8 is described by the following equation:
0.6
Gk0.4 ρ(m,τ)= τ τ + 1 d2w˜ (m )(m m )2 −2 , (32)
0.2 c− 2dm2 0 − 0
(cid:20) (cid:21)
0 where the time of singularity τ = w˜(m ). The singu-
c 0
20 40 60 80 100
larity structure, as described by Eq. (32), coincides with
kl
s that exhibited by a family of exact solutions of the full
1 set of ideal hydrodynamic equations [that is, Eqs. (2)-
b
0.8 (4) without the viscous and heat diffusion terms], re-
0.6 ported in Ref. [17, 18]. At τ = τc the density blows
Gk0.4 up as ρ(m,τc) (m m0)−4. Going back to the Eule-
∼ −
rian coordinate, we obtain a finite-mass density blowup
0.2
ρ(x,τ ) x x 4/5, where x is the Eulerian coordi-
0 c ∼ | − 0|− 0
nate of the singularity. We refer the reader to Ref. [18]
-1-0.5 0 0.5 1 1.5 2
foradetailedanalysisofthe structureofthissingularity,
log Hkl L
10 s asobservedin the gasdensity, temperature andvelocity.
Notably, the pressure field does not have any singularity
FIG. 1: (Color online.) The growth exponent Γ of theclus-
k inthe exactsolutions [17, 18], andis approximatelycon-
teringmodeversustherescaledwavenumberkl (a)andver-
s
stant in a narrow region around the density singularity.
suslog(kl )(b)aspredictedfromthecompletelinearstability
s
Thatis, the density blowup,as featured by the exactso-
analysis of a HCS [5] (the thick black line) and from the re-
duced equation (19) (the thin red line). The physical (not lutions of ideal granular hydrodynamics [17, 18], locally
rescaled) units are used, and the parameters are γ = 2 and obeys an isobaric scenario, as was noticed in Ref. [18].
k l =100. At kl <1 the isobaric growth rate (19) under- Thisprovidesthereasonwhythesametypeofsingularity
cr s s
estimatestheactual∼growthrate,butintheregionofkls 1 appears in our reduced low Mach number theory.
≫
excellent agreement between thetwo results is observed. NowwepresentacompletesolutionofEq.(31). First,
weobtainaclosedevolutionequationforthe(necessarily
positive) quantity χ(τ) = w(m,τ) by integrating the
V. NEGLECTING HEAT DIFFUSION CAUSES h i
both sides of Eq. (31) over m from 0 to :
A DENSITY BLOWUP L
dχ(τ)
= 1+χ2. (33)
Aswementionedearlier,theonlygoverningparameter dτ −
in Eq. (19), except parameters introduced by the initial We consider the solutionofthis equationwiththe initial
condition, is the rescaled system length/mass . In the condition
L
limit of 1, and for a sufficiently large-scale initial
conditionL, o≫ne can drop the diffusion term in Eq. (19). χ0 = w(m,0) = ρ(m,0)−1/2 1.
h i ≤
This approximation is valid as long as the solution re-
D E
mains large-scale. At the level of linear stability analy- The solution can be written as
sis this (intermediate-wavelength)approximationis fully χ tanh(τ)
0
jinutsetrifieestde.d HinerteheEnqo.n(l2in8e)abredceovmeleospmΓˆken≃t o1f,thaendgroonweinigs χ(τ)= 1−−χ0tanh(τ), (34)
perturbations. With the diffusion term dropped we ob- Now we can rewrite Eq. (31) as
tain ∂w
χ(τ)w(m,τ) = 1, (35)
∂w ∂τ − −
= 1+w w . (31)
∂τ − h i where χ(τ) is given by Eq. (34). Equation (35) is easily
soluble:
This nonlinear integro-differential evolution equation is
exactly soluble for any initial data w(m,0). The com- w(m,0)+χ [cosh(τ) 1] sinh(τ)
0
w(m,τ)= − − . (36)
plete solution is presented below. The main result here cosh(τ) χ sinh(τ)
0
is that, for any inhomogeneous initial condition, the so- −
The presence of the factor χ [cosh(τ) 1] sinh(τ) in
lution of Eq. (31) develops a zero w (hence an infinite 0
− −
the numeratorofEq.(36) causes,for any(non-constant)
density) in a finite time. Let us first discuss the proper-
initial data w(m,0), a singularity w 0 ina finite time.
ties of the solution in a close vicinity of the singularity
→
ThesingularityoccursattheLagrangianpointm where
w 0. In the leading order we can neglect the inte- 0
→ the function w(m,0) has its minimum, at time
gral term in Eq. (31) and obtain ∂w/∂τ = 1, so that
−
w(m,τ) = w˜(m) τ, where w˜(m) is a smooth function.
− ∆w+ ∆w2+1 χ2 1/2
The singularity occurs in the Lagrangian point m0 that τc =ln − 0 (37)
correspondstotheminimumofw˜(m). Theleadingorder " (cid:0) 1−χ0 (cid:1) #
7
where ∆w w(m0,0) χ0. Note that ∆w = 2.5 a
≡ −
(1/L) 0L[w(m0,0)−w(m,τ)]dm ≤0. 2
Now we compute the (rescaled) pressure p(τ) from 1.5
Eq. (2R0) [note that the right hand side is simply χ(τ) 1€€€€€Ρ 1
given by Eq. (34)], 0.5
0
p(τ)=(coshτ χ sinhτ)2γ , (38)
0 0 20 40 60 80 100
−
m
and use this result in Eq. (21) for the rescaled physical
time: 100
b
τ dτ 80
′
t=γ . (39)
(coshτ χ sinhτ )γ 60
Z0 ′− 0 ′ x
40
For γ =2 (a 2D gas of disks) this integralis elementary,
20
and the result is
0
2 0 20 40 60 80 100
t= . (40)
cothτ χ m
0
−
Now we can express τ through t, 2.5 c
2
2
τ =arccoth +χ0 (41) 1.5
(cid:18)t (cid:19) 1€€€€€Ρ 1
and rewrite Eqs. (36) and (38) (for γ =2) as 0.5
0
1
w(m,t) = [w(m,0) χ ] 4+4χ t (1 χ2)t 0 20 40 60 80 100
2{ − 0 0 − − 0
x
+ 2χ (1 χ2)t .q (42)
0− − 0 }
and FIG. 2: The density history of a freely cooling gas of inelas-
ticharddisksina2Dchannelinthezero-heat-diffusionlimit.
t2 −2 Therescaledinitialdensityρ(m,0)=[1+0.1cos(2πm/ )]−1.
p= 1+χ t 1 χ2 . (43) L
0 − − 0 4 The rescaled system length/mass = 100. Panel a: the
(cid:20) (cid:21) rescaledinversedensityofthegas,1L/ρ,versustheLagrangian
(cid:0) (cid:1)
So, the solution for γ = 2 is surprisingly simple. We mass coordinate m at times τ = 0, 1.5, 2.5 and the time of
remind that, in view of the chosen rescaling, the initial singularity τ 2.8755. Panel b: the rescaled Eulerian coor-
c
≃
condition w(m,0) must obey w2(m,0) = 1. To return dinate x versus m at the same times. Panel c: 1/ρ versus x
to the HCS and Haff’s law in Eqs. (42) and (43) one atthesametimes. Thesequenceofcurvesisself-explanatory.
(cid:10) (cid:11)
shouldputtherew(m,0)=χ =1. Equation(43)shows
0
thatHaff’slawisanupperboundforthethermalenergy
loss rate: any deviation from homogeneity brings about
χ <1 and a slower thermal energy decay.
0
Let us note that the solution (34) for χ(τ) vanishes
at τ = (1/2) ln[(1+χ )/(1 χ )] and becomes nega-
0 0
tive∗at larger τ. This is in ap−parent contradiction with
the positivity of w that dictates χ(τ) = w(m,τ) 0.
h i ≥
The contradictionis resolvedby noting that τ is always where E(...) is the complete elliptic integral of the sec-
greaterthanthesingularitytimeτc,beyondwh∗ichtheso- ond kind, see e.g. [30]. Figure 2a shows, at different
lution does not apply. [To see that τc τ one can use, times,therescaledinversedensity1/ρ(m,τ),asobtained
inEq.(37),that ∆w+ ∆w2+1 χ2 ≤1/2∗ 1 χ2 1/2 from Eq. (36), for δ = 0.1 and = 100. Figure 2b de-
for any ∆w 0.] Similarly, the p−ress0ure as≤pred−icte0d by picts,atthesametimes,therescLaledEuleriancoordinate
Eq.(38)orE≤q.(43)wo(cid:0)uldstartincrea(cid:1)singat(cid:0)someti(cid:1)me. x = 0mw2(m′,τ)dm′ versus the Lagrangian coordinate
Atphysically meaningfultimes τ <τ ,however,we have m. Figure 2c shows the rescaled inverse density in the
c R
χ(τ) > 0, and the pressure always decreases in accord rescaled Eulerian coordinates and illustrates the emer-
with Eq. (20). gence of the cusp density singularity at x = /2. The
L
As a simple illustration of our solution (36), let inversedensitybehaveslike(m m0)4 atsmallm m0in
us chose the following initial condition: w(m,0) = the Lagrangiancoordinate, and−like (x x0)4/5 a−t small
−
[1+δcos(2πm/ )]1/2, 0<δ <1. In this case x−x0 in the Eulerian coordinate. This simple example
L is instructive as, for δ 1, this initial condition corre-
≪
1 2δ 2δ sponds to a small single-mode density perturbation, so
χ = √1 δ E +√δ+1E ,
π − δ 1 δ+1 the initial evolution is describable by the linear theory.
(cid:20) (cid:18) − (cid:19) (cid:18) (cid:19)(cid:21)
8
VI. HEAT DIFFUSION ARRESTS THE that, in virtue of Eq. (46), is obeyed automatically for
DENSITY BLOWUP the periodic or no-flux BCs.
Equation(46) has appearedin numerous applications,
A centralresult of this work is in that, no matter how and its solutions are well known. It is convenient to in-
smallinitially,the heatdiffusionterminEq.(19)arrests terpretf as a coordinate of a Newtonian particle of unit
the density blowup. An emerging balance of the inelas- mass, moving in a potential U(f) = f3/3 f2/2. The
−
ticcoolingandheatdiffusionleadstoexistenceofsteady “total energy” E is conserved:
statesolutionsofEq.(19). Thesesolutionsdescribenovel
cooling states of the granular gas, where the (inhomoge- 1 df 2 f3 f2
E = + . (49)
neous)densityprofileistime-independent,whilethe(ho- 2 dm 3 − 2
(cid:18) (cid:19)
mogeneous) pressure continues to decay with time. We
Forthebounded(spatiallyoscillating)solutions, 1/6
found that, in our rescaled variables, the density profile
− ≤
E 0, we can write
of the novel cooling state is uniquely defined by the pa-
≤
rameter . For sufficiently large values of the rescaled
L f3 f2 (f a)(f b)(f c)
length/mass, 1, the maximum gas density of the E = − − − , (50)
L ≫ 3 − 2 − 3
novel cooling state is exponentially large in . In the
L
low Mach number theory, considered in this work, the where a > b > c are the real roots of the cubic poly-
novel cooling states represent global attractors, as they nomial. Then the bounded solutions of Eq. (46) can be
develop for any inhomogeneous initial conditions. Fi- written as
nally, the novel cooling states represent exact solutions
of the complete, unreduced set of hydrodynamic equa- a c
f(m)=c+(a c)dn2 − m, s , (51)
tions (2)-(4). − r 6 !
where
A. Steady state density profiles
a b
s= − , (52)
a c
Steady-statesolutionsofEq.(19)aredescribedbythe
−
equation
anddnisoneoftheJacobiellipticfunctions,seee.g. [30].
TherearetwolimitswhenEq.(51)simplifies. Inthelimit
d2w
=w w w2. (44) ofE = 1/6+δE,0<δE 1,thesolution,f(m)=1+
dm2 −h i − ≪
√2δEcosm,correspondstoasmall-amplitudesinusoidal
Notice that, although obtained from our reduced, low modulationoftheHCSw(m)=1. InthelimitofE 0,
→
Mach number theory, Eq. (44) also follows from the full we have a=3/2 and b=c=0, so that
set of hydrodynamic equations (2)-(4), if one assumes a
3 m 3 m
homogeneouspressureand zerofluid velocity,andtrans- f(m,E 0)= dn2 ,1 = cosh−2 , (53)
forms to the Lagrangiancoordinates. → 2 2 2 2
(cid:16) (cid:17) (cid:16) (cid:17)
Equation (44) is defined on the interval 0 m ,
Using Eqs. (47) and (51), we rewrite the steady state
≤ ≤ L
at the ends of which we demand either periodic, or no-
solutions in terms of w(m):
flux (zero first derivative) BCs. The solutions we are
interestedinmustobeytheconservationlaw(22). Toget
c+(a c)dn2 a cm,s
rid of the (a priori unknown) factor w , we introduce a − −6
new variable h i w(m)= c+(a (cid:16)cq) E(s) (cid:17), (54)
− K(s)
f(m)= w w(m) (45) r
h i
where K(s) is the complete elliptic integral of the first
and obtain
kind. The lagrangian spatial period, or wavelength, of
d2f the solution (54) is
f +f2 =0. (46)
dm2 −
24
Π= K(s). (55)
Once f is found, one can restore w via a c
r −
f Inthelimitofs 0(orE 1/6),thewavelength(55)
w = . (47) → →−
reaches its minimum value 2π. If the rescaled channel
f
h i length is less than 2π (for the periodic BCs), or less
L
The conservation law (22) epnforces a normalization con- than π (for the no-flux BCs), the only possible steady
dition state is the constantdensity state w(m)=1 correspond-
ingtoHaff’slaw. Thisresultisinfullagreementwiththe
f2 = f (48) linear stability analysis of Eq. (19), see Eq. (28). When
h i h i
9

exceeds2π (fortheperiodicBCs),orπ (fortheno-flux
LBCs), the HCS bifurcates into an inhomogeneous steady 4 =0 4 =2
state (54). In general, the rescaled channel length/mass w w
2 2
mustbeequal,byvirtueoftheBCs,toanintegernum-
L
ber of Π (for the periodic BCs), or to an integer number 0 0
0 10 20 30 40 50 0 10 20 30 40 50
of Π/2 (for the no-flux BCs). For sufficiently large value
m m
of , therefore, a whole family of steady state density
proLfiles exists. Which of the steady state solutions is 4 =4 4 =72
w w
selected by the dynamics of Eq. (19)? 2 2
0 0
0 10 20 30 40 50 0 10 20 30 40 50
B. Selected steady-state solutions: the m m
9
inhomogeneous cooling states 2 10
w =72
9
/1 10
We performed extensive numerical simulations with 1
Eq. (19), using a specially developed numerical scheme 0
describedinAppendixA.Bothperiodic,andno-fluxBCs 0 10 20 30 40 50
were used. We observed that, when 0< <2π (for the m
L
periodic BCs), or0< <π (forthe no-flux conditions), FIG.3: Numericalw-profilesattimesτ =0,2,4and72,and
the HCS appears, as eLxpected. When exceeds 2π (for 1/w-profileattimeτ =72,for =50whenstartingfromthe
the periodic BCs), a weakly inhomogenLeoussteady state initialcondition(57). ThetwopLanelsforτ =72alsoshow,by
circles,thesingleholeasymptotes(60)and(61),respectively.
densityprofilesetsin. As increasesfurther,theweakly
L
inhomogeneousstatesdevelopsintoastronglyinhomoge-
neous states. The simulations showed that the rescaled
length/mass of the gas, , uniquely selects the emerging 1.4 a
L
steady state density profile, while the initial w-profile 1.2
does not play any role in the selection. For a given 1
L w
the dynamics always selects, out of the family of steady 0.8
state solutions (54), the one with the maximum possible
0.6
wavelength Π:
0 1 2 3 4 5 6 7
Π for the periodic BCs,
= (56) m
L Π/2 for the no-flux BCs.
(cid:26)
4 b
3.5
Snapshots from a typical simulation (one of many that
3
we performed) for the periodic BCs are shown in Fig. 3. 2.5
Ρ
The initial condition is this example was 2
1.5
1
w2(m,0) = 1 0.1 cos(2πm/ ) 0.15 sin(2πm/ ) 0.5
− L − L
+ 0.2 cos(4πm/ ) 0.05sin(4πm/ ). (57) 0 1 2 3 4 5 6 7
L − L
x
The rescaledsystemlength/mass =50wassufficiently
L
large to fit in steady state solutions with several oscil-
FIG. 4: The inhomogeneous cooling state for = 7.025.
lations. Nevertheless, the dynamics selected the solu- L
Panel a: the Lagrangian steady state solution w(m) as pre-
tion with the spatialperiod equal to the rescaledsystem
dicted by Eq. (54). Panel b: the rescaled steady state gas
length . density ρ versusthe rescaled Eulerian coordinate x.
L
Figures 4 - 6 depict our analytical solutions (54) in
the Lagrangian coordinate, and the corresponding den-
sity profiles in the Eulerian coordinates, for three dif-
ferent values of the parameter . Here we assumed the One can see that, as the parameter increases, the
L L
periodic BCs and (arbitrarily) chose the position of the maximumgasdensityintheclustergrowsveryfast[note
minimum of w(m) to be in the middle of the channel. that Fig. 6b shows the density in logarithmic scale]. Let
The maximum (rescaled) gas density versus the usconsidertheasymptoticformofthesolutionat 1
L≫
rescaled channel length , predicted by Eqs. (54) and in some detail. The density maximum in this case is
L
(55), is shown in Fig. 7. This dependence can serve as exponentially large [31]. This is due to the behavior of
a bifurcation diagram of the system. One observes, at the s 1 asymptotics of the steady-state solution, see
→
> 2π, a supercritical bifurcation from the HCS to an Eq. (53). In this case the “energy” E is very small, and
L
ICS. can be expressed through the rescaled system length as
10
1.75 a 12 a
1.5 10
1.25
8
w 1 max 6
0.75 Ρ
0.5 4
0.25 2
0
0 2 4 6 8 2 4 6 8 10121416
m rescaledlength
30 1.2
b b
25
1.15
20
x
Ρ 15 ma 1.1
10 Ρ
5 1.05
0
1
0 2 4 6 8
6.1 6.2 6.3 6.4 6.5
x
rescaledlength
FIG. 5: Same as in Fig. 4, but for =8.886.
L FIG.7: Thebifurcationdiagramofthefreelycoolinggranular
gas in a channel. Shown is the maximum (rescaled) steady
2.5 a statedensity of thegas versustherescaled channellength ,
L
2 predictedbyEqs.(54)and(55). Panelbfocusesonavicinity
of the supercritical bifurcation point =2π.
w 1.5 L
1
0.5
solution (54) at m /2 is
0 | |≪L
0 5 10 15 20
3 m
m w0(m) L cosh−2 , (60)
≃ 8 2
r
(cid:16) (cid:17)
6
b where, for convenience, we have written the solution on
5
L 4 the interval /2 < m < /2 and used the approxi-
Ρ −L L
H 3 mate equality f 6/ . To calculate the asymptotics
g10 2 ofEq.(54)at mh i ≃1,wLecandealdirectlywithEq.(46)
o 1 and neglect th|e |f≫2 term. The solution of the resulting
l 0
elementary equation is a linear combination of em and
-1
0 5 10 15 20 e−m. The two arbitrary constants can be determined
from the two conditions at m = /2: df/dm = 0 and
x | | L
w f/ f =w , where w is given by Eq. (59).
0 min min
≡ h i
We obtain
FIG. 6: Same as in Fig. 4, but for = 19.869. Notice the p
L
logarithmic scale in panel b. w (m) √24 e /2cosh( /2 m), m 1. (61)
0 −L
≃ L L −| | | |≫
Note that the asymptotes (60) and (61) coincide in their
E 72 exp( ). The maximum value of w(m) is common region 1 m /2, where each of them
| |≃ −L ≪ | | ≪ L
yields
3
wmax ≃r 8L. (58) w0(m)≃√6Le−m. (62)
Toobtainthe minimumvalueofw(m) (thatcorresponds Notethat w (m) 6/ isdeterminedbytheasymp-
0
h i≃ L
tothemaximumvalueofthedensity),itisconvenientto tote (60). We compared the asymptotes (60) and (61)
p
use the exact relation w =b/ f and calculate the with the numerical solution, shown in Fig. 3, at a late
min
h i
asymptotic value of b at E 1, or 1. The result time τ = 72. Employing the periodic BCs, we shifted
| | ≪ p L ≫
is the numerical solution in m so that the maximum of
w(m,τ = 72) is at m = 0. One can see that the agree-
w √24 e /2, (59) ment is excellent.
min −L
≃ L
As higher w corresponds to a lower gas density, the
ByvirtueofEq.(53),theasymptoticsofthesteadystate regionofthe maximumofw correspondstoa hole inthe