Energy decay in a granular gas collapse

An inelastic hard ball bouncing repeatedly off the ground comes to rest in finite time by performing an infinite number of collisions. Similarly, a granular gas under the influence of external gravity, condenses at the bottom of the confinement due to inelastic collisions. By means of hydrodynamical simulations, we find that the condensation process of a granular gas reveals a similar dynamics as the bouncing ball. Our result is in agreement with both experiments and particle simulations, but disagrees with earlier simplified hydrodynamical description. Analyzing the result in detail, we find that the adequate modeling of pressure plays a key role in continuum modeling of granular matter.


Introduction
Force free granular gases were intensively studied using the tool of statistical mechanics and kinetic theory and by now there is a large body of knowledge in this field (e.g. [1,2] and many references therein). One of the first and at the same time most significant results concerns the kinetics of the energy dissipation in a homogeneous granular gas, T t T 1 t 0 2 = + t ( ) ( ) , known as Haff's cooling law [3], describing the decay of granular temperature, which can be defined as the mean kinetic energy fluctuations of the particles, analogous to the ordinary temperature. Here, T 0 is the temperature at time t=0 and τ is a time scale whose value depends on particle number density, T 0 and the coefficient of restitution, ε. The coefficient of restitution, in turn, characterizes the degree of inelasticity for a pair of colliding particles, v v e º -¢ , where v and v ¢ are the normal components of the relative velocities of the particles before and after the collision, respectively.
With the assumption of a constant coefficient of restitution, the homogeneous Boltzmann equation has a scaling solution [4], that is, the distribution, p(c) of the velocities reduced by the thermal velocity, c v v th º , is invariant in time. Thus, phenomenologically, in the homogeneous cooling state the granular gas is invariant, except that velocities decrease due to T t t 2 -( ) [3]. The behavior changes drastically when the granular gas evolves under the influence of an external gravity field: starting from a homogeneous initial state, the gas collapses rapidly and assumes a dense state at the bottom of the confinement. A careful analysis reveals that the collapse is a more complicated process [5][6][7][8]. In particular, it was found that the process is a sequence of shock front scenarios whose spatiotemporal appearance is attributed to the field of Mach number and its relation to the fields of density and temperature [8], characterizing the collapse as a supersonic phenomenon. The last stage of its evolution is a diffusive scaling regime [5,6], characterized by a critical time, t c , when the granular material comes to rest.
The evolution of a granular gas near the collapse, that is, the last stage of the sedimentation process, was first addressed in [5]. By means of hydrodynamic equations and particle simulations it was found that the energy of the gas decays according to E t t t c -b ( ) ( ) , where 2 b = , that is, the energy of the gas vanishes in finite time.
This scaling is familiar from the decay of energy of a hard ball bouncing inelastically off a rigid ground. With initial velocity, v 0 , and gravity, g, after s impacts the velocity decays to v s 0 e . The sth impact takes place at time That is, according to [5], the kinetics of the cooling of a granular gas is the same regardless of whether the particles bounce in an isolated way on the ground or whether the particles perform gas-like molecular chaos dynamics during the process of sedimentation, characterized by the collision frequency according to the current temperature and density and the corresponding cooling due to Haff's law. Thus, according to [5], the interaction between the granular particles seems to be irrelevant for the decay of energy, which is a surprising result. However, the analysis in [5] is based on several strong assumptions: the used hydrodynamic equations assume the limit of small inelasticity, ε  1. The final height of the sediment was assumed large, characterized by a dimensionless parameter L N 1 with L x being the container width, N p the number of particles and σ the particle diameter. A second parameter, L  . Moreover, [5] assumes a dilute gas, that is, the hydrostatic pressure does not include any correction due to the volume fraction occupied by the particles. As the authors mention, this assumption is invalid when the system approaches the collapse, t t c  , that is, at the time when the scaling of energy is addressed. The argument of [5] for its justification is that for t t c  , the kinetic energy of the condensate is negligible as compared to the contribution of the gas particles.
More recently, Son et al [6] performed a quasi-2d experiment to study the energy decay of a granular gas under gravity. They confirmed the power-law decay, though, with a much larger exponent, 2 b > . In [6], it is argued that the friction among the particles and the walls might cause the energy of the system to decay faster than predicted in [5] where smooth interaction was assumed. However, by means of molecular dynamics (MD) simulations disregarding boundary effects, Kachuck et al [7] confirmed the experimental results reported in [6] quantitatively.
In the present paper, we reconsider the decay of energy due to the collapse of a granular gas under the influence of gravity, that is, the final stage of the sedimentation process by means of hydrodynamic simulations based on the Jenkins-Richman constitutive relations [9]. For the numerical integration, we use a high-order shock-capturing method which was shown to deliver reliable results up to high density close to jamming, also in the presence of shocks. For a validation of this method see [10]. Without using the simplifying assumptions of [5] discussed above, for various values of the coefficient of restitution, we recover the power-law decay of the energy reported in [5]. For the exponent we find, however, a larger value, 2 b > , in contrast to [5] but in agreement with experiments [6] and particle simulations [7].

Hydrodynamic model
We consider a granular gas of smooth inelastic hard disks of mass m and diameter σ. The hydrodynamic fields of number density, n r t ,  ( ), flow velocity, u r t ,   ( ), and temperature, T r t ,  ( ), obey [11]: and g P q , , and z   denote gravity, pressure, heat flux, and cooling rate due to dissipative collisions, respectively. The constitutive relations for pressure and heat flux are with p, , , and h g k standing for hydrostatic pressure, shear and bulk viscosity and thermal conductivity. Their explicit forms are provided by the Jenkins-Richman theory [9]: with the freezing packing fraction 0.69 f f = and the random close packing fraction 0.82 During the process of sedimentation, the system shows simultaneously regions of low and high density which make the numerical solution of the compressible Navier-Stokes equations, (3), a numerical challenge requiring a highly precise (but numerically expensive) weighted essentially non-oscillatory (WENO) scheme [13,14]. This type of finite volume scheme applies a linear combination of lower order fluxes to obtain a higher order approximation to the numerical flux in the conservation equations and employs the idea of adaptive stencils to automatically achieve high order accuracy and non-oscillatory behavior near discontinuities [13,14]. This numerical scheme has been shown to give reliable results for similar problems [10].
Note that as the sedimentation progresses in time, the top region of the system would reach zero density while the sediment at the bottom approaches dense packing. Numerical hydrodynamic codes cannot handle these limits, therefore, our code restricts density to 0.0001 0.9999 max   f f while preserving the total mass.
Initially, the gas is homogeneously distributed in a container of size L L 40 y x = , avoiding any influence on the bottom plate that could come from the top wall. The container is modeled with square uniform mesh using 2560 nodes. We apply reflecting boundary conditions at all boundaries, that is, there is no loss of energy due to impact against the walls.
In order to quantitatively compare our results with data from the literature, we performed simulations for values of ε used in the references discussed above. Moreover, we studied also further values of ε as well as different initial conditions so as to explore the influence of these parameters on the results and to point out relations between them.

Simultaneous collapse scenario
Consider the decay of energy due to the process of sedimentation for the system described above and consider the system layer-wise. To this end, we define the Lagrangian mass coordinate y t k ( ) such that the fraction k of the total mass is below this vertical coordinate. Figure 1 shows the evolution of internal energy, E t nT , for layers of particles corresponding to four different Lagrangian mass coordinates, k 0.6, 0.7, 0.8, 0.9 = [ ] , using 0.98 e = , with corresponding parameters 10 2 e =and 10 L = inside the validity rage described in [5]. All curves start at the same value, E 0 k ( ), due to homogeneous initial conditions with initial density 0.15 f = , zero initial velocity field and initial temperature T = 10 6 mm s -2 . At the final stage of the sedimentation process, when the sequence of shocks described in [8] has passed, t 0.5 s  (inset in figure 1), we confirm that each layer cools down according to a power law. The black dashed lines in figure 1 show the best fits with the expression with optimal parameters given in table 1, which is done using an interval of points close to the collapse time, all weighted equally and made on the linear scale, with the smaller standard deviation of the fitting parameter. The value of t k c indicates the time when the energy contained in the strip at a certain vertical position characterized by the Lagrangian coordinate y k would vanish entirely. From table 1 we see that t 1.1 c k » s is nearly independent of k which means that all layers of the system come to rest simultaneously. This finding which is in agreement with theory [5] and the experiment [6] may surprise at a first glance, given that particles initially starting at different heights may undergo a rather different and complex history, see [8]. It may be understood from the finite value of the heat conductivity: consider two adjacent strips, at y k and y k k This gradient would cause a heat flux such that E 0 k > . Consequently, the collapse must occur for all layers simultaneously. This is true for arbitrary initial conditions, as we checked numerically for several cases. That is, while the collapse time certainly depends on the initial conditions, it is independent of the Lagrangian mass coordinate. All regions in the system come to rest simultaneously as previous studies have found. However, the role of the conductive regime precursor of the granular collapse was not properly associated with this property in those studies.

Global kinetics of the collapse
Consider now the evolution of the total kinetic energy of the system. In order to compare the results with previous theoretical, experimental, and numerical work [5][6][7], we chose the parameters of the system accordingly. In particular, we match the parameters used in [5], , thus our range of parameters covers the parameters used in [5,7]. Figure 2 shows the evolution of the total internal energy obtained from a hydrodynamic simulation together with best fits to E t t t c -b ( ) ( ) at the final stage. Obviously, the transient shock scenario appears less pronounced than in figure 1 since the sharp fronts visible in figure 1 do not occur at the same time for different Lagrangian mass coordinate. Therefore, by averaging over the system and thus, over the Lagrangian mass coordinate, the shocks appear smoother in ] was found. The system studied in [5] concerns the elastic limit 1  e . This setup is less realistic since a coefficient of restitution 0.995 e = employed in [5] is difficult to achieve in experiments. However, the restriction to 1  e is necessary due to the limitations of the computational method used. Figure 3 (blue line) shows the evolution of  the internal energy for the system parameters used in [5] together with the best fit obtained for t c = 2.05 s and In an attempt to reproduce the prediction made in [5], Kachuck et al [7] performed a particle simulation using 0.993 e = . They obtained ] which agrees with our result for the same system (see figure 3), where we obtain Thus the results of the particle simulation [7] are consistent with our results in both cases, 1  e and moderare inelasticity. This agreement is not trivial since both simulation methods, particle simulations (MD) and computational hydrodynamics, are fundamentally different approaches. By this good agreement, thus, both approaches support mutually.
We wish to point out that the decay rate, β, for 1  e found independently by means of particle simulations [7] and computational hydrodynamics (figure 3) is significantly larger than that predicted in [5] using both hydrodynamic and particle simulations.

Dependence on inelasticity
It is plausible that the time of collapse, t c , depends on the coefficient of restitution characterizing the loss of mechanical energy due to particle collisions. This assertion is compatible with the plots of figures 2 and 3. For a more systematic evaluation of the dependence of the collapse characteristics, t c and β, we performed a series of hydrodynamic simulations for 0.95, 0.995 e Î [ ] whose results are shown in table 2. These results clearly  confirm the dependence of the time of collapse on the inelasticity of the particles: the smaller ε the shorter t c . A similar (but much weaker) dependence can be seen for b e ( ). We emphasize, that for all systems investigated, the results agree with data from the literature employing both numerical particle simulations (MD) [7] and experiments [6].

Non-trivial pressure dependence
For all systems studied, the results of computational hydrodynamics employing sophisticated numerical WENO methods capable to accurately capture shock dynamics deliver energy decay rates significantly larger than 2 b = , in agreement with experiments and particle simulations [6,7]. In contrast, simplified hydrodynamics Apart from the above hand-waiving argument stating that 2 b = would correspond to non-interacting bouncing balls, thus, particle-particle interactions are unimportant, the true contribution of particle-particle interactions remains obscure. Therefore, let us have a more systematic look to understand the contradiction: the strongest simplification made in [5] is the assumption of a dilute gas, which implies the perfect gas equation of state, p = nT. While this assumption may be justified to describe the early stages of the process of sedimentation, it is questionable for larger density. As a justification, in [5] it is argued that the kinetic energy in regions of high density approaching the close packing limit, can be neglected as compared to the contribution from particles above the condensate, in the dilute region of the system. In our description, we use more complex hydrodynamic equations [15], with the hydrostatic pressure given by (6), which takes corrections due to the volume occupied by the particles into account.
For a more quantitative argument, in figure 4 we plot the ratio of the pressure used in our model (6) and the hydrostatic pressure considered in [5], ], from extreme dilution to the random close packing limit. From figure 4 we see that the corrected pressure p d , exceeds the perfect gas pressure, p 0 by one order of magnitude for 0.6 f » . For still larger density, close to dense packing, the p d may exceed p 0 by about three orders of magnitude. Taking into account that for t t c  when the system approaches the collapse, 0.82 f  (which is the situation considered), the vast majority of particles is contained in the condensed part of the system. Based on these arguments, we believe that the increase of pressure compared to the perfect gas pressure may the reason for the result 2 b = [5] of the analytical limit.

Conclusions
The process of sedimentation of a granular gas under the action of an external gravity field is a complex multistage process. In the present paper we focused onto the final stage starting after the supersonic shock scenario [8] has passed. This final stage comes along with the decay of kinetic energy according to E t t t c -b ( ) ( ) , which implies that the system comes to rest at finite time, t c , by completely dissipating its kinetic energy. This process was investigated before by means of simplified hydrodynamic equations [5] predicting the collapse scenario for t t c  , which was confirmed experimentally [6] and by means of particle simulations [5,7]. While there is qualitative agreement about the kinetics of the collapse, the quantitative description of the energy decay characterized by the exponent β is controversial: the analysis of simplified hydrodynamic equations [5] predicts 0.9-1.0 b = which is the same value as for isolated dissipative particles bouncing on a solid plane, that is, collective particle-particle interaction in the gas seems irrelevant for the collapse scenario. This prediction is, however, in conflict with experimental results [6] as well as particle simulations (MD) [7] that find always significantly larger values.
In the present paper, we investigated the collapse by means of an advanced hydrodynamic description [9, 11] using highly precise numerical WENO algorithms which allow to reliably simulate systems where dilute (gas-like) and dense (almost condensed) regions co-exist. Our hydrodynamic simulations reproduced the collapse scenario in agreement with the results in the literature. By analyzing the system layer-wise using the Lagrangian mass coordinate, we confirmed and explained the theoretical prediction [5] and experimental finding [6] that the entire system collapses simultaneously, independent of the vertical coordinate.
Considering the global kinetics of the collapse, we studied the dependence of the exponent β on the inelasticity of the particles, characterized by the coefficient of restitution, ε. For simulations with parameters taken from the literature [5][6][7] we could quantitatively reproduce the results from experiments [6] and particle simulations [7] but could not confirm the results from a simplified hydrodynamic analysis [5] predicting 2 b = .
In a systematic investigation of the values of exponent β as a function of the coefficient of restitution ε, we noticed a dependence of the time of the collapse, t c , and β on the value of ε. However, in all cases studied, β significantly exceeded the predicted value of 2. In order to explain this finding, we analyzed the ratio, p p d 0 , of the pressure in the Jenkins-Richman theory used here and the ideal gas pressure employed in [5] and find that depending on the momentary density, this ratio may reach values of p p 1 d 0  (see figure 4). Note that the Jenkins-Richman kinetic coefficients employed in our numerical simulations contain dense gas corrections, in particular the cooling coefficient and the coefficient for thermal diffusion, entering the terms that govern the latest stages of the energy decay. Therefore, we believe that the reason of the differences between the prediction in [5] as compared to the results in [6,7] is due to the dense gas corrections. This finding suggests that the correction of pressure [15] plays a key role in hydrodynamic modeling of granular matter.
While we see a strong dependence of p p d 0 on density, the dependence on inelasticity, ε is much weaker (figure 4), which may explain the weak dependence of the parameters β and t c on the coefficient of restitution shown in table 1.