next up previous
Next: References

tex2html_wrap_inline372 Institute of Physics, University of Potsdam, Am Neuen Palais, D-14469 Potsdam, Germany tex2html_wrap_inline384 Institute for Plasma Research, University of Maryland, College Park, Maryland 20742

Length-scales of clustering in granular gases

Olaf Petzschmanntex2html_wrap_inline372, Udo Schwarztex2html_wrap_inline372, Frank Spahntex2html_wrap_inline372, Celso Grebogitex2html_wrap_inline378, and Jürgen Kurthstex2html_wrap_inline372

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 tex2html_wrap_inline386, which is the ratio between the normal components of the rebound velocity and the impact velocity of the colliding particles. In general, tex2html_wrap_inline386 is a function of the impact velocity tex2html_wrap_inline390. Intuitively, a lower tex2html_wrap_inline386, 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 tex2html_wrap_inline386-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 tex2html_wrap_inline386, 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 tex2html_wrap_inline398, the particle diameter is 2 cm and tex2html_wrap_inline402 cm.

The reason for cluster formation is that the fluctuations of the homogeneous cooling state (HCS) are unstable [6]. 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 [6] 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 tex2html_wrap_inline404 identical spherical particles of diameter tex2html_wrap_inline406. Using periodic boundary conditions, the particle ensemble is simulated in a cubic box of side length tex2html_wrap_inline408, chosen to correspond to a filling factor of about 0.02. Initially, the particles are randomly placed in the configuration space having initial velocities tex2html_wrap_inline412 with randomly chosen directions according to an initial distribution function tex2html_wrap_inline414. 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 tex2html_wrap_inline418 between the collision partners is projected onto the direction of the particle centers, represented by the unit vector tex2html_wrap_inline420. The relative velocity after the collision is then given by tex2html_wrap_inline422. We simulate the assemblies of different tex2html_wrap_inline386 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[5] or Fig. 3 in[4]), 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 tex2html_wrap_inline430 in units of the smallest grid size tex2html_wrap_inline432, with tex2html_wrap_inline434. Second, we characterize the cluster formation quantitatively, using the generalized dimension, based on the Rényi entropies [7] defined as tex2html_wrap_inline436 with tex2html_wrap_inline438, the probability tex2html_wrap_inline440 for particles to be in the tex2html_wrap_inline442 box labeled by tex2html_wrap_inline444. Recently, we have shown that the generalized dimension tex2html_wrap_inline446, defined in terms of tex2html_wrap_inline448, is an appropriate tool to characterize the time evolution of the degree of clustering [5, 8]. The generalized dimension is defined by tex2html_wrap_inline450 for each time step t and restitution coefficient tex2html_wrap_inline386. As q increases from zero, tex2html_wrap_inline446 accentuates more and more the degree of clumping (high density regions) while, as q decreases from zero, tex2html_wrap_inline446 accentuates better the degree of depletion (low density regions). To characterize the increase of inhomogeneity with time, we calculate from simulation data how tex2html_wrap_inline446 changes with q, i.e., tex2html_wrap_inline468. Then, to get a representative measure of the change of inhomogeneity with time, we introduce the quantity (for tex2html_wrap_inline470)
where tex2html_wrap_inline472 and tex2html_wrap_inline474 are the times at the beginning and at the end of the simulation, respectively. The value tex2html_wrap_inline476, based on the change of the generalized dimension with q reflects the tex2html_wrap_inline386-dependent structure formation. Figure 2 shows the dependence of tex2html_wrap_inline476 on tex2html_wrap_inline386 indicating that an initially homogeneous distribution of particles changes the greatest the smaller the coefficient of restitution tex2html_wrap_inline386 is. For elastically colliding granules (tex2html_wrap_inline488), one gets tex2html_wrap_inline490 corresponding to an unchanged homogeneous particle distribution during the time evolution. =

Figure:   tex2html_wrap_inline476 vs. tex2html_wrap_inline386 [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 tex2html_wrap_inline496 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 [9]. The test statistic of the index-of-dispersion is given by tex2html_wrap_inline498 with the particle density tex2html_wrap_inline500 in the tex2html_wrap_inline442 box, the standard deviation tex2html_wrap_inline504, and the number of degrees of freedom m (total number of boxes). This test statistic is that of a tex2html_wrap_inline508 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 tex2html_wrap_inline510 vs. the sub-box diagonal L and the simulation time t for different restitution coefficients tex2html_wrap_inline386. The value of tex2html_wrap_inline518 is color coded as indicated in the color bar. Black indicates a Poissonian distribution (tex2html_wrap_inline520), 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 tex2html_wrap_inline522, which is characteristic for a Poissonian distribution of particles. A higher index-of-dispersion tex2html_wrap_inline524 points towards inhomogeneities, while tex2html_wrap_inline526 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 tex2html_wrap_inline528. The significance level of the tex2html_wrap_inline508-statistics is chosen to be tex2html_wrap_inline532. Hence, tex2html_wrap_inline518 is a measure of the inhomogeneity, i.e., of the degree of clustering. The Heaviside function tex2html_wrap_inline536 ensures that purely regular features are excluded from the statistical treatment. The time evolution of tex2html_wrap_inline518 in the critical region (tex2html_wrap_inline524) for different values of tex2html_wrap_inline386, obtained from the 3-D simulations, is shown in Fig. 3. Here the quantity tex2html_wrap_inline518 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 tex2html_wrap_inline518 as a function of L, tex2html_wrap_inline386 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 tex2html_wrap_inline558.

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 tex2html_wrap_inline508-test statistics (Fig. 3), we find the following expression for the critical length scale tex2html_wrap_inline562 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 tex2html_wrap_inline566 based on the numerical simulations. Figure 4 shows tex2html_wrap_inline568, Eqs. (3), as a function of tex2html_wrap_inline386 for a finite simulation time.

To interpret and better understand the results of our computer simulations, 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 tex2html_wrap_inline572; ":" and "tex2html_wrap_inline574" denote the double inner product and the dyadic product, respectively. The quantities tex2html_wrap_inline576, tex2html_wrap_inline578, and T are the density, the average particle velocity, and the temperature, respectively. The heat flux is tex2html_wrap_inline582, and the pressure tensor is denoted by tex2html_wrap_inline584 with the shear tensor tex2html_wrap_inline586, where p is the pressure and tex2html_wrap_inline590 is the unity tensor. The constitutive relations in the dilute limit, based on an expansion of the distribution function tex2html_wrap_inline592 around the Maxwell-Boltzmann distribution, read for a constant coefficient of restitution [2] tex2html_wrap_inline594 for the viscosity, tex2html_wrap_inline596 for the heat conductivity, and tex2html_wrap_inline598 for the cooling due to inelastic collisions. Here, tex2html_wrap_inline406 is the particle diameter and the mass and the Boltzmann constant are tex2html_wrap_inline602.

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 tex2html_wrap_inline604 with tex2html_wrap_inline606. For a force-free granular gas, tex2html_wrap_inline608 is constant in time and space, tex2html_wrap_inline610, and the ground state temperature tex2html_wrap_inline612 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 tex2html_wrap_inline614, where tex2html_wrap_inline616 is the wave vector and tex2html_wrap_inline618 is the growth rate. Without loss of generality, one can choose tex2html_wrap_inline620. The resulting eigenvalue problem gives a system of linear equations tex2html_wrap_inline622, where tex2html_wrap_inline624 is the coefficient matrix with the solubility condition tex2html_wrap_inline626, determining a characteristic fifth order polynomial, which can easily be reduced to third order. If one of its roots takes tex2html_wrap_inline628, an instability is found. The linear stability analysis yields a condition, for which an initially homogeneous granular ensemble becomes monotonically unstable [5],
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 tex2html_wrap_inline576 as an independent variable tex2html_wrap_inline634. Using this, we get tex2html_wrap_inline636, tex2html_wrap_inline638, tex2html_wrap_inline640, and tex2html_wrap_inline642. This yields an expression for the critical wave number tex2html_wrap_inline644. According to this equation, unstable wave-like oscillations with a characteristic length scale tex2html_wrap_inline646 only occur if tex2html_wrap_inline648, which quantifies the pressure instability mentioned above. It has been shown that clusters are formed under this condition either for a constant [4] or variable [5] tex2html_wrap_inline386.

In the dilute limit, the temperature evolution of an initial homogeneous granular gas with tex2html_wrap_inline652 and an initial temperature tex2html_wrap_inline654 is given by tex2html_wrap_inline656 with tex2html_wrap_inline658, and therefore, tex2html_wrap_inline660, which, together with the constitutive relations for tex2html_wrap_inline662 and tex2html_wrap_inline664 yields the critical wave number
Equation (6) indicates that clustering sets in after a certain time tex2html_wrap_inline666, which is the order of the average collision time for tex2html_wrap_inline668. On the other hand, for a given finite time t, there is a characteristic coefficient of restitution tex2html_wrap_inline386, where tex2html_wrap_inline674 has a pole and above which no clustering is expected. Obviously, there is no clustering for tex2html_wrap_inline488. Especially, for tex2html_wrap_inline678 (tex2html_wrap_inline680) one gets tex2html_wrap_inline682 [dotted line in Fig. 4], which is the result found by Goldhirsch & Zanetti [4], where the derivative tex2html_wrap_inline684 has been assumed to be a negative constant. Note that the term tex2html_wrap_inline686 has almost no influence on the shape of tex2html_wrap_inline568. Figure 4 shows tex2html_wrap_inline568 given by Eq. (6) as a function of tex2html_wrap_inline386 for a finite time and for tex2html_wrap_inline678.

In Fig. 4, the wave number tex2html_wrap_inline568 derived from numerical simulations (stars and solid line for the fit) is compared with the corresponding theoretical relation for tex2html_wrap_inline568 [Eq. (6)] in dependence on tex2html_wrap_inline386 (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 tex2html_wrap_inline568 of the homogeneous cooling for tex2html_wrap_inline678, indicating that no dramatic change of cluster formation can be expected under homogeneous cooling state conditions.


Figure:   Critical wave number tex2html_wrap_inline568 vs. tex2html_wrap_inline386. 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 tex2html_wrap_inline678.

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 tex2html_wrap_inline568 [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 tex2html_wrap_inline386, 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 tex2html_wrap_inline716. 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".

next up previous
Next: References

Udo Schwarz
Tue Mar 6 08:27:03 MET 2001