ebook img

Attempted density blowup in a freely cooling dilute granular gas: hydrodynamics versus molecular dynamics PDF

0.65 MB·
Save to my drive
Quick download
Download
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.

Preview Attempted density blowup in a freely cooling dilute granular gas: hydrodynamics versus molecular dynamics

Attempted 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

See more

The list of books you might like

Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.