Granular compaction and stretched exponentials Experiments and a numerical stochastic model

We present a stochastic model to investigate the compaction kinetics of a granular material submitted to vibration. The model is compared to experimental results obtained with glass beads and with a cohesive powder. We also propose a physical interpretation of the characteristic time τ and the exponent β of the stretched exponential function widely used to represent the granular compaction kinetics, and we show that the characteristic time is proportional to the number of grains to move. The exponent β is expressed as a logarithmic compaction rate.


Introduction
The compaction of a non-cohesive granular material (glass beads, sand, rice for example) under vibrations or vertical tapping is often well represented by a stretched exponential function [1], inspired by the Kohlrausch-Williams-Watts (KWW) relaxation model: In this expression, the packing fraction φ evolves between an initial value φ 0 and a saturated value φ ∞ , t is the number of applied taps or vibrations, τ is a characteristic time, and β an exponent.While other empirical expressions are proposed in the literature [2][3][4], we will only focus on the expression (1) throughout this paper with aim of providing a physical meaning to τ and β using a stochastic model and some experiments.
It is easy to show that τ is the time where the curve φ(ln t) exhibits an inflexion point in a log-lin space.At this particular time t = τ, the packing fraction is φ * = φ ∞ − (φ ∞ − φ 0 )/e, and the exponent β is proportional to the slope of the φ(ln t) curve: Hence β can be interpreted as a logarithmic rate of compaction.A graphical meaning of β can also be found when plotting Y as a function of ln t, In these axis, the expression ( 1) is then a straight line of slope β.We suggest in this paper a simple one-dimensional stochastic model of compaction to understand the physical origins of τ and β.Moreover, the « universality » of the KWW relaxation expression may be questioned when the granular material exhibits a cohesive property.Widely used in the industry, cohesive powders are often difficult to handle and transport.A cohesive powder usually presents a high angle of repose and a low bulk volume fraction.The cohesiveness of a powder may be measured through its ability to flow under gravity, and this ill-defined "flowability" is often described by the Hausner ratio I H , the ratio of the the tapped bulk density of the powder over the freely settled bulk density [5].A Hausner ratio greater than 1.25 is considered to be an indication of weak flowability.
It is thus interesting to trial the role of cohesion on the compaction curve, with an extension of the stochastic model, and also with some experiments on a cohesive powder under vibrations.

Stochastic model 2.1 Non cohesive model
We first describe here a simple stochastic model to simulate the compaction of a non-cohesive granular material.the model is a set of N unit grains shared out on a discretized one-dimensional space of size H 0 bounded with a static grain at the bottom z = 0 (see Figure 1a).The initial linear fraction is φ 0g = N/H 0 and we write φ ∞g = 1 the maximum packing fraction.The free spaces between two consecutive grains model the pore space between physical grains and are of the same order of magnitude as the grain volume [6].From t = 0 to t 1 , the green-labeled grain is allowed to move downwards whereas the red-labeled grain can not move.(b) Sketch of the cohesive model with 6 clusters containing 9 grains each.The blue-labeled cluster may move down of a distance of its size (from t = 0 to t 1 ) or of the available space (from t 1 to t 2 ).
At each time step, all the grains are tested in a random order.For each grain a random number r determines its ability to move: if r p g , it may move down of a space unit only if the space below is free.The grain motion probability p g is governed by the packing fraction as which is the ratio of the free volume [7] at time t by the free volume at time t = 0.Many other expressions of this probability (also named mobility) are available in the literature [7,8] are derived from statistical physics principles, but we prefer an expression which expresses a decrease of this probability from the initial state (φ 0g ) to the final state (φ g∞ ) in the simplest way.We checked that the order of tests of the N grains has no influence on the global dynamics.At the end of the loop on the N grains, the global packing fraction φ(t) = N/H(t) is simply computed with the height H(t) of the highest grain of the set at time step t.The computation stops after a predefined number of time steps.The system obviously does not evolve anymore when the packing fraction has reached its maximum limit φ ∞g = 1.To avoid random fluctuations on the results, several runs were averaged before presenting the results.
A first example of result of this model is shown on Figure 2 where the height of the packing is plotted in a log-lin space (blue curve).During the process, a compaction front is uprising.This front is located at H c (t), and the grains below H c are at the maximum packing fraction φ ∞g .If we assume that the packing fraction above H c is a constant value φ m = φ 0 , the front location must be Looking at the dashed curve in Figure 2, this assumption seems to be checked at any time.
Simulated compaction curves averaged on 20 000 runs are shown in Figure 3 for different initial packing fractions  φ 0 and for various representations.From these data, we can compute the characteristic time T and the logarithmic compaction rate b defined by These compaction curves are compared with the KWW expression, and we can notice that the numerical result (ln t, Y) is not a straight line at all time, the stretched exponential expression does not completely fit the numerical data, especially at short time.However, around t = T (indicated by the stars), the numerical data are fairly well approximated by the KWW expression, showing that τ is indeed a good approximate of T and β ≈ b.Varying both φ 0 and N, this model shows that the characteristic time is expressed by with a fitting parameter A = 0.6, and the logarithmic compaction velocity b = 0.34 The characteristic compaction time is thus governed by the number of grains N. The ratio φ 0g /φ ∞g is also a governing parameter for T and β.

Cohesive model
The previous model can be extended for a cohesive granular material.the cohesive granular system is modeled as a set of N unit grains shared out between N c clusters, each cluster containing n grains.As previously, the grains and clusters are located on a discretized one-dimensional space of size H 0 bounded with a static grain at the bottom z = 0.The initial state is prepared first by placing randomly the clusters of even size n/φ 0g without overlap, with a linear fraction of clusters φ 0c .Then the grains are randomly placed inside each cluster with a linear fraction φ 0g .The initial global packing fraction is then φ 0 = φ 0c φ 0g .
At each time step, the particles and the clusters may move according to the motion probability laws where subscript c is for clusters, g for individual grains, φ c (t) is the cluster linear fraction, and φ g (t) is the linear fraction of grains inside the clusters.
A representative result of the simulation is given in Figure 4.The compaction curve shows an obvious twostages evolution.The first stage corresponds to the fast compaction of the clusters while the second stage corresponds to the compaction of the individual grains.This approach clearly shows that two separate characteristic times T c and T g are present, where the subscripts c and g stand for clusters and grains.As for the non-cohesive model, the characteristic times are defined through the inflexion points of the φ(log t) curve, and scale on the number of clusters or individual grains, and Whatever the model parameters (φ 0c , φ 0g , N, n), the numerical results can be well fitted by an extension of the stretched exponential function expression with two exponentials: with two time-scales τ c and τ g , two exponents β c and β g and a plateau packing fraction φ p (See Figure 4a) .The fitted characteristic times τ c and τ g are again very close to the computed times T c and T g , and the exponents β c and β g may be interpreted as logarithmic compaction rates for clusters and grains respectively.

Experiments
The results from the stochastic model are compared with compaction experiments with an experimental setup based on a horizontal vibration of a vertical tank of square section 15×15 mm 2 .A first set of experiments was made with glass beads of 130 μm, a non-cohesive granular material (I H = 1.08).The compaction curves are shown in Figure 5 for 5 different initial heights.While it is difficult to compare the 3 dimensional experimental results with 1 dimensional numerical results, we were able to test one of the main result from the stochastic model.By varying the mass of grains in the shaken tank, hence varying the num- ber of grains in the system, the fitted characteristic time over the initial packing height τ/H 0 varies linearly with φ ∞g /φ 0g , as predicted by equation (7).
A second set of experiments was conducted with a cohesive UO 2 powder.This powder is made of grains of d = 30 μm diameter with rough surfaces and a Hausner ratio I H = 1.53.The figure 6 presents three sets of data for three different filling heights of the tank.Starting at a low initial volume fraction (the UO 2 grains have an intrinsic porosity), a first increase occurs at t ≈ 10 3 cycles of vibrations.A direct observation of the system during the beginning of the vibration process shows the existence of macro-cavities rising upwards at least near the front plate of the tank, and these large void structures may exist also in the bulk (Figure 7).For a larger time, a second stage of compaction occurs at t ≈ 10 4 cycles.
Despite a long experimental time (10 6 cycles), we did not observe a saturation of the volume fraction, the limit φ ∞ seems to be ill-defined in our experiment.However, the double exponential (12) fits well the experimental data, and we were able to extract the characteristic times τ c and τ g as a function of the number of particles in the system.The insert of Figure 6 seems to show a linear trend between the characteristic times and the height of the initial packing.

Figure 1 .
Figure 1.(a) Sketch of the non-cohesive model with individual grains only.From t = 0 to t 1 , the green-labeled grain is allowed to move downwards whereas the red-labeled grain can not move.(b) Sketch of the cohesive model with 6 clusters containing 9 grains each.The blue-labeled cluster may move down of a distance of its size (from t = 0 to t 1 ) or of the available space (from t 1 to t 2 ).

Figure 2 .
Figure 2. Height of the modeled packing H(t) (red curve) normalized by the initial height H 0 , height of the compacted grains H c (t) (blue curve), compared with expression 5 (black dashed line).The model parameters are N = 250, φ 0 = 0.9 from an average of 20 000 runs.

Figure 3 .
Figure 3. Simulation results (colored continuous lines) for N = 250 grains, averaged on 20 000 runs.Each curve is fitted by Eq. 1 (dashed black line) (a) Compaction curves for three different initial conditions.The star symbols indicate the inflexion points at time t = T , (b) Same data in the (ln t, Y) space.

Figure 4 .
Figure 4. Simulation result for a modeled cohesive granular material.The parameters are N = 1250, n = 50, φ 0c = φ 0g = 0.5.(a) volume fraction as a function of time in a log-lin space (blue curve) and the fitted expression (12) (black dashed line).The stars indicate the T c and T g characteristic times.(b) Same data plotted as Y = ln(ln(1/C)) as a function of ln t.

Figure 5 .
Figure 5. Experimental results of the compaction of a glass bead assembly for different initial height.The red lines are the KWW expression (1).Insert: the fitted characteristic time τ as a function of the initial volume fraction.

Figure 6 .
Figure 6.Experimental results of the compaction of a cohesive powder under vibrations (100 Hz, 7g).The red lines are the KWW2 expression (12).Insert: fitted characteristic times vs packing height.

Figure 7 .
Figure 7. Image of the UO 2 powder during the vibration process.The arrows indicate macro-cavities near the front plate.Only the top of the tank is shown here. )