Institute of Physics, University of Potsdam, Am Neuen Palais, D-14469 Potsdam, Germany Institute for Plasma Research, University of Maryland, College Park, Maryland 20742
Olaf Petzschmann, Udo Schwarz, Frank Spahn, Celso Grebogi, and Jürgen Kurths
March 6, 2001
Clustering of granular gases is caused by the dissipation occurring in inter-particle collisions. We use the generalized dimension to characterize the structural evolution and the index-of-dispersion in order to quantify the typical length scale of the density clusters of force-free granular assemblies obtained by 3-D numerical computer simulations. We derive then an expression for the expected length-scales of cluster formation from stability analysis in the hydrodynamical approximation. We find a good agreement between our theoretical prediction and numerical experiments.
74.80.-g, 47.20.-k, 05.20.Dd, 02.50.Fz
In recent studies much attention has been devoted to the dynamics of granular gases [1, 2, 3, 4]. The difference of such systems to conventional gases is their dissipative nature, i.e., the energy loss occurring in inter-particle collisions. This dissipation causes the granules to form clusters, which results from an increase in the collision frequency in regions of fluctuative density enhancements. As more and more mechanical energy is dissipated, there is a corresponding increase in the effective cooling and a decrease in pressure. The consequence is that the granular matter migrates into regions of smaller pressure and higher densities, resulting in the cluster formation [4, 5].
In this Letter, we argue and demonstrate that the decisive quantity for the rate of granular clustering is the coefficient of restitution , which is the ratio between the normal components of the rebound velocity and the impact velocity of the colliding particles. In general, is a function of the impact velocity . Intuitively, a lower , corresponding to larger dissipation per collision, leads to a more intense clustering of the granular gas. We validate this notion by numerical and analytical investigations of the cooling process in three-dimensional configuration space. We characterize and quantify the clustering of granular assemblies through a statistical analysis of the phase-space data obtained from event-driven particle simulations. We introduce two statistical measures to accomplish that:
(i) To characterize the increase of inhomogeneity in the particle density with time, we define a quantity in terms of the generalized dimensions. This quantity reflects the -dependence of clustering formation. We find a monotonous dependence of the the rate of clustering on the dissipation.
(ii) To quantify the characteristic size of inhomogeneities, we define another quantity in terms of the index-of-dispersion. We find that the dependence of cluster size on , as obtained from our 3-D simulations, agrees well with our theoretical predictions derived from a stability analysis using a hydrodynamical approximation [e.g. Eq. 6]. =
Figure: Two-dimensional projection of the 3-D configuration space slice obtained after about 15 collisions per particle, where , the particle diameter is 2 cm and cm.
The reason for cluster formation is that the fluctuations of the homogeneous cooling state (HCS) are unstable . In this case, the molecular chaos assumption breaks down, the assumption which serves as the foundation for the formulation of a kinetic theory of Boltzmann-type. Numerical simulations have shown  that the beginning of the cooling process is well described by the HCS, while, after a certain critical time, depending on the number of collisions per particle, the cooling process is slowed down by velocity fluctuations. After this critical moment, the cooling is almost no longer determined by particle-particle but by cluster-cluster interactions. This would imply that the clusters have already been formed during this HCS-phase.
In our simulations, we use a 3-D event-driven code with identical spherical particles of diameter . Using periodic boundary conditions, the particle ensemble is simulated in a cubic box of side length , chosen to correspond to a filling factor of about 0.02. Initially, the particles are randomly placed in the configuration space having initial velocities with randomly chosen directions according to an initial distribution function . We ensure that the initial conditions are homogeneous in density and isotropic in the velocities. For convenience, we restrict ourselves to the case of constant coefficients of restitution fixing its values at 0.1, 0.2, ..., 0.9. In our simulations only translational degrees of freedom are considered and, thus, only the normal coefficient of restitution is taken into account. According to this, the relative velocity between the collision partners is projected onto the direction of the particle centers, represented by the unit vector . The relative velocity after the collision is then given by . We simulate the assemblies of different for the same absolute time, taking care that the system is still in the HCS-phase.
Compared to the 2-D case (see e.g. Fig. 1 in or Fig. 3
in), a 3-D granular assembly shows a much
weaker clustering, when we project, for visualization purpose the 3-D
particle locations onto a
two-dimensional plane, as illustrated in Fig. 1.
Thus, the inhomogeneities in Fig. 1 are
difficult to distinguish from those ones of an ensemble with purely
elastic particle interactions.
Therefore, we use measures from nonlinear dynamics and
statistics to characterize and quantify the spatial inhomogeneities
in our simulations.
First, we split the box in m = 2, 4, 8, ..., 4096 sub-boxes.
The equal length L of the diagonal of each sub-box is given by
in units of the smallest grid
size , with .
Second, we characterize the cluster formation quantitatively, using
the generalized dimension, based on the Rényi entropies
 defined as
the probability for
particles to be in the box labeled by
Recently, we have shown that the generalized dimension , defined
in terms of , is an appropriate tool to
characterize the time evolution of the degree of clustering
[5, 8]. The
generalized dimension is defined by
for each time step t and restitution coefficient .
As q increases from zero, accentuates more and more the degree
of clumping (high density regions) while, as q decreases from zero,
accentuates better the degree of depletion (low density regions).
To characterize the increase of
inhomogeneity with time, we calculate from simulation data
how changes with
q, i.e., .
Then, to get a representative measure of the change of inhomogeneity
with time, we introduce the quantity (for )
where and are the times at the beginning and at the end of the simulation, respectively. The value , based on the change of the generalized dimension with q reflects the -dependent structure formation. Figure 2 shows the dependence of on indicating that an initially homogeneous distribution of particles changes the greatest the smaller the coefficient of restitution is. For elastically colliding granules (), one gets corresponding to an unchanged homogeneous particle distribution during the time evolution. =
Figure: vs. [Eq. (1)]. The diamonds are the simulation data points, the solid line is the result of a quadratic regression fit.
This behavior agrees with predictions [4, 5] that, as the dissipation rate increases, the clustering becomes more intense. However, it is not possible to quantify, using this method, the length-scales of the inhomogeneities in the particle density of the simulations. Therefore, we introduce a new measure defined in terms of the index-of-dispersion . The test statistic of the index-of-dispersion is given by with the particle density in the box, the standard deviation , and the number of degrees of freedom m (total number of boxes). This test statistic is that of a goodness-of-fit test for the hypothesis that the particles are independently and uniformly distributed in the whole simulation box. =
Figure: The evolution of the index-of-dispersion [Eq. (2)] in the critical region vs. the sub-box diagonal L and the simulation time t for different restitution coefficients . The value of is color coded as indicated in the color bar. Black indicates a Poissonian distribution (), values higher than 0 (colored) indicate a significant lumpy distribution of the particles in the configuration space. White corresponds to the highest deviation from a Poissonian particle distribution.
For a homogeneous particle distribution, the index-of-dispersion satisfies
which is characteristic for a Poissonian distribution
A higher index-of-dispersion points towards
inhomogeneities, while is typical for a
regular distribution of particles, as it is found, e.g., in lattices.
Therefore, we define the normalized quantity
where we relate the index-of-dispersion to the critical value . The significance level of the -statistics is chosen to be . Hence, is a measure of the inhomogeneity, i.e., of the degree of clustering. The Heaviside function ensures that purely regular features are excluded from the statistical treatment. The time evolution of in the critical region () for different values of , obtained from the 3-D simulations, is shown in Fig. 3. Here the quantity is color coded and plotted vs. time t and length scale L. The clumping starts at the smallest scales L after a characteristic time. At later times, inhomogeneities at larger scales evolve. Figure 3 shows the tendency of the values as a function of L, and time. Due to the finite number of particles used in the simulations, one observes, in addition to the density inhomogeneities caused by the dissipation, fluctuations .
The most prominent length-scale of the
inhomogeneities, under these conditions, is obtained as follows.
From the index-of-dispersion, which is in
the critical region of the corresponding -test statistics
(Fig. 3), we find the following expression for the critical length
scale by averaging the value L of space and time, as
Note that the size of the inhomogeneities detected by this procedure [Eqs. (2)-(3)] are half of the size of the typical wave length. Using Eq. (3), we derive the typical wave number based on the numerical simulations. Figure 4 shows , Eqs. (3), as a function of for a finite simulation time.
To interpret and better understand the results of our computer
we use a stability analysis in the hydrodynamical
approximation, which describes well the physics of a granular gas
during the HCS-phase [1, 2, 3, 4, 5, 6].
The equations of hydrodynamics for a force-free granular gas read
[10, 11] as
where the material derivative is given by ; ":" and "" denote the double inner product and the dyadic product, respectively. The quantities , , and T are the density, the average particle velocity, and the temperature, respectively. The heat flux is , and the pressure tensor is denoted by with the shear tensor , where p is the pressure and is the unity tensor. The constitutive relations in the dilute limit, based on an expansion of the distribution function around the Maxwell-Boltzmann distribution, read for a constant coefficient of restitution  for the viscosity, for the heat conductivity, and for the cooling due to inelastic collisions. Here, is the particle diameter and the mass and the Boltzmann constant are .
In order to analyze the stability of an initially homogeneous granular
gas, we linearize the hydrodynamic Eqs. (4) and
expand the state variables about the ground state
with . For a force-free granular gas,
is constant in time and space, , and the ground state temperature is
a function of time, due to the permanent cooling.
In a periodic or infinite system, the eigenfunctions of the linearized
Eqs. (4) are plane waves , where is the wave vector and is the growth rate. Without loss of
generality, one can
choose . The resulting eigenvalue problem
gives a system of linear equations , where is the
with the solubility condition , determining a characteristic fifth order polynomial,
which can easily be reduced to third order.
If one of its roots takes , an
instability is found.
The linear stability analysis yields a condition, for which an
initially homogeneous granular ensemble becomes monotonically
The subscript 0 denotes the ground state about which the linear stability analysis is carried out. To quantify the derivatives in Eq. (5), we apply the equation of state of an ideal gas with as an independent variable . Using this, we get , , , and . This yields an expression for the critical wave number . According to this equation, unstable wave-like oscillations with a characteristic length scale only occur if , which quantifies the pressure instability mentioned above. It has been shown that clusters are formed under this condition either for a constant  or variable  .
In the dilute
limit, the temperature evolution of an initial homogeneous
granular gas with
and an initial temperature is given by
with , and therefore,
, which, together with
the constitutive relations for and
yields the critical wave number
Equation (6) indicates that clustering sets in after a certain time , which is the order of the average collision time for . On the other hand, for a given finite time t, there is a characteristic coefficient of restitution , where has a pole and above which no clustering is expected. Obviously, there is no clustering for . Especially, for () one gets [dotted line in Fig. 4], which is the result found by Goldhirsch & Zanetti , where the derivative has been assumed to be a negative constant. Note that the term has almost no influence on the shape of . Figure 4 shows given by Eq. (6) as a function of for a finite time and for .
In Fig. 4, the wave number derived from numerical simulations (stars and solid line for the fit) is compared with the corresponding theoretical relation for [Eq. (6)] in dependence on (dashed line). We find a surprisingly good agreement between the theoretical approach and the numerical experiments. This confirms the validity of the theory to explain cluster-formation, including the sensitive quantitative properties such as the typical cluster-size. The rather small deviation of the theoretical curve from the data points is due to fluctuations associated with the finite number of particles, and the averaging procedure. The dotted line corresponds to the wave number of the homogeneous cooling for , indicating that no dramatic change of cluster formation can be expected under homogeneous cooling state conditions.
Figure: Critical wave number vs. . The stars are the data points obtained from the simulations, the solid line shows the curve for the quadratic regression fit. The dashed line presents the result of the stability analysis [Eq. (6)] using the same parameters as in the simulations, while the dotted one is the result for .
In conclusion, our theoretical investigations, based on a stability analysis of the hydrodynamical equations, Eqs. (4), yield typical cluster sizes depending on the dissipation, represented by the wave number [Eq. (6)]. Using the generalized dimension and the index-of-dispersion, we characterize and quantify the structural evolution of density inhomogeneities. Statistical analyses of 3-D numerical simulation data, based on the index-of-dispersion, yield typical length-scales of cluster formation [Eq. (3)] as a function of the restitution , which is in good agreement with our theoretical predictions [Eq. (6)]. However, the homogeneous cooling state is a transient phase for the evolution of an initially homogeneous force-free granular gas. The linearized hydrodynamic approximation, as well as the constitutive relations, might not hold for granular assemblies with very inhomogeneous particle distributions and for . Therefore, the analysis of the cooling regime beyond the homogeneous cooling state, both through theoretical investigations and numerical simulations, is needed.
This work was supported by the "Deutsche Forschungsgemeinschaft" (grant Sp 384/5-1), the "Deutscher Akademischer Austauschdienst", and the "Alexander von Humboldt Stiftung".