Fractal aggregation in charged granular gases

Fractal aggregation in charged granular gases

Chamkor Singh Max Planck Institute for Dynamics and Self-Organization (MPIDS), 37077, Göttingen, Germany Georg-August-Universität Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany    Marco G. Mazza Max Planck Institute for Dynamics and Self-Organization (MPIDS), 37077, Göttingen, Germany Interdisciplinary Centre for Mathematical Modelling and Department of Mathematical Sciences, Loughborough University, Loughborough, Leicestershire LE11 3TU, United Kingdom
August 6, 2019

The empirical observation of aggregation of small solid particles under the influence of electrostatic forces lies at the origin of the theory of electricity. The growth of clusters formed of small grains underpins a range of phenomena from the early stages of planetesimal formation to aerosol and industrial processes. However, the collective effects of Coulomb forces on the nonequilibrium dynamics and aggregation process in a granular gas, a representative of the above physical systems, are unclear. Here, we establish the premise of a hydrodynamic theory of aggregating granular gases that exchange charges upon collisions and interact via the long-ranged Coulomb forces. We derive the governing equations for the evolution of granular temperature, charge variance, and number density for a homogeneously aggregating gas. We find that, once the aggregates are formed, the system obeys a physical constraint of nearly constant dimensionless ratio of electrostatic energy to granular temperature ; this constraint on the collective evolution of charged clusters is confirmed both by the theory and the detailed molecular dynamics simulations. The inhomogeneous aggregation among monomers and clusters in the mutual electrostatic field statistically proceeds in a self-similar manner. Our work illustrates the defining role of long-ranged interactions and can lead to novel designing principles in particulate systems.

The electrostatic aggregation of small particles is among the oldest scientific observations. Caused by collisional or frictional interactions among grains, large amounts of positive and negative charges can be generated. These clusters have far-reaching consequences: from aerosol formation to nanoparticle stabilization Castellanos (2005); Schwager et al. (2008), planetesimal formation, and the dynamics of the interstellar dust Wesson (1973); Harper et al. (2017); Brilliantov et al. (2015); Blum (2006). The processes accompanying granular collisions and charge buildup can also lead to catastrophic events such as silo failure, or dust explosions.

Experimental investigations of the effects of tribocharging date back to Faraday, and recent in situ investigations have revealed important results Lee et al. (2015); Yoshimatsu et al. (2017); Poppe et al. (2000); Haeberle et al. (2018). However, careful experiments are plagued by a myriad of technical difficulties that often impede unambiguous interpretation Spahn and Sei (2015) of the experiments. A consequence of these difficulties is the lack of consensus about whether electrostatics facilitate or hinder the aggregation process of a large collection of granular particles Spahn and Sei (2015). The aggregation of micron-size particles has been described by mean-field rate equations without considering dissipation or charge transfer mechanisms Ivlev et al. (2002); Dammer and Wolf (2004). Monopolar dilute systems without charge distribution have also been studied Müller and Luding (2011). Intensive studies have been conducted on aggregating systems with short-range potentials Ulrich et al. (2009); Brilliantov et al. (2018).

Despite its far-reaching significance, a statistico-mechanical description of aggregation in a granular system with a mechanism of charge transfer is still lacking. The theoretical treatment requires both dissipation of kinetic energy described by a monotonic dependence of the coefficient of restitution on velocities , and also inclusion of long-range electrostatic forces due to the dynamically-changing charge production. Understanding the growth of charged aggregates requires a statistical approach due to the different kinetic properties and aggregate morphology.

In this work, we derive a Boltzmann equation for the inelastic collisions of grains that interact via Coulomb forces, and produce charges upon collision. We derive the hydrodynamic equations for the number density of the aggregates, the granular temperature , and the charge fluctuations . We find that the dynamics of the charged granular gas approach, but do not overcome, a limiting behavior marked by the value of the dimensionless ratio of electrostatic energy to thermal energy.

To bolster our results, and explicitly consider microscopic fluctuations and morphological structures, we also use three-dimensional molecular dynamics (MD) simulations that explicitly include Coulombic interactions and a charge-exchange mechanism. We find that the granular dynamics agree quantitatively with the predictions of the Boltzmann equation. The cooling gas undergoes a transition from a dissipative to an aggregative phase marked by a crossover in the advective transport.

.1 Kinetics

Figure 1: Schematic of the collision and aggregation process. When particles collide two pathways are considered in the theory depending on their velocities and charges: (i) a typical inelastic collision (dissipation only) with charge exchange, or (ii) an inelastic collision that leads to the aggregation and merging of charges. Primed variables represents post-collision or post-aggregation values.

To obtain a general description, we consider the single particle charge-velocity phase space probability density function , where the particle velocity , charge , and position are the phase space variables, and denotes time. We consider a homogeneously aggregating (homogeneous and quasi-monodisperse) granular gas with charge exchange and dissipation upon collisions. In this limit, the spatial and particle-size dependence of drops out and its time evolution is given by the simplified Boltzmann equation Pitaevskii and Lifshitz (2012); Brilliantov and Pöschel (2010)


where we define as the generalized collision integral which includes dissipation as well as charge exchange during particle collisions. We will now elucidate how the charge exchange mechanism modifies the collisions.

Let us consider contact collisions of particles and with pre-collision velocity-charge values and , respectively. In the ensemble picture, particle collisions will change in the infinitesimal phase-space volumes and , centered around and , respectively. The number of direct collisions per unit volume which lead to loss of particles from the intervals and in time are


where , is the unit vector at collision pointing from the center of particle towards particle , and is the differential collision cross-section. The Heaviside step function selects particles coming towards , while we use to ensure that a contact with an approaching particle takes place only when the Coulomb energy barrier can be overcome, where , F m is the vacuum permittivity, and is the particle diameter. If the interaction is repulsive, is positive, and only if . In case of attractive interaction, is negative and thus always. Essentially, filters repulsive interactions which do not lead to a physical contact between particles.

Consider now particles with initial velocity-charge values and in the intervals . The number of particles per unit volume which, post-collision, enter the interval and in time are

Figure 2: Time evolution of hydrodynamic quantities. Evolution of a. temperature of cluster population , b. charge variance of cluster population , c. number density of cluster population , and d. Bjerrum number (insets are for different filling fraction ). The granular temperature, charge variance and average size of the cluster population during aggregation evolve in such a manner that their non-dimensional combination . This is not captured in the kinetic theory if only restitutive (no aggregation) collisions are considered. The granular MD simulations confirm the analytical results. The number density evolution, however, is highly dynamic and exhibits a non-monotonic behavior due to emergence of macroscopic flow, quantified by (c. inset, see also section: Interplay between fractals and macroscopic flow) in the text


The net change of number of particles in time per unit volume, then reads (see Methods for details)


where and are Jacobians of the transformations for and respectively (see Methods), which lump together the microscopic details of the collision process, namely dissipation and charge exchange in the present study.

Integrating over all incoming particle velocities and charges from all directions, dividing by and taking the limit , we obtain the formal expression for the collision integral

Here we assume that the differential collision cross section, and the contact condition specified by , retain their form for direct and inverse collisions. At this point the particle encounters which do not lead to a physical contact have been excluded using . While taking moments of , a fraction of those contact collisions that lead to aggregation is accounted for by taking the limit for certain conditions on the relative velocity , and by considering the charge transferred to particle equal to the charge on particle [Methods, Eq. (24)-(29)]. In , distant encounters, which do not lead to a contact between particles (glancing collisions) are neglected and the charge exchange and dissipation is considered only during the contact. The long-range effect is incorporated via the collision cross-section [see Eq. (34)].

After setting up the collision integral, we derive the macroscopic changes of number density , granular temperature , and the charge variance (due to aggregation, dissipation and charge exchange, respectively) by taking the moments of the Boltzmann equation (see Methods: splitting restitution and aggregation). The particles are initially neutral and the charge on them is altered either by collisions or during aggregation. However, due to charge conservation during collisions and aggregation, the system remains globally neutral and the mean charge variation is zero. An alternative choice is thus . In order to obtain closed form equations, and for analytical tractability, we assume quasi-monodispersity and homogeneity of the aggregating granular gas at any given time, as illustrated in Fig. 1. This means that during aggregation the mass of the clusters is assumed to grow homogeneously, while their numbers decreases in a given volume. In the following we show that even under these assumptions our hydrodynamic equations provide predictions close to the MD results. We also show that the time evolution of , however, is more involved.

After integration under the assumption of uncorrelated, Gaussian charge and velocity distributions (see Methods) we find the governing equations


where the coefficients , , , , (listed in Table 1) are functions of the time dependent ratio


between electrostatic energy and granular temperature. We term the ratio as Bjerrum number. The variables are also time-dependent functions of aggregating cluster mass , size , and the material constants (Table 1).

First of all, we note the consistency of the above set of equations with the classical Haff’s law in the absence of collisional charging. In this limit , , and we obtain , and , reducing to Haff’s law.

Figure 3: Aggregation process. Snapshots of the aggregating gas at different non-dimensional times (a. , b. , and c. ). Here clusters containing 10 particles or more () are shown, and each color represents a different cluster. Clusters are identified on the basis of the monomer distances: if the centers of two particles are separated by a particle diameter, or less, they belong to the same cluster. Here is the non-dimensional time.

For non-zero charge exchange, the coupled equations (5)-(7) are solved numerically. Figure 2 shows the temporal evolution of the hydrodynamic quantities , , , as well as the dimensionless parameter . Initially, the collisional velocities remain larger than the threshold (see Eq. (23) and (39) in Methods for the role of the threshold ). In this time regime, the collisions are purely restitutive leading to charge exchange and dissipation. The charge exchange increases while the dissipation reduces . The number density , and thus size of the effective particles, remain unchanged as there is no aggregation. The Bjerrum number , however, increases. As approaches the threshold , , and the aggregation initiates, the concentration begins to decrease and thus the size of the aggregates (now considered as new “‘effective” particles; see Fig. 1) increases. The individual particles, or monomers, cluster in such a way that the charge variance on the cluster population now begins to reduce. The theory shows that once aggregation sets in, the aggregating granular gas obeys the constraint


This physical limit is also confirmed by granular MD simulations (see Methods, and Singh and Mazza (2018), for details on MD). At long times, adjusts itself such that the right hand sides of the equations for number density, temperature and charge variance [Eq. (5)-(7)] remain real valued.

To understand the origin of this approximate conservation law in the theory, we solve Eq. (5)-(7) for a system with only dissipative collisions and without aggregation; hence, the cluster size remains unchanged. The results are shown in Fig. 2 (dashed lines). In this limit of only restitutive kinetics, increases continuously, while deviates from Haff’s law. If the complete aggregation dynamics are considered (i.e., dissipation and aggregation) we find that Haff’s law for aggregating clusters is closely obeyed. This is also confirmed by our MD simulations.

Figure 2(b) shows that, after the initiation of aggregation, the charge variance of the aggregating clusters also follows a power law decay. It is notable that a different charge exchange model might provide a different charge buildup rate during the purely dissipative (restitutive) phase, however, the decay of charge variance during the aggregation of the monomers into the clusters is not expected to be influenced by a change in the charge exchange mechanism. The reason is that aggregation sets in at relatively low temperature where the relative thermal motion of monomers inside the clusters is insignificant, and thus collisional charge-exchange is expected to be negligible. In our theory, the charge transferred to a particle during aggregation with particle is equal to the charge on particle as they are considered merged into one single particle (see Fig. 1). Thus the charge variance of cluster population is reduced by the aggregation process, rather than by the collisional charge exchange in the aggregation phase. In the MD simulations, the decay of during the aggregation phase agrees well with the theory, even though we observe that the charge exchange in the simulations leads to a symmetric but non-Gaussian charge distribution during the restitution phase (see Methods, Fig. 5). Although neglected in the current model, the discharging of the monomers inside the aggregates at very low collision speeds is neglected, which might change the charge variance decay rate.

The key finding here is that the decay of the charge variance and the growth of the size of the aggregates is balanced by the decay of temperature during the aggregation, resulting in the stationary value of . The steady state of is apparent in the theory. The granular MD simulations suggest and confirm the upper limit of . We simulate the above results under different initial filling fractions and find that the constraint persists with variation of . The temperature of the cluster population also closely follows Haff’s law for different (see Fig. 2 insets).

The number density evolution from the MD simulations is more involved than the prediction of the quasi-monodisperse kinetic theory. The theory indicates the decay of cluster density only in an average sense (Fig. 2(c) and insets). To explore the mechanism behind the non-monotonic evolution of , we first explore the cluster dynamics and nature of the structures from the full spatially-heterogeneous MD simulations.

coefficient expression
Table 1: Expressions of the coefficients in (5),(6) and (7). The constant is from and (30) while and are from (27), and influence the viscoelastic properties of the particles and the threshold aggregation velocities. Other notations are as described in the main text.

.2 Inhomogeneous aggregation and fractal growth

To gain access to the spatial structure formation in the gas, we perform a detailed cluster analysis of the results from granular MD simulation (see Methods and Singh and Mazza (2018) for more details); our cluster analysis code is made publicly available (see Support Information).

Figure 3 shows snapshots of the clusters at different times from a typical run.

Figure 4: Dynamics of fractal dimension. a-c. The scaling between cluster mass and their radius of gyration , at different times, and d. the time evolution of , thus obtained, for different filling fractions. The average fractal dimension in the aggregating charged gas varies across reported average values for ballistic cluster-cluster aggregation (BCCA, ) and diffusion-limited particle-cluster aggregation (DLPCA, ) Blum (2006); Smirnov (1990). .

The morphology of the aggregates is studied by computing the average fractal dimension Mandelbrot (1977); Jullien (1987) for the cluster population from the scaling law between cluster mass and radius of gyration , where the index runs over all monomers in a given aggregate. Figure 4(a-c) shows scatter plots for versus at different times. Figure 4(d) shows the time evolution of the exponent for varying filling fraction. The average fractal dimension lies in between the average values reported for the ballistic cluster-cluster aggregation (BCCA, ) and the diffusion-limited particle-cluster aggregation (DLPCA, ) Blum (2006); Smirnov (1990) models. In time, is dynamic and changes across the two model limits. These results indicate that the aggregate structures retain their fractal nature over time.

The BCCA and DLPCA mechanisms of aggregation have been reported in neutral dust agglomeration (e.g hit and stick, ballistic motion, Blum (2006)), wet granulate aggregation (sticking due to capillary bridges and ballistic motion, Ulrich et al. (2009)), colloidal aggregation (van der Waals and repulsions Lebovka (2012)), and hit and stick agglomeration in Brownian particles under frictional drag Kempf et al. (1999). The observation that lies between the reported average values of for BCCA and DLPCA indicates the presence of mixed characteristics from both of these simplified models. The size distribution in an aggregating charged granular gas (see Singh and Mazza (2018)) tends to resemble a DLPCA-like behavior where the smaller size aggregates are larger in number, in contrast to a BCCA-like model where the size distribution is typically bell-shaped Blum (2006). On the other hand, the monomer motion is found to be highly sub-diffusive Singh and Mazza (2018) in agreement with the BCCA model. In addition, the Coulombic interactions will cause considerable deviations from the short-ranged or ballistic propagation typical of the BCCA or DLPCA models. The long-range forces due to a bipolar charge distribution lead to the intermediate value of between what is observed in the above two physical aggregation models.

.3 Interplay between fractals and macroscopic flow

Apart from the long-range effects, the morphology of the aggregates is also altered by additional mechanisms. We discuss two physical processes that are not captured in the analytical theory, but that we investigate via our MD simulations.

First, in our Boltzmann kinetic description, the collisions between aggregates at any given time are considered as collisions between two spheres with sizes equal to the average size of the aggregate population. This is a quasi-monodisperse assumption typically used in cluster-cluster aggregation (CCA) models. The quasi-monodispersity however neglects the morphology and surface roughness of the colliding aggregates. Collisions between aggregates with large size difference, and with individual monomers, are also simplified.

Secondly, granular gases are characterized by the emergence of a convective flow Hummel et al. (2016) which we find, in the present case, inducing the non-monotonicity in the temporal evolution of the number density [Fig. 2(c)]. To quantify the emergence of the macroscopic flow, we calculate the Mach number


Here the system is divided into equal sized cubic boxes and is the center of mass velocity in a box (or local advective velocity), scaled by , the thermal velocity of the monomers in the gas. The time evolution of is shown in Figure 2(c) (inset), which indicates that grows at a higher rates when the number density evolution becomes non-monotonic (), indicating a coupling between aggregation and the generation of macroscopic flow. Due to the macroscopic flow (see also the Supplementary Media), aggregates which are weakly connected are prone to fragmentation. This results in an intermediate regime where the concentration increases instead of decreasing.

Excluding the two mechanisms above explains the shortcoming of our quasi-monodisperse Boltzmann theory in capturing the non-monotonic behavior of occurring at , and only exhibits the general downward trend of the number density [Fig 2(c)].

.4 Conclusions

We have derived the rate equations for the evolution of the number density , granular temperature , and charge variance of the cluster population in a charged, aggregating granular gas. In contrast to well-known Smoluchowski-type equations, we have explicitly coupled to the decay of and charge variance. We have compared the results with three-dimensional molecular dynamics simulations and the outcomes of a detailed cluster analysis, and have explored the morphology of the aggregating structures via fractal dimension.

Taken together, our results indicate that the aggregation process in a charged granular gas is quite dynamic, while respecting some physical constraints. The growth process obeys , while morphologically, the clusters exhibit statistical self-similarity, persistent over time during the growth. The fractal dimension and growth of structures is intermediate between the BCCA and DLPCA models. We also demonstrate that the application of a purely dissipative kinetic treatment is not sufficient to make predictions about global observables such as and in an aggregating charged granular gas.

Finally, we believe that our kinetic approach can be applied to study aggregation processes in systems such as wet granulate collections with ion transfer mechanism Lee et al. (2018); Zhang et al. (2015), dissipative cell or active particle collections under long-range hydrodynamic and electrostatic effects Yan et al. (2016); Friedl and Gilmour (2009), and charged ice-ice collisions Dash et al. (2001).

.5 Acknowledgements

We gratefully acknowledge highly constructive comments by Jürgen Blum, Eberhard Bodenschatz, Stephan Herminghaus, and Reiner Kree. We thank the Max Planck Society for funding.

.6 Methods and supporting information

Kinetics and generalized collision integral. Let us consider the number of particles per unit volume having initial velocity-charge values and in the intervals and which, post-collision, enter the in the intervals and in time are


and thus the net increase of number of particles per unit time and volume is . We can relate the primed variables to the unprimed via


where is the Jacobian of the transformation for viscoelastic particles. Here is a material constant. To obtain the transformation , we consider the ratio of relative charges after and before the collision


and in addition we impose charge conservation during collisions


The above two relations finally provide the transformation


where, for example, for a constant . This means that the differential charge-space volume element shrinks or expands by a factor of . In general, for velocity and charge-dependent case the expressions of and can be quite complicated as it depends on how the charge exchange takes place during collisions and its dependence on myriad factors (such as size, composition, and crystalline properties). Incorporating the above phase-space volume transformations due to collisions, the net change of number of particles per unit phase-space volume and in time reads


where we assume that the differential cross-section and the contact condition specified by are the same for direct and inverse collisions. Finally, dividing by , and integrating over all incoming particle velocities and charges from all directions in the limit , we obtain a formal expression for the collision integral


At this point the particle encounters which do not lead to a physical contact have been excluded using , however, collisions that lead to aggregation have not been explicitly accounted. We do so by taking the limit for certain conditions on the relative velocity , and by considering the charge transferred to particle equal to the charge on particle ((24)-(29)). In , distant encounters, which do not lead to a contact between particles (glancing collisions) are neglected and the charge exchange and dissipation is considered only during the contact. The long-range effect is incorporated via collision cross-section in Eq. (34).

Splitting restitution and aggregation.

The time rate of change of the average of a microscopic quantity is obtained by multiplying the Boltzmann equation for by and integrating over , i.e.


It can be shown that


where and is the change of during the collision between pair , and the prime denotes a post collision value. We note that the transformations in Eq. (12) and (15) are reversed back while integrating weighted with quantity of interest . We consider the number density, the kinetic energy or granular temperature, and the charge variance (the system is globally neutral and the mean charge variation is zero), respectively


At this point we differentiate the restitutive or dissipative collisions from aggregative ones by splitting as


where . If , the particles collide and separate after the collision irrespective of the sign of (attractive or repulsive). This leads to dissipation of energy with finite non-zero , and charge exchange according to a specified rule. The aggregative part is zero in this case. If and (attractive encounters at low velocities), the particles collide and aggregate with , and with charge exchange to particle equal to . If and (repulsive encounters at low velocities), no physical contact takes place between the particles which leads to neither dissipation nor aggregation ().

The expressions for and are obtained as follows. The particle number does not change during a dissipative collision but reduces by one in an aggregative collision, i.e.


For the granular temperature,


where we take the limit for the aggregation.

The change in the charge variance is obtained as


where equals the charge transferred to particle during its collision with particle , and is the mean charge on the pair. For the charge transfer, we consider


which together with Eq. (26), gives


For aggregation, the charge tranferred to particle equals the charge on the merging particle , i.e, , which gives


Putting Eq. (24),(25),(28), and (29) in (23), and then (23) in (19), the resulting integrals are solved (see Appendix for details), assuming the statistical independence of charge-velocity distribution function, i.e., , and assuming that they are Gaussian. In addition to the charge exchange, the coefficient of restitution is taken as velocity dependent, i.e


while the long-range effects due to Coulomb interactions are taken into account by the change in collision cross section.

Granular MD simulations. The equation of motion of the form


is solved for a setup with periodic boundary conditions in a cubic box of size , where is monomer diameter. Here the elastic force strength is , and a dissipative force strength of is set to obtain clustering in a finite size system (See Singh and Mazza (2018) for more details). Also . The Coulomb force strength is varied between . The long-range Coulomb forces for a setup with periodic boundary conditions is challenging and conditionally convergent as it depends on the order of summation. We employ the Ewald summation that converges rapidly, and has a computational complexity . The algorithm is partially parallelized and highly optimized on graphics processing unit (GPU). In our simulations, the total computing time to reach simulation time for a typical simulation with , including the long-range electrostatic forces, is of the order of weeks. See Singh and Mazza (2018) for more details.

Codes and media. The following computer codes and media are made public:

  1. MATLAB program to solve Eq. (5)-(7).

  2. MATHEMATICA framework to solve the restitutive and aggregative parts of the kinetic integrals.

  3. The cluster analysis code in MATLAB to obtain the fractal dimension, cluster size distribution, average cluster size, and other statistical quantities, is made available at

  4. Movie of the evolving and aggregating structures in the charged granular gas.

Data availability. All data are available from the authors.

Dimensionless quantities to laboratory values. Typically non-Brownian growth in planetary dust is dominant for monomer sizes of several Zsom et al. (2010) and the growth barrier problem Spahn and Sei (2015) starts to arise near . The mass of silica particles in this range of size is . Thus setting gives a velocity reference value . This fixes the time reference . Thus in our results the growth over units of non-dimensional time approximately implies growth over , or for such size and mass particles. The average size in the growth period grows approximately by one order of magnitude (e.g. the growth of the largest cluster is from to in time under these reference values, and for initial filling fraction ).

.7 Restitutive and aggregative parts of (19).

To solve the integrals in Eq. (19) for different from (20)-(22), we assume that the normalized velocity as well as charge distribution in the gas at any time remain Gaussian, and the two are uncorrelated, i.e


In Figure 5, we show charge distribution of monomers obtained from typical simulation runs. In (32), we assume that shape of the scaled charge distribution remains close to Gaussian during size growth.

Figure 5: The scaled charge distribution of individual particles obtained from typical MD simulation runs (dots) in the aggregated granular gas. The solid line is a Gaussian fit. Here .

The coefficient of restitution is velocity dependent, i.e


The attractive or repulsive long-range effects are emulated through an effective differential cross-section for a binary collision, which changes depending upon the sign and magnitude of charges on the particle pair and their relative velocity, according to


where is the differential cross-section per unit solid angle . The expression is independent of , and thus the total cross section is . For neutral particles, , and thus . For (repulsive encounters), , while for (attractive encounters), .

Below we explain the solution procedure for the restitutive part of the equation for . Similar procedure can then followed for the aggregative part of , except the fact that in this case and integration limits are selected by , where the shorthand . Plugging Eq. (25) into (23) and then Eq. (23) into (19), we find