Temperature of a granular gas with regard to the stochastic nature of particle interactions

Granular gases, that is, ensembles of dissipatively interacting macroscopic particles may be considered as a model for granular matter in the dilute limit. The dynamics of granular gases is characterized by ballistic trajectories of the particles, interrupted by instantaneous dissipative collisions. In physical terms, the duration of collisions are considered negligible as compared with the time of mean free flight. This paradigm may serve as a definition of the term 'granular gas' in contrast to other granular systems where the particles may interact via permanent or long-lasting contacts. Besides of its practical importance for the description of natural phenomena, such as clouds of dust or the rings of the large planets of our solar system, granular gases have been intensively studied for their rich phenomenology.

From the point of view of theoretical physics, of particular interest is force-free granular gas, that is, an spacially infinite system in the absence of external gravity and boundaries due to confinement. For such systems, statistical physics provides powerful tools like the Liouville- and Boltzmann equation which may be adopted to account also for dissipative particle interactions [3, 4]. Subsequently, the full hydrodynamics can be derived, including the corresponding transport coefficients, e.g. by exploiting the Chapman–Enskog approach [4–6].

Because of the dissipative nature of particle interaction, granular gases are always in non-equilibrium, giving rise to many interesting phenomena, such as non-Maxwellian velocity distribution, e.g. [7], correlations of the velocities, e.g. [8], anomalous diffusion [9] and finally the instability of the homogeneous state in the long-time evolution [10] which may be transient, depending on the details of the particle interaction [11].

Starting from a homogeneous state at the time of initialization the evolution of a granular gas follows certain stages: (a) in a fast first stage, the distribution function of the gas rapidly converges to its native velocity distribution: on a time scale corresponding to a few collisions per particle, the distribution looses memory of the initial distribution and becomes essentially near Maxwellian, except in the high speed tail, where convergence is much slower [12] (mathematically it takes infinite time). (b) In the second stage of its evolution, the granular gas stays homogeneous, while the absolute velocity of the particles decays in time which may be expressed by a decay of temperature [13]. Therefore, this stage is called homogeneous cooling state (HCS). The velocity distribution and in particular temperature (its second moment), describe the HCS exhaustively, that is, the time dependence of the one-particle distribution function is solely a function of temperature. The distribution differs from the Maxwellian in two main points: the main part of the distribution, that is for v ∼ v_T , where v_T is the thermal velocity, is distorted which may be described by a low-order Sonine expansion [14]. For large velocities, v ≫ v_T , the distribution reveals an overpopulated high-energy tail, $f(v)\sim \exp (-v)$ [7, 15] in difference to $f(v)\sim \exp (-v^2)$ expected for molecular gases. The high-energy tail was found first from kinetic arguments [15] and numerical simulations and later from a high-order Sonine expansion [16]. The properties of granular gases in the HCS have been very intensively studied, e.g. [7, 17, 18], revealing in particular the breaking of equipartition in frictional [19, 20] or polydisperse [21–24] granular gases. The HCS has also been employed in derivations of Green–Kubo relations for granular gases [25, 26]. (c) The third stage of the evolution is characterized by the instability of the HCS and the formation of dense clusters, due to localized fluctuations of the density, associated with the instability of the shear mode [10, 27], for large enough systems [10, 27–29]. Those fluctuations, by locally modifying the collisions frequency (and thus energy dissipation), break the pressure balance and generate a positive feedback [10]. (d) Depending on the details of particle interaction, the cluster state may become unstable [11] because for small impact velocity, defined as the component of the relative velocity of two colliding particles along the direction joining their center, the collisions become less and less dissipative expressed by the coefficient of restitution approaching unity. In this case, the clusters resolve and the granular gas turns back to a homogeneous state at very low temperature. Resolution of the clusters may be observed for gases of viscoelastic spheres [11] as well as for charged particles and similar cases [30–32].

Whether the stages of the evolution of a granular gas occur on separate time scales depends mainly on the extend of energy loss due to dissipative collisions, expressed by the coefficients of restitution. In this paper we assume that the dissipation is small enough such that the HCS (stage (b)) exists for a certain time interval. That is, we assume that there exists an interval in time where the temperature is the only system characteristics which varies in time, before at stage (c) correlations of density and velocity become important.

Apart from particle number density and initial temperature, the key parameter, characterizing the kinetics of granular gases is the coefficient of restitution, ε, quantifying the loss of energy of the relative particle velocity due to collisions. While in most references it is assumed that ε is a material constant, detailed analysis of the collision process as well as experiments show clearly that the coefficient depends on the impact velocity which was mentioned in the literature already long ago, e.g. [33–36]. Moreover, even for virtually perfect spheres like ball bearings, in experimental measurements one observes significant scatter because of microscopic surface asperities [37–40]. This scatter may be expressed by considering the coefficient of restitution as a fluctuating quantity. The aim of this paper is to address the question how the stochastic nature of dissipative particle interactions influences the decay of temperature of a force-free granular gas as a function of time.

Throughout this paper we assume that the particle size and material properties as well as the thermal velocity and the number density of the gas are chosen such that the dynamics can be described by a purely repulsive coefficient of restitution. In particular, we require that the characteristic impact velocity is large enough to neglect van der Waals cohesive forces and small enough to avoid plastic deformation and brittle of the particles.

The structure of the paper is as follows. In section 2, some general results concerning the decay of temperature in the HCS are discussed. Section 3 presents a numerical and analytical analysis of the temperature behavior in the case of a stochastic, velocity-dependent coefficient of restitution. Section 4 is devoted to the study of the velocity distribution function, employing a slightly simplified (Laplacian) distribution for the coefficient of restitution. Finally, some conclusions are presented in section 5.

2.1. Constant coefficient of restitution and Haff's law

The inelasticity of a collision between two particles may be characterized by the coefficient of restitution, ε, relating the normal component of the particles' relative pre-collisional velocity, v'_ij ≡ v'_i − v'_j, whose norm defines the impact velocity, to that of the post-collisional relative velocity, v_ij,

$$\begin{equation} v_{ij} = -\varepsilon v_{ij}^{\,\prime}\,,\quad 0\leqslant \varepsilon\leqslant 1. \end{equation} \tag{ 1 }$$

where

$$\begin{equation} v_{ij}^{\prime} = \left(\vec{v}_i^{\,\prime}-\vec{v}_j^{\,\prime}\right)\cdot \vec{e}_{ij}\,,\quad v_{ij}= \left(\vec{v}_i-\vec{v}_j\right)\cdot \vec{e}_{ij}, \end{equation} \tag{ 2 }$$

$\vec {v}_i^{\,\prime }$, $\vec {v}_j^{\,\prime }$ are the pre-collisional velocities, $\vec {v}_i$, $\vec {v}_j$ are the post-collisional velocities, and

$$\begin{equation} \vec{e}_{ij}=\frac{\vec{r}_i-\vec{r}_j}{\left|\vec{r}_i-\vec{r}_j\right|} \end{equation} \tag{ 3 }$$

is the unit vector at the instant of contact. Using these definitions we obtain the post-collisional velocities via

$$\begin{equation} \vec{v}_{i}=\vec{v}_i^{\,\prime}-\frac{\left( 1+\varepsilon \right)}{2}v_{ij}^{\,\prime }\vec{e}_{ij}. \end{equation} \tag{ 4 }$$

Because of the dissipative nature of the collisions, characterized by ε < 1, the absolute velocities of the particles decay on average, which may be expressed by a decay of temperature. If we define temperature in analogy with the temperature of molecular gases,

$$\begin{equation} T=\textstyle{\frac{1}{3}}m \overline{v^{2}}, \end{equation} \tag{ 5 }$$

where m denotes the mass of a particle, and an overline denotes an average over the particles, the evolution of temperature is given by

$$\begin{equation} \frac{\mathrm{d}T}{\mathrm{dt}} \sim -nd ^{2}(1-\varepsilon ^{2})T^{3/2}, \end{equation} \tag{ 6 }$$

where n and d are the number density of the gas and the diameter of the particles, respectively. If the coefficient of restitution is assumed independent of the impact velocity, ε = const, we obtain Haff's law [13],

$$\begin{equation} T(t) = \frac{T_0}{\left(1+t/\tau_0\right)^2}\quad {\rm with}\quad \tau_0^{-1} \propto n d^2\left(1-\varepsilon ^{2}\right) \sqrt{T_0}\, \end{equation} \tag{ 7 }$$

where T₀ is the initial temperature, that is asymptotically T(t) ∼ t⁻². Haff's law is the most simple expression for granular cooling. The asymptotic form has been confirmed experimentally, where gravity was compensated by diamagnetic levitation [41]. Note that Haff's law is valid after the relaxation of the initial (Maxwell-) distribution to the native distribution; during this very short relaxation, the cooling is dominated by exponential terms [42]. Haff's law seems to be rather general; it was shown to hold true for smooth spheres in dimension 2–6 [43] as well as for a gas of needles [44]. Even beyond the HCS, that is at later times when significant correlations emerge, Haff's law is a fair approximation [18].

Somewhat unclear is the influence of the shape of the particles on the cooling law. While Haff cooling is observed for needles in 2d [44] and ellipsoids in 3d [45], in simulations of 'elongated' sharp edged particles T ∼ t^−5/3 was found [46].

2.2. Temperature decay due to an impact-velocity dependent coefficient of restitution

While the assumption ε = const simplifies the kinetic theory significantly, for realistic particles, the coefficient of restitution is a function of the impact velocity, ε = ε(v). For non-adhesive viscoelastic spheres, ε(v) can be obtained by integrating Newton's equation of motion describing the dynamics of the collision when the (normal component of the) interaction force is given by a generalized Hertz law [47]:

$$\begin{eqnarray*} &&F=\min \left(0, -\frac{Y\sqrt{d}}{3\left(1-\nu_{\mathrm{p}}^2\right)}\xi^{{3}/{2}}-\frac{3}{2}A\frac{Y\sqrt{d}}{3\left(1-\nu_{\mathrm{p}}^2\right)}\sqrt{\xi}v_{12} \right), \end{eqnarray*}$$

where $\xi \equiv 2 - \left |\vec {r}_1 -\vec {r}_2 \right |$ is the compression of the colliding spheres ( $\vec {r}_1$ and $\vec {r}_2$ denoting their positions), Y is the Young modulus, ν_p is the Poisson ratio and the dissipative constant A is a function of the elastic and viscous parameters. The coefficient of restitution is then given by $\varepsilon =- \frac {v_{12}(t_{\mathrm {c}})}{v_{12}(0)}$, where t_c is the duration of the collision, determined by requiring $\frac {\mathrm {d}^2 \xi (t_{\mathrm {c}})}{\mathrm {d} t^2}=0$. This complicated problem was solved rigourously [48] with the solution

$$\begin{equation} \varepsilon ({v})=1+\sum_{i=0}^{\infty }h_{i}\left( \sqrt{\frac{3A}{2}}\left(\frac{2Y\sqrt{d}}{3m(1-\nu_{\mathrm{p}}^2)} \right) ^{{1}/{5}} v^{{1}/{10}} \right) ^{i}, \end{equation} \tag{ 8 }$$

h_i being analytically known, universal constants. For realistic materials, equation (8) agrees well with experiments, e.g. [49], while a constant coefficient of restitution disagrees even with a dimension analysis [50]. From equation (8) follows that ε depends on the temperature in equation (6). In this case, the temperature decays due to

$$\begin{equation} \frac{\mathrm{d}T}{\mathrm{dt}}\propto -n d^2T^{{8}/{5}}\,, \end{equation} \tag{ 9 }$$

with the solution

$$\begin{equation} T(t) = \frac{T_0}{\left(1+t/\tau_0^\prime\right)^{{5}/{3}}}, \end{equation} \tag{ 10 }$$

that is, asymptotically T(t) ∼ t^−5/3.

In many cases, temperature decays algebraically with time, but this is not always the case: for ε due to

$$\begin{equation} \varepsilon(v)= \left\{\begin{array}{@{}ll@{}} \varepsilon^*\quad & {\rm for}\quad v>v^*,\\ 1 & {\rm for}\quad v\leqslant v^*, \end{array}\right. \end{equation} \tag{ 11 }$$

temperature decays logarithmically slow [31]:

$$\begin{equation} T(t) \sim 1/\log(\kappa t), \end{equation} \tag{ 12 }$$

where κ = const. The same behavior was found for particles with long range (electrostatic) repulsive interaction [30–32] while for attractive long range interaction enhanced cooling (compared with Haff's law) was reported [32].

Note that for the derivation of equations (7), (10) and (12) it was assumed that the state of the system is described by the second moment of the velocity distribution. In general, however, due to the non-Maxwellian distribution, the kinetic state depends also on higher (even) moments [51, 52], that is,

$$\begin{equation} \frac{\mathrm{d}T}{\mathrm{dt}}=F\left(T, \overline{v^4}, \overline{v^6}, \ldots\right)\,. \end{equation} \tag{ 13 }$$

For the cases discussed above, it was shown [53, 54] that the deviations of the distribution function from the Maxwellian, expressed by the Sonine coefficients, changes the decay rate, τ and τ', but not the functional form of the function T(t). The other type of deviation from the Maxwellian, that is, the overpopulated high-energy tail [15] leaves also its fingerprints in the temperature decay law, but its effect is very small [12].

2.3. Decay of temperature beyond the homogeneous cooling state

3.1. Stochastic coefficient of restitution

Realistic granular particles are certainly not perfect spheres. Even nearly spherical particles like glass spheres or precision steel spheres as used in ball-bearings, reveal surface asperities. These asperities lead to a surface roughness and in turn to a stochastic behavior of the particles in collisions which may be expressed by a stochastic coefficient of restitution. Even for microscopically small asperities, it was shown that the fluctuations of the coefficient of restitution are significant [40]. Based on experiments consisting of 2 × 10⁶ single impacts, performed by means of a robot [40] it was found that the functional form of the stochastic coefficient of restitution is strongly non-Gaussian. Instead, its probability distribution is Laplace-like, that is, composed of two exponential functions. Notably, the fluctuations observed in experiments exhibited a remarkable asymmetry:

$$\begin{equation} \rho_v(\varepsilon)=\left\{\begin{array}{@{}ll@{}} \displaystyle\frac{1}{\sigma_{p}(v)+\sigma _{m}(v)}\mathrm{e}^{-\frac{\varepsilon -\varepsilon _{{\rm max}}(v)}{\sigma _{p}(v) }} & {\rm if~} \varepsilon \geqslant \varepsilon _{{\rm max}}(v),\\ \displaystyle\frac{1}{\sigma _{p}(v)+\sigma _{m}(v)}\mathrm{e}^{\frac{\varepsilon -\varepsilon _{\max }(v)}{\sigma_{m}(v) }} & {\rm else}\,, \end{array}\right. \end{equation} \tag{ 14 }$$

where ε_max(v) is the value of ε where ρ_v(ε) adopts its maximal value. It was found that the function ε_max(v) is the same as for viscoelastic spheres, equation (8).

The probability distribution, equation (14), was experimentally found and confirmed by simulations where the surface roughness was modeled by a large number (∼ 10⁶) of microscopic spheres covering the surface of the sphere [40]. Figure 1 is obtained from experimental data [40] where the coefficient of restitution of a particle in contact with a flat plate is investigated as a function of the impact velocity. It shows the normalized frequency ρ(ε) of measurements of a certain value of ε for several small intervals of the impact velocity v.

Zoom In Zoom Out Reset image size

The mathematical form of equation (14) is a Laplace-probability density

$$\begin{equation} \rho^{\mathrm{L}}(\varepsilon)=\frac{1}{2\sigma }\mathrm{e}^{\frac{-\left|\varepsilon-\left<\varepsilon\right>\right|}{\sigma}} \end{equation} \tag{ 15 }$$

with two generalizations:

• (a)  
  The mean value $\left <\varepsilon \right >$ is replaced by a decaying function of the impact velocity, ε_max(v), to account for the physical fact that the coefficient of restitution increases with decreasing impact velocity. This behavior is well known for dry viscoelastic particles and for impact velocity small enough to avoid plastic deformations and large enough to neglect surface effects like adhesion, e.g. [57, 58].
• (b)  
  The Laplace distribution, equation (15), is symmetric with respect to $\left <\varepsilon \right >$, while the distribution density, equation (14) is asymmetric, characterized by two (impact velocity dependent) standard deviations σ_p(v) and σ_m(v), for ε > ε_max and ε < ε_max, respectively. The asymmetry can be explained considering the mechanism of energy transfer between the (external) linear degree of freedom and internal degrees of freedom during the collision [40]. Due to this asymmetry, strictly speaking, the maximum of the distribution function, ε_max(v) is different from the expectation value $\left <\varepsilon \right >$, however, both values are very close.

In the subsequent analytical and numerical calculations we will assume that the coefficient of restitution obeys the probability density function equation (14).

Fluctuating coefficients of restitution were also introduced in other contexts such as to model heating through collisions in a vertically vibrating granular system [59, 60], to describe a gas of one-dimensional particles accounting for internal degrees of freedom [61, 62] and in gases on Maxwellian molecules [63] to investigate the high velocity tail of the distribution function. More general random collision rules were considered in [64]. All these models are rather different from our collision model, equation (14), in that the functional form of the coefficient of restitution was not based on considerations of the physical processes during collisions.

3.2. Decay of temperature—kinetic theory

At the kinetic level, the system is described by a distribution function for the velocities, $f(\vec {v})$, by mean of which one can express the number density field, and temperature, as

$$\begin{equation} n=\int f\left( \vec{v}\,\right) \mathrm{d}\vec{v}, \end{equation} \tag{ 16 }$$

$$\begin{equation} T=\frac{1}{n}\int f\left( \vec{v}\,\right) \frac{m v^2}{3} \mathrm{d}\vec{v}. \end{equation} \tag{ 17 }$$

The equation of motion for T is given by

$$\begin{equation} n\frac{\partial T}{\partial t}=-\gamma , \end{equation} \tag{ 18 }$$

where γ is the cooling coefficient, describing the rate of decay of the kinetic energy (multiplied by ${2}/{3}$), due to the inelasticity of the collisions. Using the velocity transformations, equation (4), the change of energy ΔE in a collision between particles with velocities $\vec {v}_1$ and $\vec {v}_2$, along the direction spanned by the unit vector $\vec {e}$, and characterized by a coefficient of restitution ε, is given by

$$\begin{equation} \Delta E =-\frac{\left(1-\varepsilon^2\right)}{2} m\left( \vec{v}_{12}\cdot\vec{e}\,\right)^2\,. \end{equation} \tag{ 19 }$$

In addition, the collision frequency ν₁₂ of such collisions is

$$\begin{equation} \nu_{12}=d^2 f\left(\vec{v}_1\right) f\left( \vec{v}_{2}\right) \left|\vec{v}_{12}\cdot\vec{e}\,\right|\mathrm{d}\vec{v}_{1}\mathrm{d}\vec{v}_{2}\mathrm{d}\vec{e}, \end{equation} \tag{ 20 }$$

where d is the diameter of the particles. Therefore, the cooling coefficient is given in terms of the distribution function by

$$\begin{equation}\fl \gamma =\frac{m}{12}d^{2}\int \int \int \int \rho_{\left| \vec{v}_{12} \cdot \vec{e}\,\right| } \left(\varepsilon\right) f\left(\vec{v}_1\right) f\left( \vec{v}_{2}\right) \left(1-\varepsilon ^{2}\right) \left( \vec{v}_{12}\cdot\vec{e}\,\right)^{3} H\left( \vec{v}_{12}\cdot \vec{e}\,\right) \mathrm{d}\vec{v}_{1}\mathrm{d}\vec{v}_{2}\mathrm{d}\vec{e}\,\mathrm{d}\varepsilon, \end{equation} \tag{ 21 }$$

where a factor 1/2 has been included to avoids double counting of collisions, $\rho _{v } \left (\varepsilon \right )$ denote the conditional probability for a coefficient of restitution ε, given a normal relative velocity v, H is the Heaviside step function and the integration range of ε is [ − ∞,+∞]. Note that the distribution function depends on the distribution $\rho _{v} \left ( \varepsilon \right )$. In order to take into account this dependence, one needs to solve the Boltzmann equation pertaining to the system. For the moment, we will for sake of simplicity neglect this dependence and assume that the deviation of the velocity distribution function from the Maxwell distribution

$$\begin{equation} f^{\left(0\right)}=n\left( \frac{m}{2 \pi T}\right)^{{3}/{2}} \mathrm{e}^{-\frac{m v^2}{2T}} \end{equation} \tag{ 22 }$$

is small. This assumption is reasonable since the range of interest for the coefficient of restitution is of the region $\varepsilon\gtrapprox 0.97$ due to the experimental results on which this paper is based [40], see figure 1. Recall that the value of the second Sonine coefficient (corresponding to the inelastic correction to the Maxwellian distribution, cf section 4) for a granular gas with ε = const is a₂ < 0.01 in this region [14, 54]. In section 4 we will consider a slightly simplified distribution $\rho \left ({\varepsilon }\right )$ and show that the deviations from the Maxwellian are small indeed.

If we approximate the distribution function by a Maxwellian, $f\left (\vec {v}\right ) \approx f^{(0)}\left (\vec {v}\right )$, the cooling coefficient γ reads

$$\begin{equation}\fl \gamma =\frac{m}{12} d^2 \int \int \int \left(1-\left\langle\varepsilon ^{2}\right\rangle\right) f^{(0)}\left( \vec{v}_{1}\right) f^{(0)}\left( \vec{v} _{2}\right)\left( \vec{v}_{12}\cdot\vec{e}\,\right)^{3}H\left( \vec{v}_{12}\cdot \vec{e}\,\right) \mathrm{d}\vec{v}_{1}\mathrm{d}\vec{v}_{2}\mathrm{d}\vec{e}, \end{equation} \tag{ 23 }$$

where brackets denote an average over ε. It appears that the evolution of the temperature deviates from the decay corresponding to the case of non-fluctuating coefficient of restitution due both to the variance of the distribution, and to the asymmetry of the distribution function, ρ_v(ε) (i.e. 〈ε²〉 ≠ ε²_max). While the former has negligible effect, the latter has a qualitative influence on the asymptotic evolution of the temperature. By using the coefficients σ_m(v) and σ_p(v) obtained from simulation data (cf section 3.3 below), the average coefficient of restitution corresponding to the distribution (14):

$$\begin{equation} \left\langle \varepsilon \right\rangle(v) =\varepsilon _{\max }(v)-\left[ \sigma _{m}(v)-\sigma _{p}(v)\right] \end{equation} \tag{ 24 }$$

is found to be well described by a function of the form

$$\begin{equation} \left\langle\varepsilon\right\rangle =1-\delta^{\ast }v ^{\eta}\quad{\rm with}\quad \delta^{\ast }=0.053\,17\,;\quad \eta =\frac{4}{25} \end{equation} \tag{ 25 }$$

(see figure (2)). A straightforward calculation then yields

$$\begin{equation} \gamma =\frac{288}{625} 2 ^{\frac{33}{50}}\sqrt{\pi } \Gamma \left( \frac{2}{25}\right) n^2d^2\frac{1}{m^{\frac{29}{50}}}T^{\frac{79}{50}} \delta^{\ast }, \end{equation} \tag{ 26 }$$

i.e. solving equation (18):

$$\begin{equation} T\left( t\right) =\frac{T_{0}}{\left[ 1+\frac{4176}{15625}\left( 2^{ \frac{33}{50}}\right) \Gamma \left( \frac{2}{25}\right) \sqrt{\pi }nd ^{2}\delta^{\ast }\left( \frac{T_{0}}{m} \right) ^{\frac{29}{50}}t\right] ^{\frac{50}{29}}}, \end{equation} \tag{ 27 }$$

that is,

$$\begin{equation} T=\frac{T_{0}}{(1+\frac{t}{\tau^\rho} )^{{50}/{29}}} \ \ {\overset{t\to \infty}{\longrightarrow}} \ \sim t^{-{50}/{29}} \end{equation} \tag{ 28 }$$

with the constant τ^ρ defined by equation (27). Comparing this result with the decay of temperature of a gas of viscoelastic particles [48, 65], $T\sim t^{- {5}/{3}}$, equation (10), we obtain that the fluctuating coefficient of restitution lets the gas cool noticeable faster but slower than gases of particles interacting via ε = const, where $T\sim t^{- {1}/{2}}$, according to Haff's law [13].

Zoom In Zoom Out Reset image size

3.3. Decay of temperature—numerical simulation

3.3.1. Numerical description of the stochastic coefficient of restitution

For the direct simulation Monte Carlo (DSMC) simulation we need random numbers distributed according to $\rho _v\left (\varepsilon \right )$ as given by equation (14) where ε_max(v) in turn is given by equation (8). Using the inverse transformation method [66], these random numbers have to be computed from standard random numbers, x, equally distributed in the interval x∈(0,1), available on the computer. The cumulative distribution function

$$\begin{equation} \varrho\left(\varepsilon\right) = \int\nolimits_{-\infty}^\varepsilon \rho(\epsilon) \mathrm{d}\epsilon \end{equation} \tag{ 29 }$$

maps each value of ε to a number in the interval (0,1). So if ϱ can be inverted, ϱ⁻¹(x), obeys the distribution ρ if x is an equally distributed random number, x∈(0,1).

For later use, we write the distribution function more general,

$$\begin{equation} \rho_v(\varepsilon )=\left\{\begin{array}{@{}ll@{}} \mathrm{e}^{a\varepsilon +b}, & \varepsilon <\varepsilon_{\rm max}(v),\\ \mathrm{e}^{c\varepsilon +d}, & {\rm else}, \end{array}\right. \end{equation} \tag{ 30 }$$

where the coefficients a,b,c and d follow from equation (14). Its cumulative distribution function reads

$$\begin{equation} \varrho_v(\varepsilon )=\int\nolimits_{-\infty}^{\varepsilon} \rho_v(\epsilon) \mathrm{d}\epsilon= \left\{\begin{array}{@{}ll@{}} \displaystyle\frac{1}{a}\mathrm{e}^{a\varepsilon +b} & \varepsilon <\varepsilon _{\max }(v)\\ \displaystyle\frac{1}{a}\mathrm{e}^{a\varepsilon_{\max }+b}-\displaystyle\frac{1}{c}\mathrm{e}^{c\varepsilon _{\max }+d}+\displaystyle\frac{1}{c}\mathrm{e}^{c\varepsilon+d} & {\rm else} \end{array}\right\} \end{equation} \tag{ 31 }$$

with the inverse

$$\begin{equation} \varrho_v^{-1}(x)=\left\{\begin{array}{@{}ll@{}} \displaystyle\frac{-b}{a}+\displaystyle\frac{1}{a}\ln (ax), & \varrho^{-1}(x)<\varepsilon _{\max }(v),\\ \displaystyle\frac{-d}{c}+\displaystyle\frac{1}{c}\ln \left(cx-\frac{c}{a}\mathrm{e}^{a\varepsilon _{\max }+b}+\mathrm{e}^{c\varepsilon_{\max }+d}\right), & {\rm else}. \end{array}\right. \end{equation} \tag{ 32 }$$

In this way, we obtain random numbers distributed according to ρ(ε) from equidistributed random numbers x∈(0,1).

Note that the functions ϱ⁻¹_v depend on the impact velocity through $\varepsilon _{\max }(v)$. To obtain a random number according to ρ_v(), thus, we first determine $\varepsilon _{\max }(v)$ and then by means of a random number x∈(0,1), via equation (32) we obtain the desired random number. In principle, $\varepsilon _{\max }(v)$ could be computed via equation (8) by truncating the series at a certain order. However, this series in $v^{ {1}/{10}}$ converges very slowly [48]. Therefore, we prefer the Padé approximation of order [1/4] provided in [67]

$$\begin{equation} \varepsilon _{\max }(v)\approx\frac{1 + c_1 v_*}{1+b_1v_*+b_2v_*^2+b_3v_*^3+b_4v_*^4}\,,\quad v_{\ast }=\beta ^{1/2}v^{1/10} \end{equation} \tag{ 33 }$$

with the universal (material independent) constants c₁ = 0.501 086, b₁ = 0.501 086, b₂ = 1.153 45, b₃ = 0.577 977, b₄ = 0.532 178.

For demonstration, figure 3 shows a histogram of 10⁶ random numbers obtained by application of the described algorithm. The material constant β = 0.0467 was used, corresponding to the experimental values given by Montaine et al [40] for the collision of stainless steel spheres.

Zoom In Zoom Out Reset image size

3.3.2. Direct simulation Monte-Carlo details

To study the decay of temperature in the HCS, we apply DSMC [68, 69]. DSMC simulations of dilute granular gases have been first performed by Brey et al [17]. DSMC is by orders of magnitude faster than molecular dynamics (MD) (event-driven as well as force controlled MD), however, its main advantage compared to MD is that DSMC is not a particle simulation [70] but a tool to integrate the Boltzmann (or Boltzmann–Enskog) equation by means of a Monte-Carlo procedure. This allows us to compare the results of a DSMC simulation directly with the results obtained from the kinetic theory. When investigating a granular gas in the HCS, MD is restricted to almost elastic particles, $\varepsilon \lesssim 1$, and low density, to avoid formation of correlations of different kind, conflicting the homogeneous state. In contrast, DSMC works for all values of the coefficient of restitution and density by avoiding correlations. We wish to remark here that large part of the literature confirms that influence of the correlation of particle velocities due to inelastic collisions on the cooling rate is small, e.g. [71–75], which has, however, also been questioned [76].

In our simulation, at each time step a pair of (quasi-) particles is selected randomly with equal probability. A collision according to the collision rule, equation (4), was executed if

$$\begin{equation} \left|\vec{e}\cdot\vec{v}_{12}\right| >{\rm RND}\,7v_{T}, \end{equation} \tag{ 34 }$$

where RND is a uniformly distributed random number from the interval RND∈(0,1) and v_T is the thermal velocity. To compute the post-collisional velocities using equation (4) we determined a random unit-vector $\vec {e}$. The term 7v_T is an estimate for the smallest upper bound of the average normal velocity, see [69] for justification of this term. The simulation was performed with 10⁵ quasi-particles, initialized with a Maxwellian velocity distribution according to the thermal velocity T₀ = 1.

3.3.3. Decay of temperatures obtained from DSMC simulations

We simulated a granular gas of particles interacting via a fluctuating coefficient of restitution due to equation (14) in the HCS and computed the granular temperature as a function of time via equation (5). For large times, the temperature decays according to a power law, T ∼ t^−k, with the best fit k = 1.71 being in excellent agreement with the result obtained from kinetic theory, equation (28), $k= {50}/{29}\approx 1.7241$. Figure 4 shows the DSMC simulation data with a best fit of the function

$$\begin{equation} T=\frac{T_{0}}{(1+t/\tau )^{1.71}} \end{equation} \tag{ 35 }$$

with a fit value for τ together with the analytical result, equation (28). The asymptotic decay due to Haff's law [13] (for ε = const) and for a gas of viscoelastic particles [65] is also shown.

Zoom In Zoom Out Reset image size

The result deviates from Haff's law (k = 2) [13] as well as from the evolution of temperature for a gas of viscoelastic particles ($k= {5}/{3}$) [65]. Note that the deviation from the case of viscoelastic particles originates from the asymmetry of the probability function, ρ_v(ε). Although the median of the coefficient of normal restitution equals the expression for viscoelastic spheres [48, 65], the mean value deviates considerably because of the skewness of the functions.

The cooling coefficient γ depends on the distribution ρ_v(ε) explicitly, as shown by equation (21), but also implicitly, since the velocity distribution function, $f(\vec {v})$, obviously depends on the choice of ρ_v(ε). In section 3.2, we assumed that the deviation of the distribution function from a Maxwellian is small, based on (a) our numerical observations, section 3.3.3, and (b) the fact that for ε = const the second Sonine coefficient quantifying the deviation of $f(\vec {v})$ from the Maxwell distribution, is very small ($a_2\lesssim 0.01$) in the range of ε ≈ 0.98 being of interest here. In this section we wish to quantify the claim that the distribution function is close to a Maxwell distribution to—à posteriori—justify the above assumption. In this section we are not interested in the absolute value of the coefficient of restitution (characterized by ε_max), nor in the skewness of its distribution, ρ_v(ε), being both of major importance for the decay of temperature. Instead, here we consider the fluctuations of the coefficient of restitution, ε. To this end, we reduce the full distribution function, equation (14), to the Laplace distribution, equation (15), by neglecting the velocity dependence and the skewness of the distribution. Apart from these details, the Laplace distribution, equation (15), is of the same characteristic shape as equation (14), that is, it consists of two exponentials

$$\begin{equation} \rho^{\mathrm{L}}(\varepsilon)=\frac{1}{2\sigma }\mathrm{e}^{\frac{-\left|\varepsilon-\left<\varepsilon\right>\right|}{\sigma}} . \end{equation} \tag{ 36 }$$

The Boltzmann equation describing a granular gas in the HCS at the kinetic level reads, in the case of a velocity-independent (fluctuating) coefficient of restitution:

$$\begin{equation} \frac{\partial f}{\partial t}=\int B\left(f,f,\varepsilon \right) \rho \left(\varepsilon \right) \mathrm{d}\varepsilon, \end{equation} \tag{ 37 }$$

where $\rho \left (\varepsilon \right )$ is the (velocity independent) distribution of the coefficient of restitution, and $B\left (f,f,\varepsilon \right )$ is the Boltzmann collision operator corresponding to a fixed coefficient of restitution

$$\begin{equation}\fl B\left(f,f,\varepsilon \right) =d^{2}\int \int \left[ \frac{1}{\varepsilon^2 }f\left(\vec{v}_{1}^{\,\prime }\right)f\left(\vec{v}_{2}^{\,\prime }\right)-f\left(\vec{v}_{1}\right)f\left(\vec{v}_{2}\right)\right] \left(\vec{v}_{12}\cdot \vec{e}_{12}\,\right) H\left(\vec{v}_{12}\cdot \vec{e}_{12} \right) \mathrm{d}\vec{v}_{2}\mathrm{d}\vec{e}_{12} . \end{equation} \tag{ 38 }$$

The unit vector $\vec {e_{12}}$, the sets $\left \{ \vec {v}_{1}^{\,\prime },\vec {v}_{2}^{\,\prime }\right \}$, $\left \{ \vec {v}_{1},\vec {v}_{2}\right \}$ of the pre- and postcollisional velocities and the relative velocity, $\vec {v}_{ij}$, were introduced in equations (2)–(4).

In order to take into account the inelastic corrections to the Maxwellian distribution, we consider a normal solution of the Boltzmann equation, i.e. a solution depending on time through its (functional) dependence on the hydrodynamic fields. In the case of a homogeneous system, the only hydrodynamic field which is time dependent is the temperature, and we seek a similarity solution in the time-independent scaled velocity,

$$\begin{equation} c\equiv\sqrt{\frac{mv^{2}}{2T}}, \end{equation} \tag{ 39 }$$

that is, we scale the velocity by the thermal velocity to compensate the decay of temperature. While this is the standard procedure in the physical literature, the pure existence of the time-independent similarity solutions is not mathematically proven [77]. Assuming it exists, the Boltzmann equation, equation (37), in this case reduces to

$$\begin{equation} \frac{\gamma }{2nT}c\frac{\partial f}{\partial c}=\int B\left(f,f,\varepsilon \right) \rho \left(\varepsilon \right) \mathrm{d}\varepsilon. \end{equation} \tag{ 40 }$$

The standard way to solve the above equation is to represent the distribution function, $f(\vec {v})$, in form of a truncated series of Sonine polynomials [78], whose coefficients characterize the deviation of $f(\vec {v})$ from the Maxwell distribution:

$$\begin{equation} f\left(\vec{v}\right) =n\left(\frac{m}{2\pi T}\right) ^{\frac{3}{2}}\mathrm{e}^{- c^2} \sum_{i=0}^N a_i S_{\frac{1}{2}}^{i}\left(c^2\right), \end{equation} \tag{ 41 }$$

where N is the order of truncation, and the Sonine polynomials are defined by [78]

$$\begin{eqnarray*} &&S_{l}^{i}\left(x\right) =\sum\limits_{p=0}^{i}\frac{\Gamma\left( l+i+1\right) }{\left( i-p\right) !p!\Gamma\left( p+l+1\right) }\left( -x\right) ^{p}. \end{eqnarray*}$$

This method was first applied by Goldshtein and Shapiro [14, 54] to granular gases. Note that there are different truncation schemes to compute the Sonine coefficients a_i [16, 79–81]. Using the expansion equation (41) in equation (40) and in equation (21) for γ, and projecting back on the Sonine polynomials, one obtains a system of equations for the coefficients $\{a_{0},a_{1,}a_{2,}\ldots\}.$ This system is to be solved together with the orthogonality conditions expressed by the requirements equations (16) and (17) that the number density and temperature are given by the appropriate moments of the distribution function. Those conditions impose constraints on the coefficients {a_i}. In particular, due to the orthogonality properties of the Sonine polynomials, these requirements imply a₀ = 1 and a₁ = 0. Truncating the series, equation (41), to the first non-trivial order (N = 2) and assuming that the coefficient a₂ is small enough so that orders higher than linear can be neglected, one readily obtains the following expression for the correction:

$$\begin{equation} a_{2}=-\frac{32\left\langle \varepsilon ^{4}\right\rangle -48\left\langle \varepsilon ^{2}\right\rangle +16}{30\left\langle \varepsilon ^{4}\right\rangle +13\left\langle \varepsilon ^{2}\right\rangle -64\left\langle \varepsilon \right\rangle -81}, \end{equation} \tag{ 42 }$$

where brackets denote the average over the values of the coefficients of restitution. In the case of constant ε, the standard expression for a₂ is recovered [14, 54, 82]. Using the above expression, together with the expansion (41) in equation (21), the cooling coefficient reads, to linear order in a₂

$$\begin{equation} \gamma =4\sqrt{\pi }\left( \frac{(1-\left\langle \varepsilon ^{2}\right\rangle )\left( 32\left\langle \varepsilon ^{4}\right\rangle -48\left\langle \varepsilon ^{2}\right\rangle +16\right) }{30\left\langle \varepsilon ^{4}\right\rangle +13\left\langle \varepsilon ^{2}\right\rangle -64\left\langle \varepsilon \right\rangle -81}\right) \frac{nd ^{2}T^{{3}/{2}}}{\sqrt{m}}. \end{equation} \tag{ 43 }$$

In the case of a Laplace distribution, equation (15), γ is given by

$$\begin{eqnarray*} &&\gamma \left( \left\langle \varepsilon \right\rangle ,\sigma \right) =\frac{4 \sqrt{\pi }n^{2}d ^{2}T^{{3}/{2}}}{\sqrt{m}}\left( 1-\left\langle \varepsilon\right\rangle ^{2}-2 \sigma^2 \right)\left(1+\gamma^{(1)} \left( \left\langle \varepsilon \right\rangle ,\sigma \right) \right), \end{eqnarray*}$$

where

$$\begin{eqnarray*}&&\fl \gamma^{(1)}\left( \left\langle \varepsilon \right\rangle ,\sigma \right) =\frac{3}{2}\frac{ 32\left\langle \varepsilon \right\rangle ^{4}-48\left\langle \varepsilon \right\rangle ^{2}+16 +96\sigma ^{2}\left( 4\left\langle \varepsilon \right\rangle ^{2}+8\sigma ^{2}-1\right) }{30\left\langle \varepsilon \right\rangle ^{4}-13\left\langle \varepsilon \right\rangle ^{2}-64\left\langle \varepsilon \right\rangle -81+2\sigma ^{2}\left( 180\left\langle \varepsilon \right\rangle ^{2}+360\sigma ^{2}-13\right) }, \end{eqnarray*}$$

corresponds to the correction due to a₂ to the cooling coefficient. It can be checked that this correction is negligible (cf figure 5). For standard deviations σ < 0.1, this correction is lower than 4% even for low values of the mean coefficient of restitution, i.e. highly inelastic systems. For the value 〈ε〉 = 0.98 considered here, the correction is lower than 1%. In order to characterize the influence of the fluctuations, the ratio

$$\begin{equation} \Delta \gamma\equiv\frac{\gamma(\langle \varepsilon \rangle, \sigma)-\gamma((\langle \varepsilon \rangle, 0)}{\gamma((\langle \varepsilon \rangle, 0)} \end{equation} \tag{ 44 }$$

is plotted in figure 6 as a function of σ, for a mean coefficient of restitution 〈ε〉 = 0.98, together with simulation results for the same quantity, obtained by fitting the cooling curve to Haff's law. Agreement between simulation and theory is excellent, and the fluctuations are seen to have a rather weak influence, yielding a deviation of about 10% with respect to the non-fluctuating case, even in the case of relatively large fluctuations σ = 0.05. (Notice that in this case, the coefficient of restitution can exceed unity at very low probability.)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We considered the decay of temperature of a force-free homogeneous granular gas of particles interacting by instantaneous binary collisions, where we took the stochastic nature of particle interaction into account. To this end, we employed the velocity dependent conditional distribution for the coefficient of restitution described recently in [40], based on both large scale experiments and numerical simulations.

By means of kinetic theory we find a novel type of asymptotic decay of temperature following the law $T\sim t^{- {50}/{29}}$. This function is in contrast to both standard cases of temperature decay in dilute gases of purely repulsively interacting spheres, namely Haff's law T ∼ t⁻² corresponding to a gas where the particles interact via a constant coefficient of restitution, and to $T\sim t^{- {5}/{3}}$ corresponding to a gas of viscoelastic spheres. The result was checked by numerical DSMC simulations and very good agreement was obtained. Although the deviation of the new exponent from the exponent for gases of viscoelastic spheres is rather small, both cases can be clearly distinguished in DSMC simulations. That is, the stochastic nature of particle interaction leaves a clear footprint in the macroscopic properties of the gas.

The deviation of the exponent from the standard cases was shown to be associated with the asymmetry of the distribution function of the fluctuating coefficient of restitution. This asymmetry results in a difference between the (impact velocity dependent) mean value and the (viscoelastic) most probable value of the coefficient of restitution. Despite the obtained qualitative difference in the decay of temperature, the effect is quantitatively small, due to the small variances of the distribution for ε found in simulation and experiments [40].

In order to à posteriori justify the assumption of a Maxwellian velocity distribution employed for the analytical calculation (but not for the DSMC simulations!), we considered a simplified (velocity-independent) probability distribution of the coefficient of restitution of Laplace type. For this case, we solve the Boltzmann equation via a truncated Sonine expansion to obtain the velocity distribution of the gas particles. We could show that the correction to the Maxwellian distribution due to the inelasticity is negligible indeed, justifying the overmentioned assumption.

For non-spherical particles of more complex shape, stronger effects can be expected as a result of a wider distribution of the fluctuating coefficient of restitution.

This work was supported by the German Science Foundation (DFG) through the Cluster of Excellence 'Engineering of Advanced Materials'. DS thanks the Humboldt Foundation for support via the Alexander von Humboldt post-doctoral fellowship. We acknowledge support by Deutsche Forschungsgemeinschaft and Friedrich–Alexander–Universität at Erlangen–Nürnberg within the funding programme Open Access Publishing.

Substituting equation (22) into equation (23), the cooling coefficient reads

$$\begin{equation} \fl \gamma =\frac{m}{12}n^{2}d^{2}\left( \frac{m}{2\pi T}\right) ^{3}\int \int \int \left( 1-\left\langle \varepsilon ^{2}\right\rangle \right) \mathrm{e}^{-\frac{m}{2T}\left( v_{1}^{2}+v_{2}^{2}\right) }\left( \vec{v}_{12}\cdot \vec{e}\,\right) ^{3}H\left( \vec{v}_{12}\cdot \vec{e}\right) \mathrm{d}\vec{v}_{1}\mathrm{d}\vec{v}_{2}\,\mathrm{d} \vec{e}. \end{equation} \tag{ A.1 }$$

Using the form equation (25) in the above equation, neglecting the variance of $\rho \left ( \varepsilon \right )$ (the latter being of the order of the square of the small standard deviations σ_m and σ_p) and quadratic terms in δ* yields

$$\begin{equation} \fl \gamma =\frac{m}{6}n^{2}d^{2}\left( \frac{m}{2\pi T}\right) ^{3}\delta ^{\ast }\int \int \int \mathrm{e}^{-\frac{m}{2T}\left( v_{1}^{2}+v_{2}^{2}\right) }\left( \vec{v}_{12}\cdot \vec{e}\,\right) ^{3+\eta }H\left( \vec{v}_{12}\cdot \vec{e}\right) \mathrm{d} \vec{v}_{1}\mathrm{d}\vec{v}_{2}\mathrm{d}\vec{e}+\mathcal{O}\left( \delta ^{\ast 2}\right). \end{equation} \tag{ A.2 }$$

First, the integral over $\vec {e}$ is readily performed:

$$\begin{equation} \int \left( \vec{v}_{12}\cdot \vec{e}\right) ^{3+\eta }H\left( \vec{v}_{12}\cdot \vec{e}\,\right) \mathrm{d} \vec{e}=2\pi \frac{\left\vert \vec{v}_{12}\right\vert ^{3+\eta }}{4+\eta}. \end{equation} \tag{ A.3 }$$

Next, upon changing the variables $\vec {v}_1$ and $\vec {v}_2$ to

$$\begin{equation} \begin{array}{@{}ll@{}} \displaystyle \vec{u_{c}}=\frac{m}{2T} \left( \vec{v}_{1}+\vec{v}_{2}\right),\\ \displaystyle \vec{u}_{12}=\frac{m}{2T}\left( \vec{v}_{1}- \vec{v}_{2}\right), \end{array} \end{equation} \tag{ A.4 }$$

$$\begin{equation} \begin{array}{l} \gamma\!=\!\displaystyle\frac{m}{3\pi^{2}}n^{2}d^{2}\left(\displaystyle\frac{m}{2T}\right)^{-\displaystyle\frac{3\!+\!\eta}{2}}\delta^{\ast}\displaystyle\frac{1}{4\!+\!\eta }\int \int \int \mathrm{e}^{-\left( u_{c}^{2}+\displaystyle\frac{u_{12}^{2}}{2}\right) }\left\vert \vec{u}_{12}\right\vert ^{3+\eta }\mathrm{d}\vec{u}_{c}\mathrm{d}\vec{u}_{12}, \\ \hspace{0.9pc}=\displaystyle\frac{32\sqrt{2}}{3}\displaystyle\sqrt{\pi }n^{2}d^{2}\delta ^{\ast }\displaystyle\frac{2^{\eta}}{4+\eta} \Gamma \left( 3+\displaystyle\frac{\eta }{2}\right) \displaystyle\frac{T^{\displaystyle\frac{3+\eta }{2}}}{m^{\displaystyle\frac{1+\eta }{2}}}. \end{array} \end{equation} \tag{ A.5 }$$

For $\eta =\frac {4}{25}$, one obtains equation (26).