Dilute Wet Granulates: Nonequilibrium Dynamics and Structure Formation

Dilute Wet Granulates: Nonequilibrium Dynamics and Structure Formation

Stephan Ulrich ulrich@theorie.physik.uni-goettingen.de Universität Göttingen, Institute of Theoretical Physics, Germany    Timo Aspelmeier Max-Planck-Institut für Dynamik und Selbstorganisation, Dept. Dynamics of Complex Fluids, Göttingen, Germany    Annette Zippelius Universität Göttingen, Institute of Theoretical Physics, Germany Max-Planck-Institut für Dynamik und Selbstorganisation, Dept. Dynamics of Complex Fluids, Göttingen, Germany    Klaus Roeller Max-Planck-Institut für Dynamik und Selbstorganisation, Dept. Dynamics of Complex Fluids, Göttingen, Germany    Axel Fingerle Max-Planck-Institut für Dynamik und Selbstorganisation, Dept. Dynamics of Complex Fluids, Göttingen, Germany    Stephan Herminghaus Max-Planck-Institut für Dynamik und Selbstorganisation, Dept. Dynamics of Complex Fluids, Göttingen, Germany
July 5, 2019

We investigate a gas of wet granular particles, covered by a thin liquid film. The dynamic evolution is governed by two-particle interactions, which are mainly due to interfacial forces in contrast to dry granular gases. When two wet grains collide, a capillary bridge is formed and stays intact up to a certain distance of withdrawal when the bridge ruptures, dissipating a fixed amount of energy. A freely cooling system is shown to undergo a nonequillibrium dynamic phase transition from a state with mainly single particles and fast cooling to a state with growing aggregates, such that bridge rupture becomes a rare event and cooling is slow. In the early stage of cluster growth, aggregation is a self-similar process with a fractal dimension of the aggregates approximately equal to . At later times, a percolating cluster is observed which ultimately absorbs all the particles. The final cluster is compact on large length scales, but fractal with on small length scales.

45.70.-n, 47.57.-s, 61.43.Hv
preprint: APS/123-QED

I Introduction

Granular materials are systems of macroscopic particles which interact only when they are in mutual contact, and the interaction is dissipative. In spite of this simple definition, collective phenomena arising in such sytems are of utmost complexity, and have inspired strongly increasing research activities in recent years. The particular interest in granular systems is mainly due to the fact that their importance spans from technology and applied research to very fundamental questions of interdisciplinary relevance. On the one hand, storage and handling of bulk solids is among the most significant tasks in industrial technology, and still poses a large number of unsolved problems Jaeger et al. (1996); de Gennes (1999); Duran (2000). On the other hand, granular systems provide a comparatively simple, experimentally accessible model for physics far from equilibrium Brilliantov and Pöschel (2004); Kudrolli (2004); Umbanhowar et al. (1996). This is at the heart of self-organization and pattern formation processes, so that granular systems have been considered as genuine model systems for structure formation on various length scales, including the formation of planetesimals from interstellar dust and the formation of planets and stars from accretion discs ast ().

In most studies so far, models were inspired by dry granular systems, where the dissipative contact interaction consists in the loss of a certain fraction of the kinetic energy in every impact. Adding a small amount of liquid to the granular system changes its properties dramatically: while dry sand can flow freely similar to a liquid, wet sand has properties of a plastic solid. This difference in the macroscopic behavior is reflected in a corresponding difference in particle interactions Herminghaus (2005). The collisions of dry granulates are typically purely repulsive and characterized by the coefficient of restitution which specifies which fraction of the kinetic energy is dissipated. Wet granular particles are covered by a thin liquid film. When two particles come into contact, the films merge and a capillary bridge is formed, exerting an attractive force on the particles. As the particles separate from each other again, the bridge stays intact up to a critical distance . At this point the bridge ruptures Willet et al. (2000) and a fixed amount of energy is dissipated. Thus wet granular particles are characterized by a hysteretic attractive interaction and a well defined energy which is dissipated when a capillary bridge ruptures.

The existence of a well defined energy scale (and corresponding time scale), which is absent in dry materials, is the essential microscopic ingredient not only of wet granulates but also of cohesive gases. In fact the liquid bridge can be thought of as a particular realisation of a more general cohesive force. A particularly important aspect of free cooling in cohesive gases is the aggregation process which sets in, when the kinetic energy falls below the bond breaking energy. Wet granular systems may provide a realisation of various aggregation models and so-called sticky gases Carnevale et al. (1990), where particles move diffusively or ballistically until they collide and get stuck to an aggregate which is thereby growing. Such models have attracted a lot of interestLiang and Kadanoff (1985); Jiang and Leyvraz (1993, 1994); Carnevale et al. (1990); Alves and Jr. (2006); van Dongen and Ernst (1985); Westbrook et al. (2004); Jullien and Kolb (1984); Trizac and Krapivsky (2003); Trizac and Hansen (1995), due to a wide range of applications ranging from the formation of dust filaments, snowflakes and clouds to the size distribution and impact probability of planetasimals in accretion discs.

Kinetic properties of granular gases have been discussed mainly for dry materials. In particular, free cooling has been studied extensively Ben-Naim et al. (1999); Nie et al. (2002), and it was shown that the dissipative interactions are responsible for many novel phenomena, unexpected from the kinetic theory of molecular gases: The particles’ velocities are not distributed according to a Maxwell- Boltzmann distribution Goldshtein and Shapiro (1995), equipartition does not hold Huthmann and Zippelius (1997); Garzo and Dufty (1999); W. Losert and Gollub (1999), a spatially homogeneous state is generically unstable Goldhirsch and Zanetti (1993), and linear and angular motion are correlated Brilliantov et al. (2007).

Much less is known about wet granular media, which have been addressed only recently Thornton et al. (1995); Lian et al. (1998); Huang et al. (2005); Herminghaus (2005); Zaburdaev et al. (2006); Fingerle and Herminghaus (2006, 2008); Fingerle et al. (2008), focussing on nonequilibrium phase transitions Fingerle et al. (2008), the equation of state Fingerle and Herminghaus (2008), agglomeration Ennis et al. (1991); Thornton et al. (1995); Lian et al. (1998), shear flow Huang et al. (2005), and cooling in one dimension Zaburdaev et al. (2006); Fingerle and Herminghaus (2006).

Structure formation in wet granulates during free cooling has hardly been studied yet and is the focus of our paper which is organized as follows. In Sec. II we introduce the model and discuss the decay of the average kinetic energy in Sec. III. Aggregation is discussed in Sec. IV, before we present conclusions in Sec. V. A short summary of our results has appeared in Ulrich et al. (2009).

Ii Models

In the present article, we are interested in the zero-gravity free cooling dynamics of wet granular gases. We assume the particles to be covered by a thin liquid film, as it is the case if the liquid completely wets the particle material Israelachvili (1992). The particles approach freely, until these surface films come into contact. The liquid then rapidly accumulates around the contact due to the interfacial forces. A capillary bridge forms at the contact, exerting an attractive force on the grains due to its negative Laplace pressure. This liquid bridge is stretched but stays intact (or even continues to grow) as the particles move apart. The attractive force thus remains until a certain critical separation is reached, where the liquid neck becomes unstable and ruptures. As mentioned above, the hysteretic formation and rupture of the bridge gives rise to a characteristic loss of energy, , which depends upon the thickness of the liquid film wetting the grains.

In order to design a suitable model, a few words on the details of this process are in order. The formation of capillary bridges is quite fast in real systems. Between typical grains of one millimeter diameter it takes less than a millisecond. It is clear, however, that this formation cannot in general be considered instantaneous if the velocity of the impacting grains, , is large. If the time scale of the impact process, which may be written as , is of the same order or even smaller than the time of capillary bridge formation, the accumulated liquid volume of the bridge, and hence , will be smaller than for slow impacts. However, this will not greatly affect the main features of the wet system, in particular as to its characteristic difference from the dry granulate. In order to see that, we compare the effective restitution coefficient of the dry and of the wet system. This is shown in Fig. 1, where the restitution coefficient for the dry system is shown as the dotted curve. It tends to be mildly depending on impact energy Brilliantov and Pöschel (2004), , with a negative slope throughout. The effective restitution coefficient of the wet system, , is shown as the solid curve, assuming constant . In strong contrast to the dry system, it has a zero at , and a markedly positive slope. This illustrates the dramatic difference between these two systems. The dashed line qualitatively accounts for the effect of finite formation time of the capillary bridge. Since must stay below one, the difference between the solid and the dashed curve is very limited, and the qualitative picture concerning the comparison of dry and wet granular gases remains unchanged.

Figure 1: Restitution coefficients for dry (dotted) and wet (solid and dashed) granular systems, plotted vs the impact energy in units of the wet energy loss, . The main feature in the wet case is the zero at , which is unchanged if the finite formation time for capillary bridges is taken into account (dashed curve).

Our system consists of identical and spherical particles with diameter and mass in a three-dimensional cubic volume . The particles have a hard core interaction, such that two particles are reflected elastically, if their centers of mass reach the hard-core distance, which is the particle diameter .

To account for the liquid film, a liquid bridge is allowed to form between a pair of particles if they come close enough (“close enough” is specified later). When these particles are moving apart and their distance exceeds the bond breaking distance , the liquid bridge will break and a fixed amount of kinetic energy is dissipated; thereby, momentum is conserved and the relative velocity changes to according to


with the reduced mass . If, however, the relative kinetic energy is smaller than , the particles are elastically reflected towards each other. The effect of the capillary force, which is present in reality for distances up to , is thus solely modelled by the enrgy loss which occurs when . This has been shown before to be a very good approximation Fingerle et al. (2008), and enables event-driven simulations as discussed below. For the formation of the liquid bridge, we distinguish between two models:

In the thin film model, the liquid bridge forms when the particles touch, i.e. the distance of their centers is equal to . This model assumes that the liquid film covering the particles is infinitesimally thin and the capillary bridges form a thin liquid neck, which breaks off at the critical distance .

In the thick film model, a liquid bridge forms as soon as particles come closer than the critical bond breaking distance . This model assumes that the outer diameter of the liquid film is and its shape stays spherical and is not deformed by the particles. Although this may seem unphysical, we include this case in our study because similar assumptions have been used in many simulation studies in earlier articles. As it will turn out, the differences in most of the results are only minute. The two models are illustrated in Fig. 2.

Figure 2: Illustration of the thin film model and thick film model. In the thick film model, the liquid bridge forms, as soon as the bond breaking distances overlap. The same initial configuration in the thin film model does not create a liquid bridge, since the hard cores of the particles do not touch. Thus, the particles just pass by.

In general there is some energy being transfered to the atomic degrees of freedom of wet grains as well. In this paper we are going to neglect this dissipation mechanism because it is usually small as compared to the energy loss due to the breaking of capillary bridges, especially if the granular temperature is small. However, we want to point out that such a dissipation mechanism can easily be incorporated in the simulations, replacing the elastic reflection by incomplete normal restitution. We restrict ourselves here to perfectly smooth particles, such that translational and rotational motion are decoupled. Furthermore, we investigate free cooling only, so no gravity is present, and no energy is injected into the system.

The particular way of accounting for the liquid film used in these models makes it possible to use an event-driven simulation scheme. The possible events are the reflection of the particles at the hard core distance and the crossing of the bond-breaking distance . As mentioned above, we have previously compared event-driven simulations of the wet system with full molecular dynamics simulations integrating the equation of motion Fingerle et al. (2008). We found good quantitative agreement in the results of both methods, justifying the event-driven approach we chose exclusively for the present study.

We use dimensionless units such that , particle mass and particle diameter . The bond-breaking distance is chosen as , unless noted otherwise, and volume fraction, , is varied from up to . We use periodic boundary conditions in the - and -direction and hard walls in -direction.

Iii Cooling Dynamics

We define the granular temperature and investigate its decay in time from a given initial value . In all our simulations we choose . Simple arguments can be used to derive an analytical form of the temperature decay. In each collsion a capillary bridge ruptures with probability , giving rise to dissipation of a fixed amount of energy, the bond-breaking energy . Particles collide with frequency , so that the average loss of energy per unit time is given by:


The factor takes into account that two particles are involved in one bond-rupture.

iii.1 Early stage of cooling

In the early stage of cooling the average kinetic energhy per particle is much larger than the bond breaking energy, so that and almost every collision gives rise to dissipation by . For a dilute gas, the collision frequency


is well established, with the particle density and the pair correlation function at contact (e.g. Brilliantov and Pöschel (2004)). The two models differ only in the cross section (see Fig. 2), which is given by in the thin film model and in the thick film model.

The only temperature dependent quantity remaining on the right hand side of Eq. (2) is the collision frequency, from (3), giving rise to the following simple equation:


which is solved by . Insertig the prefactors and the initial value , one obtains, similar to Haff’s law Haff (1983), an analytical form of the decay of the temperature:


with a charecteristic time scale


Note that, in this simplified model, the assumption that every collision causes an energy loss gives rise to a time-scale after which all energy is dissipated. Even though this assumption does not hold for all times in the simulation (since the bonds do not break anymore if the relative kinetic energy is too small), the timesacle has a clear physical relevance. It sets the time after which the temperature is comparable to the bond-breaking energy and after which persistent clusters will form. In Fig. 3 the evolution of the granular temperature from the simulation is compared to (5) for different volume fractions: .

Figure 3: (color online) Decay of the granular temperature for the thick film model and volume fractions (from left to right) . particles are fixed. The corresponding solid lines show the analytic form (5) with a decay to zero at time , given in (6). At that time the temperature of the simulated granulate shows a rapid transition to a value below the bond-breaking energy . In the inset temperature data are plotted versus scaled time , such that data for different volume fractions collapse onto a single curve.

In the simplified cooling law (5), the volume fraction only enters into . Hence we try to superimpose the data by scaling time with . As can be seen in the inset of Fig. 3, the data obey the expected scaling well, except for the long time limit, which has different asymptotic behavior and is treated in the next section.

In Fig. 4, we compare data from the thin and thick film model for two volume fractions, and . The difference is solely due to different scattering cross-sections, entering in (6) and can be absorbed into the rescaling of time by .

Figure 4: (color online) Decay of the granular temperature for the thick film model () and the thin film model () for volume fractions and .

iii.2 Late stage of cooling

In the late stage of aggregation, when the system is strongly aggregated, it becomes very unlikely that a capillary bridge ruptures. Hence we observe a very slow time evolution of our system. The slow decrease of the temperature can be understood with simple arguments. The probability to break a bond is given by the probability to find a kinetic energy larger than :


We approximate the velocity distribution by a Maxwellian


and evaluate the above integral in the limit . The probability to break a bond becomes exponentially small in that limit:


The decrease of kinetic energy, as given by Eq. (2), is now dominated by the probability to break a bond. The collision frequency is not known for the clustered state, but is expected to be proportional to . Using (9) and in the rate equation (2) yields:


The prefactor is determined by the precise form of the collision frequency. Separation of variables can be used to integrate Eq.(10)


with the intital value . In the asymptotic limit and with , one finds a logarithmically slow time decay of the temperature


which is due to the very low probability to break a bond, Eq. (9). This is in strong contrast to the algebraic time decay observed for dry granular systems with coefficient of restitution . Brilliantov and Pöschel (2004)

In Fig. 5, the full solution (11) is compared to the simulation data, showing good agreement. The unknown prefactor is a fit parameter. It is noteworthy that for all densities, the temperature seems to approach a universal curve as .

Figure 5: (color online) Asymptotic time dependence for several volume fractions as in Fig. 3; data (dots) in comparison to the analytical results (lines)

iii.3 Partitioning of the energy into translational, rotational and internal degrees of freedom

After the time has passed, stable clusters emerge. For the definition of a cluster, we define particles as neighbors, if a bridge is formed and the relative kinetic energy is not sufficient to break it. This makes sure that particles which are just “passing by”, are not considered neighbors. A cluster is a set of particles connected through this neighbor-relationship. Hereby we refer to the cluster mass as the number of particles a cluster contains. Clusters defined in this way are not truly stable. Particles belonging to the cluster are occasionally kicked out, if hit by a very energetic particle.

For a more detailed understanding of the system, we investigate the cooling dynamics on the cluster level, and determine how energy is partitioned among the degrees of freedom. We split the total temperature into three constituents, the translational temperature defined via the center-of-mass velocities of the clusters, the rotational termperature defined via the angular momenta of the clusters, and the internal temperature describing the relative movement of the particles inside a cluster. These three temperatures are defined as follows.

Our definition of neighborhood relations gives rise to distinct clusters numbered by . We denote by the -th cluster with particles. Its centre of mass position and velocity are given by:


Note that single particles with are also considered as clusters.

The center of mass movement of each cluster has translational degrees of freedom, so that the total number of translational degrees of freedom of these clusters is simply . Homogeneous cluster translations are thus characterized by the translational temperature


Analogously, the rotational temperature describes the energy in homogeneous cluster rotations. The angular momentum, , of cluster is given in terms of the relative particle positions and velocities


The rotational energy of cluster with is thus given by


where the moment of inertia tensor is defined in the usual way. The case , requires special treatment, since the inertia tensor is singular. The rotational energy of a dimer can be easily calculted to , where denotes the relative velocity perpendicular to the axis of the dimer. The rotational temperature is thus


with for dimers and for larger clusters.

All the left-over kinetic energy describes the relative movement of particles inside a cluster and contributes to the internal temperature. Each cluster has a total of degrees of freedom, so that the remaining number for internal degrees of freedom is . The internal temperature is:

Figure 6: (color online) Top: Division of the total degrees of freedom into the translational, rotational and internal parts, dependent on time. Bottom: Evolution of the total (, black), translational (, green), rotational (, red), and internal (, blue) granular temperatures. Data for particles and volume fraction are shown; the behavior is qualitatively the same for all investigated system sizes. The horizontal line at corresponds to the bond breaking energy.

Fig. 6 (top) shows how the total of degrees of freedom divide up into translational, rotational, and internal degrees of freedom. The corresponding temperatures are shown in the lower half of the figure. As one might expect, for almost all degrees of freedom are translational, since most clusters are just single particles, and . Keeping in mind that two particles are only defined as neighbors if their relative velocity is not sufficient to break the bond, only stable clusters (mostly dimers) enter the internal and rotational temperatures, and therefore for . 111The factor is due to the relation between temperature and energy.

In the transitional regime , when the number of intermediate size clusters increases, the rotational degrees of freedom become important. Larger objects can have higher rotational energies without rupture 222roughly speaking, the maximum rotational energy of a cluster with radius and mass is , where the maximum rotational frequency is limited by the centrifugal force . This yields . In our case the bond breaking energy is related to the maximum force on the particles by , with the freely movable distance of a particle ., therefore the growing clusters obtain energy from caught particles, and thus increases until reaching the value of . After that, the energy of the incoming lumps is not sufficient to increase any further.

In contrast to the homogeneous cluster rotations, the internal degrees of freedom which have higher energies than will in most cases result in a bond rupture, independent of the cluster size. Therefore, decreases monotonically. At late times , large clusters have formed, thus almost all degrees of freedom are internal and .

Iv Aggregation

When the average kinetic energy per particle is comparable to the bond breaking energy, , the system starts to form aggregates, which seem to grow in a self-similar process. In the following we are going to analyze these aggregates and compare them to cluster-cluster aggregation Jullien and Botet (1987) models. As time proceeds, larger and larger clusters are formed. We observe a spanning or percolating cluster for all finite densities, and ultimately all particles and clusters have merged into a single cluster.

Figure 7: (color online) Snapshot of the system with volume fraction and particles taken at time ; the largest cluster (grey) contains 22% of the particles. Particles of the same cluster have the same color shade.
Figure 8: (color online) Same as Fig. 7 for ; the largest cluster contains 99% of the particles.

Figs. 7 and 8 show snapshots of a system at and with small volume fraction, . At the smaller time the system is not yet percolating, even though rather large clusters have already formed, the largest one (in grey) contains of all particles. The second snapshot, taken at a much longer time, shows a spanning cluster. At such large times the average kinetic energy is much smaller than the bond breaking energy (), so that bonds almost never break up. The cluster shown is already well beyond the critical time for percolation with of the particles in the cluster.

Fig. 9 shows the evolution of the cluster mass distribution , which is the number of clusters containing particles at time . One can clearly see that after some time, , which depends on volume fraction, the largest cluster emerges from the rest of the distribution. For all volume fractions a gelation transition was observed at the percolation time . The critical behavior of the gelation transition is still controversial. Since aggregation is a nonequilibrium process, there is a priori no reason that it should be in the same universality class as the corresponding equilibrium percolation transition. Yet there is some evidence in favour of this conjecture. Gimel et al. Gimel et al. (1995) observe a crossover from self-similar growth at small times and volume fractions – called the flocculation regime – to the percolation regime around . In the latter they observe critical exponents as in standard percolation theory. Kolb and Herrmann Kolb and Herrmann (1985) on the other hand obtain values for the fractal dimension of the percolating cluster, distinct from percolation theory as well as from flocculation theory. Both studies refer to diffusion limited cluster-cluster aggregation.

Figure 9: Histogram of the cluster mass distribution dependent on time, for volume fraction and . The number of clusters at the respective time and size is color coded on a logarithmic scale so that the single largest cluster is visible. At one can see the large cluster emerging, clearly distinguishable from the rest of the distribution.

In this paper we do not analyze the gelation transition in detail but defer such a discussion to future work. Instead we investigate two regimes in detail in the following:

a) The self-similar growth process, or flocculation regime, which is present for small times and volume fractions.

b) The properties of the final cluster which emerges, when (almost) all particles have aggregated to form one large cluster.

iv.1 Self-similar growth

iv.1.1 Fractal dimension of the aggregates

A central quantity of aggregation models is the fractal dimension of the aggregates. It is usually determined from the radius of gyration as a function of cluster mass. We consider a cluster of particles with positions and define its radius of gyration by (see e.g. Stanley and Ostrowsky (1986))


If the clusters are fractal we expect a scaling relation for large of the form


which yields the fractal dimension . This method is commonly used in aggregation models, where particles move diffusively, ballistically, or are interacting and stick to the aggregate once they touch it Alves and Jr. (2006); Westbrook et al. (2004); Jullien and Kolb (1984); Meakin (1991).

In Fig. 10 we show the radius of gyration for a sytem of particles at volume fraction . Several snapshots of the ensemble of growing clusters have been taken at times with the percolation time , when a spanning cluster is first observed. The data scale well according to Eq.(20), some scatter is observed for the largest masses, corresponding to times close to the percolation transition.

Figure 10: (color online) Radius of gyration as a function of cluster size for a system of 262144 particles at volume fraction ; different colors/shades correspond to simualtion times between (yellow/light gray) and (black)); The slope of the solid line corresponds to ; inset: fractal dimension as a function of time, extracted from the slope of the curves in the main figure.

In contrast to aggregation models, where the clusters are static and do not break up, we occasionally do observe the breaking of bonds. In addition there are internal deformations of the clusters during growth, so that the fractal dimension could depend on time. We have therefore checked the relation between and for many instances of time and show the fractal dimension as a function of time in the inset of Fig. 10. As can be seen from the Figure, there is no systematic dependence on time, and the fractal dimension is close to .

iv.1.2 Cluster size distribution

All information about the connectivity of the clusters is contained in the cluster size distribution , the number of clusters of size at time . In Fig. 11 we show for a system with and . The time interval has been chosen such that (for this volume fraction). In this time interval the mean cluster mass increases roughly by a factor of 30.

Figure 11: (color online) The cluster mass distribution . The different graphs represent different times, which are increasing from top to bottom (left side of the graph). The inset shows how the mean cluster mass increases during the investigated time period.

It has been suggested (e.g. Meakin (1991)) that for aggregating systems the mass distribution evolves towards a self-preserving scaling form, independent of the initial distribution:


where the time dependence is only contained in the mean cluster mass


This scaling form has been applied sucessfully to various aggregating systems Vicsek and Family (1984); Botet and Jullien (1984); Jiang and Leyvraz (1993, 1994); van Dongen and Ernst (1985); Meakin (1991); Trizac and Hansen (1995), involving fractal as well as non-fractal objects. Mass conservation requires Meakin (1991).

We plot in Fig. 12 the scaling function for the same data sets as in Fig. 11. We expect scaling to hold only in the aggregation regime, i.e. for times not too close to , where the system gels (see sec. IV.1.3). Hence we restrict ourselves in Fig. 12 to times . We have also left out the data points for , i.e. clusters consisting of single particles. As can be seen from Fig. 12 the data scale very well. Deviations occur only for times close to the percolation transition (not shown here), where they should be expected.

Figure 12: (color online) Rescaled cluster size distribution from eq. (21) versus the normalized cluster mass . The color coding as in Fig. 11 is used.

iv.1.3 Number of clusters

Another characteristic of a realisation of clusters is simply the total number of clusters , which decreases as aggregation proceeds. As long as the system is in the scaling regime (i.e. relation (21) is fullfilled), the mean cluster mass, and the number of clusters are simply related: . However, as mentioned above, the scaling relation (21) only holds in the aggregation regime and is expected to break down as the percolation transition is approached. At that point, should diverge due to the formation of a spanning cluster. On the other hand, there is still a large number of smaller clusters coexisting with the macroscopic cluster, so that remains finite at the percolation transition.

The aggregation of particles to larger objects has been investigated for various ballistic aggregation models Carnevale et al. (1990); Jiang and Leyvraz (1994); Trizac and Krapivsky (2003); Family and Vicsek (1985), where spherical particles of mass and diameter move ballistically, until two of them collide to form clusters irreversibly. In a particularly simple model, one assumes that two colliding particles form one larger spherical particle with conserved momentum and a mass equal to the sum of the two particles masses, so that is always equal to the number of initial particles contained in a given cluster. For spatial dimension , the diameter increases like , assuming the particles to be compact spheres which conserve volume when merging. For this model, a mean field theory Jiang and Leyvraz (1994) and simple scaling arguments Carnevale et al. (1990); Trizac and Krapivsky (2003) yield the dependence of the expected average mass on time like with an exponent (assuming ).

Since the aggregating clusters in our system are not compact, but fractal objects with fractal dimension , the assumption for the diameter does not hold and must be changed to . With this assumption, we follow the scaling arguments of Trizac et al. Trizac and Krapivsky (2003), and find the scaling relation between and .

We assume that the number of clusters per volume, , is reduced by one whenever two clusters collide:


The collision frequency Brilliantov and Pöschel (2004) is approximately given by with the linear dimension of the cluster and its typical velocity. The average momentum should scale as Trizac and Krapivsky (2003), and therefore


Plugging in all these scaling relations as well as , one obtains:


which is solved by


where the integration constant is the onset of cluster growth. In our context 333As one can see in the main plot of Fig. 13, the actual onset of cluster growth is not exactly at , but a little bit earlier.. This implies the following growth law for the mean cluster mass in the scaling regime:


which generalises the result for compact objects, with to fractal ones with .

In Fig. 13 we show how the number of clusters decreases over time as larger and larger aggregates form for .

Figure 13: (color online) Evolution of the mean cluster mass. Labeling and parameters as in Fig. 3. For the inset, the origin of the time-axis has been shifted to the transition point to investigate the scaling relation . The solid line has a slope of .

The inset of fig. 13 investigates the scaling behavior (26), with the origin of the time axis shifted to the transition point . One can see that the slope of , obtained from (27) for and is in good agreement with the simulation.

iv.2 Properties of the asymptotic cluster

The fractal dimension of the largest cluster – well beyond the percolation transition for most volume fractions – will be the main focus of this section. In particular we determine its fractal dimensions and coordination numbers.

iv.2.1 Fractal dimension from radius of gyration

One way to determine the fractal dimension is the radius of gyration, as was done in Sec. IV.1.1 for aggregates. Here, however, we only have one large cluster and have to find a way to obtain the function as a function of cluster size . We implement this in following way: Starting from a random particle of the cluster, we mark all particles that can be reached through neighbor-to-neighbor steps. Thus, for every , we get a partial cluster with particles and radius of gyration , which yields the scaling relation and the fractal dimension . For good statistics, we repeat this procedure 100 times (each with a different initial particle) and average over the obtained values of . Note furthermore that the procedure takes care that no particle is marked a second time, in order to make sure that one does not go through the cluster several times because of the periodic boundary conditions.

Figure 14: (color online) Radius of gyration dependent on the mass of the partial cluster at simulation time . The particle number is fixed and the volume fractions are (from bottom to top) . The lines along the data points are the respective fits. The outer solid lines have slopes (top) and (bottom) corresponding to fractal dimensions of and , respectively.

In Fig. 14 we show the results of this procedure for the radius of gyration as a function of for different densities. For high volume fractions we are well beyond the percolation transition and hence expect on the largest length scales of the cluster. This is clearly seen in Fig. 14, e.g. for and . On smaller length scales, however, we find a fractal dimension . For smaller volume fractions, the crossover to happens at larger masses and hence the “interior” region extends to larger scales.

iv.2.2 Fractal dimension from box counting algorithm

To further investigate the Hausdorff dimension of the largest cluster at intermediate length scales, we use the box counting algorithm Grassberger (1983); Hentschel and Procaccia (1983). The system is divided into sub-boxes of edge length . Then each box which contains or hits at least one particle is marked. In this way, we find the number of boxes necessary to cover the whole cluster. This number should scale with like


with the Hausdorff dimension .

On length scales much smaller than the particle diameter, , the system obviously behaves three-dimensionally. In this regime, the number of filled boxes is just the volume fraction times the total number of boxes , therefore:


Since our system is finite and contains a system-spanning cluster, the scaling behavior on large length scales should also be three dimensional. On this length scale, almost all the boxes should be filled, so that


In particlar, the relation must include the point , since a box of the system size includes all particles and will certainly be marked.

Only in the regime between these two limiting cases is it possible to observe the fractal dimension with the box-counting method. Comparing (29) and (30) shows that the interesting range is proportional to , which only depends on the volume fraction, but not on the particular choice of the system size. A schematic plot is given in Fig. 15 where the number of boxes containing particles is plotted against the edge length of a box.

Figure 15: Schematic double logarithmic plot of the box size versus the number of boxes of that size needed to cover the cluster. The negative slope is the fractal dimension. We expect three scaling regions: For small and large , the system should behave three dimensionally, and the region in between yields the non-trivial fractal dimension. If only the particle centers are considered, the algorithm simply counts the number of particles in the cluster for resulting in a horizontal line (dotted line).

For numerical reasons, it is very tedious to observe the expected slope of for small , because of the vast amount of boxes to account for. Since this regime is not relevant anyway, it has only been investigated exemplarily and is reached for . For all other runs we simplify the algorithm and only use the centers of the particles, i.e. a box is only marked, if a particle center is inside. With this definition the number of boxes needed to cover the system for small box sizes is just the particle number , resulting in a horizontal line on the left side of the graph, instead of the slope (dotted line in Fig. 15).

Figure 16: (color online) Top: versus at time and volume fraction for the box counting algorithm; particle number is varied from bottom to top, according to . The straight lines are fits to the data to the left and right of the cross-over point , which is also a fitting parameter and shown as a star (). Bottom: As top, but normalized by cluster mass; the solid lines have slopes (left) and (right); the vertical dashed line represents the particle size.

Fig. 16 (top) shows the outcome of the box-counting algorithm, at a time , where roughly all particles are inside the largest cluster. It yields the relation between the box size and the number of boxes of that size, needed to cover the cluster. The slope of that curve is the negative fractal dimension. The result for different system sizes, but with the same volume fraction are presented. As proposed, for all system sizes, there is a cross-over point , at which the slope changes. On length scales between and , the fractal dimension is roughly (the fits yield values between and ). Above the fractal dimension has a trivial value of about (fit values between and ), which means that on these large length scales all the boxes are filled and is therefore an indication that the cluster is system-spanning.

In the lower half of Fig. 16 we show the number of boxes normalized by the cluster mass. The data collapse well onto a single curve, obviously with the same slopes. Here one can see very well that for systems with the same volume fraction, the slopes as well as the cross-over point do not depend on the absolute system size.

Results of the box counting algorithm for different densities are presented in Fig. 17. We only include densities for which a spanning cluster has developed. The slopes of and in the two scaling regions are not affected by the volume fraction , but the size of the non-trivial region (with ) is seen to increase significantly as the density decreases. Even for the lowest density, the size of the scaling region is less than 2 decades, which makes it difficult to extract precise values for the fractal dimension. For the three most dense systems, the scaling region is less than one decade. As discussed earlier, this is an intrinsic feature of the “high” density systems, which can not be resolved by taking larger systems ( with constant ). As the inset shows, we can collapse all data on a single curve by rescaling with and with in agreement with the dependence of the crossover length on .

Figure 17: (color online) Result of the box counting algorithm, as in Fig. 16 at ; the particle number is fixed and the volume fraction varies from left to right according to ; the lines along the data points are the respective fits; the outer lines have slope (top) and (bottom).

iv.2.3 Coordination number

Given the definition of a neighborhood relation (two particles are neighbors if they have built a bridge and their kinetic energy is not sufficient to break it), we can extract the average number of neighbors of a particle, i.e. the average coordination number. In Fig. 18 we show histograms for the coordination number in the percolating cluster for two different bond breaking distances. As one would expect, these distributions are rather broad with coordination numbers between one and thirteen. The smaller bond-breaking distance (left) gives rise to a more asymmetric distribution with more weight for smaller coordination numbers.

Figure 18: Histogram of the coordination number for two different bond breaking distances (left) and (right); for both plots and .

In Fig. 19, we show the time evolution of the average coordination number for different bond-breaking distances . After a strong increase at the time , the coordination number continues to grow slowly. This slow increase is strongly suppressed in the thin film model (right) as compared to the thick film model (left). Within the thin film model the slow growth with time is further suppressed for decreasing bond-breaking distance .

Figure 19: (color online) Time evolution of the average coordination number of particles in the largest cluster ( and ); the critical break-off distances are and from top to bottom; thick film model (left) in comparison to thin film model (right).

As one can see in Figs. 18, 19, the coordination number becomes smaller for smaller . This is reasonable, because the particles can more easily collect neighbors for higher . As the average coordination number of the thin film model approaches 6, which is the isostatic value. This is demonstated in the right half of Fig. 20, where we plot the asymptotic coordination number as a function of . Here the asymptotic value is taken, when for the first time.

Figure 20: Influence of the bond-breaking distance on the final value of the average coordination number for the thin film model (right, and ) and on the fractal dimension of the thick film model (left, system and colors as in Fig. 19).

Naively one might expect that the increase of the coordination number with larger is caused by a compactification and therefore accompanied by an increase of the fractal dimension. However, as can be seen in the left part of Fig. 20, there is no significant influenece of on the development of the fractal dimension. Thus, we conclude that this compactification is mostly occuring on the single particle length scale and therefore increasing the average coordination number, but not influencing the stucture on larger length scales 444Note that there is also a very slight increase of the fractal dimension and therefore also a very slow compactification on larger length scales. However these effects are much less pronounced than the change of the coordination number..

V Conclusions

We have analysed a simple model of a wet granulate allowing for large scale event driven simulations. A central feature of wet granulates is the existence of an energy scale associated with the rupture of a capillary bridge between two grains. This energy scale has important consequences not only for the phase diagram Fingerle et al. (2008) but also for the free cooling dynamics investigated in this paper. The most important feature is a rather well defined transition at a time , when the kinetic energy of the particles becomes equal to .

For the particles are energetic enough to supply the bond breaking energy , so that very few collisions result in bound pairs and most particles are unbound. Cooling is very effective in this regime, but drastically different from a dry granulate. Whereas in dry granulates the dissipated energy is proportional to the energy of the colliding particles, in wet granulates the dissipated energy is , independent of the energy of the colliding particles so that . Consequently Haff’s law does not hold and is replaced by for . The simulations are in very good agreement with this cooling law for .

For , the kinetic energy of the particles is too small to provide the bond breaking energy, so that larger and larger clusters form. We call this regime the aggregation regime and analyse the properties of the aggregates. For not too long times and sufficiently small volume fractions, we observe flocculation characterized by nonoverlapping, weakly interacting clusters. The fractal dimension of the aggregates is approximately . The cluster size distribution follows a simple scaling form, , which has been applied successfully to different aggregation models before. The increase of the typical cluster size can be undestood by a simple scaling analysis: Assuming that clusters irreversibly stick together when they hit upon each other and that their radius grows with the number of particles like , yields a cluster growth . This scaling relation shows good agreement with the simulation for fractal dimension .

At larger times, a spanning cluster forms, and a gelation transition is observed for all finite volume fractions. At the gelation transition a spanning cluster coexists with many small ones, wheras at very long times almost all particles are connected to one large cluster. On the largest length scales the final cluster is no longer a fractal but compact, as one would expect for a spanning cluster in the percolating phase. On smaller length scales, however, we find fractal structures with . The range where a nontrivial fractal dimension can be observed increases with decreasing density as .

Even on the longest time scales, the temperature continues to decay. In this regime the limiting process is the breaking of a bond. The probability for this process becomes exponentially small as the temperature goes to zero. Hence the cooling law for high temperatures is replaced by in very good agreement with the data.

Several extensions of our work might be interesting. So far we have completely neglected all inelasticities except for the bond rupture. One expects the collisions at the hard core to be dissipative as they are in dry granular media. In the simplest model these could be described by normal restitution. Furthermore real wet grains experience frictional forces, coupling translational and rotational motion of the grainsBrilliantov et al. (2007). We are not aware of any such studies for wet granulates.

We gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through Grant SFB 602/B6.


  • Jaeger et al. (1996) H. M. Jaeger, S. R. Nagel, and R. B. Behringer, Rev. Mod. Phys. 68, 1259 (1996).
  • de Gennes (1999) P. G. de Gennes, Rev. Mod. Phys. 71, 374 (1999).
  • Duran (2000) J. Duran, Sands, powders and grains: An introduction to the physics of granular media (Springer, New York, 2000).
  • Brilliantov and Pöschel (2004) N. V. Brilliantov and T. Pöschel, Kinetic Theory of Granular Gases, vol. 1 (Oxford University Press, Oxford, 2004).
  • Kudrolli (2004) A. Kudrolli, Rep. Prog. Phys. 67, 209 (2004).
  • Umbanhowar et al. (1996) P. B. Umbanhowar, F. Melo, and H. L. Swinney, Nature 382, 793 (1996).
  • (7) J. Blum et al., Phys. Rev. Lett. 85, 2426 (2000); F. C. Bridges, A. Hatzes and D. N. C. Liu, Nature Lett. 309, 3333 (1984).
  • Herminghaus (2005) S. Herminghaus, Advances in Physics 54, 221 (2005).
  • Willet et al. (2000) C. D. Willet, M. J. Adams, S. A. Johnson, and J. P. K. Seville, Langmuir 16, 9396 (2000).
  • Carnevale et al. (1990) G. F. Carnevale, Y. Pomeau, and W. R. Young, Phys. Rev. Lett. 64, 2913 (1990).
  • Liang and Kadanoff (1985) S. Liang and L. P. Kadanoff, Phys. Rev. A 31, 2628 (1985).
  • Jiang and Leyvraz (1993) Y. Jiang and F. Leyvraz, J. Phys. A 26, L179 (1993).
  • Jiang and Leyvraz (1994) Y. Jiang and F. Leyvraz, Phys. Rev. E 64, 2148 (1994).
  • Alves and Jr. (2006) S. G. Alves and S. C. F. Jr., Phys. Rev. E 73, 051401 (2006).
  • van Dongen and Ernst (1985) P. G. J. van Dongen and M. H. Ernst, Phys. Rev. Lett. 54, 1396 (1985).
  • Westbrook et al. (2004) C. D. Westbrook, R. C. Ball, P. R. Field, and A. J. Heymsfield, Phys. Rev. E 70, 021403 (2004).
  • Jullien and Kolb (1984) R. Jullien and M. Kolb, J. Phys. A 17, L639 (1984).
  • Trizac and Krapivsky (2003) E. Trizac and P. L. Krapivsky, Phys. Rev. Lett. 91, 218302 (2003).
  • Trizac and Hansen (1995) E. Trizac and J.-P. Hansen, Phys. Rev. Lett. 74, 4114 (1995).
  • Ben-Naim et al. (1999) E. Ben-Naim, S. Y. Chen, G. D. Doolen, and S. Redner, Phys. Rev. Lett. 83, 4069 (1999).
  • Nie et al. (2002) X. Nie, E. Ben-Naim, and S. Y. Chen, Phys. Rev. Lett. 89, 204301 (2002).
  • Goldshtein and Shapiro (1995) A. Goldshtein and M. Shapiro, J. Fluid Mech. 282, 75 (1995).
  • Huthmann and Zippelius (1997) M. Huthmann and A. Zippelius, Phys. Rev. E 56, R6275 (1997).
  • Garzo and Dufty (1999) V. Garzo and J. Dufty, Phys. Rev. E 60, 5706 (1999).
  • W. Losert and Gollub (1999) J. D. A. K. W. Losert, D. G. W. Cooper and J. P. Gollub, Chaos 9, 682 (1999).
  • Goldhirsch and Zanetti (1993) I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 (1993).
  • Brilliantov et al. (2007) N. Brilliantov, T. Poeschel, T. Kranz, and A. Zippelius, Phys. Rev. Lett. 98, 128001 (2007).
  • Thornton et al. (1995) C. Thornton, K. K. Yin, and M. J. Adams, J. Phys. D 29, 424 (1995).
  • Lian et al. (1998) G. Lian, C. Thornton, and M. J. Adams, Chem. Eng. Sci. 53, 3381 (1998).
  • Huang et al. (2005) N. Huang, G. Ovarlez, F. Bertrand, S. Rodts, P. Coussot, and D. Bonn, Phys. Rev. Lett. 94, 028301 (2005).
  • Zaburdaev et al. (2006) V. Y. Zaburdaev, M. Brinkmann, and S. Herminghaus, Phys. Rev. Lett. 97, 018001 (2006).
  • Fingerle and Herminghaus (2006) A. Fingerle and S. Herminghaus, Phys. Rev. Lett. 97, 078001 (2006).
  • Fingerle and Herminghaus (2008) A. Fingerle and S. Herminghaus, PRE 77, 011306 (2008).
  • Fingerle et al. (2008) A. Fingerle, K. Röller, K. Huang, and S. Herminghaus, New J. Phys. 10, 053020 (2008).
  • Ennis et al. (1991) B. J. Ennis, G. Tardos, and R. Pfeiffer, Powder Technology 65, 257 (1991).
  • Ulrich et al. (2009) S. Ulrich, T. Aspelmeier, K. Roeller, A. Fingerle, S. Herminghaus, and A. Zippelius, Phys. Rev. Lett. 102, 148002 (2009).
  • Israelachvili (1992) J. Israelachvili, Intermolecular and Surface Forces (Acad. Press, San Diego, 1992).
  • Haff (1983) P. K. Haff, J. Fluid Mech. 134, 401 (1983).
  • Jullien and Botet (1987) R. Jullien and R. Botet, eds., Aggregation and Fractal Aggregates (World Scientific, 1987).
  • Gimel et al. (1995) J. C. Gimel, D. Durand, and T. Nicolai, Phys. Rev. B 51, 11348 (1995).
  • Kolb and Herrmann (1985) M. Kolb and H. Herrmann, J. Phys. A 18, L435 (1985).
  • Stanley and Ostrowsky (1986) H. E. Stanley and N. Ostrowsky, eds., On Growth and Form: Fractal and Non-Fractal Patterns in Physics, vol. 1 (Martinus Nijhoff Publishers, 1986).
  • Meakin (1991) P. Meakin, Reviews of Geophysics 29, 317 (1991).
  • Vicsek and Family (1984) T. Vicsek and F. Family, Phys. Rev. Lett. 52, 1669 (1984).
  • Botet and Jullien (1984) R. Botet and R. Jullien, J. Phys. A 17, 2517 (1984).
  • Family and Vicsek (1985) F. Family and T. Vicsek, J. Phys. A: Math. Gen. 18, L75 (1985).
  • Grassberger (1983) P. Grassberger, Physics Letters A 97, 224 (1983).
  • Hentschel and Procaccia (1983) H. G. E. Hentschel and I. Procaccia, Physica 8D, 435 (1983).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description