Stochastic Aspects of Granular Flow in Vibrated Beds †

A stochastic approach has been taken to examine the flow behaviour of particles in a vibrated granular bed. The coordinates of a single grain were determined as a function of time using the technique of positron emission particle tracking. Following a particle for extended periods allowed a detailed representation of the displacement probability density to be built up. This showed that, ini-tially, particles dispersed according to a Gaussian law, but that the variation in the solid fraction of grains, their granular temperature and the presence of walls meant that significant distortions were observed over longer timescales. The determination of the displacement probability density allowed the spatial distribution of the granular temperature and the diffusion coefficient to be measured in both mono- and bi-disperse granular systems.


Introduction
Interest in granular materials behaviour has steadily grown over the last few decades as the importance of powder technology has been increasingly realised. The challenges to understanding the underlying fundamental mechanisms of granular flow are considerable. Whether the granular materials have been fluidised or are being stored, one finds new phenomena that are not observable in the usual states of matter. For example, in vibro-f luidised granular beds, one finds that the kinetic theory approach is broadly applicable under a range of situations [1], but that important differences arise, principally due to the dissipation and friction during particle collisions [2]. In addition, it cannot be assumed a priori that a continuum approach will be appropriate [3]. Much progress has been made in predicting granular f low behaviour using the kinetic theory approach, but until recently however, there has been little experimental validation of these studies at the single particle level. Nevertheless, in the past few years significant research has been undertaken to investigate granular dynamics in two dimensions, and more recently in three dimensions.
Much of the progress in understanding granular f lows has been through theoretical and numerical approaches (e.g. see Ref. [4]). Recently though, a number of experimental techniques have been developed to investigate two-and three-dimensional vibrofluidised granular beds. In two dimensions, a series of experiments by Warr et al using high-speed digital photography, demonstrated that highly vibro-f luidised two-dimensional granular beds can operate near to equilibrium, e.g., the granular velocities broadly follow Maxwell distributions [5,6], and the local structure of the granular bed is similar to that seen in a thermal f luid [7]. Later, a novel method of calculating granular temperature was developed [8] based upon the short-time behaviour of the mean squared displacement, which does not require detection of the collision events. This method was used to measure granular temperature concurrently with determination of the self-diffusion coefficient [1], and subsequently allowed simple kinetic theory approaches to granular beds to be validated [1].
The visualisation of the internal dynamics of threedimensional granular flows is more challenging than viewing two-dimensional arrays of grains. A number of techniques have recently been developed to probe f lows in three-dimensional geometries, including diffusive wave spectroscopy (DWS) [9], magnetic resonance imaging (MRI) [10,11] and positron emission particle tracking (PEPT) [12]. The authors have recently demonstrated the suitability of PEPT for studies of highly fluidised granular beds [13]. In that article, it was shown that variables such as the mean squared displacement and self-diffusion coefficient could be measured, and that for very dilute systems (packing fractions ȁ 0.05) the granular temperature could be measured.
In this paper, we demonstrate that stochastic approaches developed for near-equilibrium systems may be employed to develop an improved understanding of granular behaviour. We show that the displacement probability density is a useful measurand and leads naturally to the determination of variables such as packing fraction, granular temperature and the diffusion coefficient.

Stochastic Motion in Fluids
If we consider a gas of similar particles undergoing rapid motion and successive collisions, and within that gas we "tag" a small number of the particles, then if the concentration of the tagged particles is negligible, the motion of the tagged particles is equivalent to that of a single particle. In this case, the conditional displacement probability density P(r,t ) and the current j(r,t ) satisfy the continuity equation [14] P • (r,t )ѿ∇.j(r,t )҃0 (1) and the Fickian constitutive relation j(r,t )҃ҀD∇P(r,t ) where r is the position vector of a tagged particle, D is the self-diffusion coefficient and t is the time. Solution of (1) and (2) in the low frequency limit, i.e. t >> t E , where t E is the mean time between collisions, leads to the result The mean squared displacement can be calculated from this probability density function in the usual way (see e.g. Ref. [14]) and leads to the well-known Einstein relation However, this analysis is only true for unbounded, homogeneous systems such as an unconstrained thermal gas or f luid. In a granular system, the presence of gravity, as well as causing acceleration between collisions, results in a varying packing fraction profile in the vertical ( y) direction, and so D is actually a position-dependent quantity. In this situa-1 6t r 2 4Dt 1 (4pDt ) 3/2 tion, the constitutive relation (2) is modified to j(r,t )҃ҀD(r)∇P(r,t )Ҁu(r)P(r,t ) (5) and the combination of (5) with (1) results in the Smoluchowski equation P • (r,t )҃∇.{D(r)∇P(r,t )ѿu(r)P(r,t )} (6) where u(r) is the terminal velocity of a particle in a resisting medium, when an external potential field is applied.
In the case of a vibro-f luidised granular bed, u(r) ҃u( y) only, as gravitational forces act only in the vertical direction. Similarly D(r)҃D( y), since the granular temperature and packing fraction distributions on which D depends are only weak functions of x and z [15]. Knowledge of the displacement probability density P(r,t ) can lead to determination of D from inverse analysis of (6), and as we shall see, a measure of the packing fraction and the granular temperature. A single particle tracking technique such as positron emission particle tracking is ideal for the determination of P(r,t ), and developments in this area will be discussed in the following section.

Experimental Procedure
In this paper we present results using the recently upgraded Birmingham PEPT facility. Although PEPT tracks only a single particle, the automated facility allows experimental data to be logged for a considerable length of time (up to 6 hours), resulting in pseudo whole-field data for steady state systems. The technique has recently been used to investigate a number of experimental situations, e.g., rotating beds [12] and paste f low [16], and more recently it has been used to successfully analyse three-dimensional vibro-f luidised granular beds [17][18][19][20][21]. The technique uses a radionuclide that decays by positron emission, and relies on detecting the pairs of back-to-back gamma rays produced when positrons annihilate with electrons. These gamma rays are very penetrating and an accurate location can be determined from detection of a small number of back-to-back pairs. Coincident detection of two gamma rays in a pair of position-sensitive detectors defines a line passing close to the point of emission without the need for collimation. The University of Birmingham Positron Camera is a Forte dual-headed gamma camera (Adac Laboratories, Ca., USA). Each head contains a single crystal of NaI scintillator, 500҂400 mm 2 , 16 mm thick, optically coupled to an array of photo-multiplier tubes. The current maximum count rate of 4҂10 4 s Ҁ1 enables an ideal temporal resolution of 2 ms with a spatial accuracy of ȁ 1 mm. In the following experiments, a tracer particle was created by irradiating a glass ballotini sphere with a beam of 3 He, which leads to a bead that is physically indistinguishable from the remaining beads within the experimental cell. This radiolabelling process converts the available O nuclei to a radioisotope of F, which decays resulting in the emission of positrons. In a dense medium such as ballotini, a positron quickly annihilates with an electron, producing two back-to-back 511-keV gamma rays. These are detected in the pair of diametrically placed camera heads (Fig. 1). Through triangulation of successive location events, the position of the tracer particle can be located in three dimensions.
A three-dimensional granular gas was generated using a Ling Dynamic Systems (LDS) vibration system. A sinusoidal signal was fed through a field power supply [LDS FPS 1] and power amplifier [LDS PA 1000] into a wide frequency band electro-dynamic transducer [LDS V651]. This system has a frequency range of 5Ҁ5000 Hz, a maximum acceleration of 100g and maximum amplitude of 12.5 mm. A cell of dimensions 145 mm in diameter and 300 mm in height was placed on the upper surface of the vibrating piston, itself placed between the photon detectors ( Fig. 1). The cell was constructed of polymethyl methacrylate (PMMA) to limit the attenuation of the gamma rays as they travelled through the experimental apparatus, and the walls were coated with conducting copper tape to reduce electrostatic charging of the grains. The cell was vibrated at a frequency of 50 Hz. The amplitude of vibration A O was 1.74 mm. Initially, glass ballotini balls of diameter d҃5.0Ȁ0.2 mm (with an inter-grain restitution coefficient, ε, measured using high-speed photography, of 0.91 and mass m҃ 1.875҂10 Ҁ4 kg) were used as the granular medium (the number of grains, N҃700), though later 4 mm grains were also employed (m҃0.96҂10 Ҁ4 kg) in experiments on binary mixtures.
At grain speeds of aboutc ҃1 m s Ҁ1 , the PEPT camera has an accuracy in the x-y plane of about Ȁ1 mm. However, the accuracy in the z direction is substantially worse as the grain needs to be located in a direction normal to the faces of the detectors. When calculating the granular temperature (Section 6), the average behaviour in the z direction was therefore assumed to be equivalent to that in the x direction, due to cylindrical symmetry. During each experiment, the motion of a grain was followed for about one hour, resulting in up to twenty million location events. Each location event was ascribed a coordinate in space and time (x,y,z,t).

Displacement Probability Density
One expects that granular particles placed into a highly agitated granular system will act very much like a Brownian particle; the buffeting of the grain by the external particles will result in a random walk, causing diffusion of the displacement probability density according to Eq. (6). The displacement probability density (DPD) is calculated experimentally by apportioning location events to discrete height intervals. The displacements of the grain at times after these location events are then binned to create the DPD. One then follows the development of the DPD as a function of time. Figure 2 shows one example of the evolution of the density function P(x,t ) and P( y,t ) over 100 ms at steps of 10 ms, with the grain starting at 22.5 mm above the base. Initially, P is a delta function (not shown), as the particle is considered to have started at the same point every time it is located in that region, but as time progresses, the particles diffuse outward from this location. One can see this behaviour when examining the function at times t0. At short times, the density function resembles the Gaussian function expected from Eq. (3), but in the y direction, as the particles move throughout the bed the inf luence of the variations in density and the base becomes apparent and distortions appear. The final positions of particles starting from different altitudes are shown overlaid onto Figure 3. This shows clearly that the long-time probability density of the grains is Position-sensitive g-ray detectors Tracer particle Back-to-back 511 keV photons from positron-electron annihilation the same regardless of their initial position, lending further credence to the assumption of ergodicidity previously used [13]. Using such an assumption can lead to a more accurate measure of the packing fraction, h, and can be obtained by examining the residence time distribution of the grain. The assumed ergodicity means that a temporal average is equivalent to a spatial average leading to the following ex-pression: where V i is the volume of the i th segment and F i ( y) is the fraction of time spent by the particle in this volume element. The packing fraction distribution determined using this method is overlaid onto one obtained from the DPD (Fig. 4) because at extremely long times, the displacement density function is proportional to the packing fraction density. It shows the equivalence of the two methods but highlights the greater accuracy obtained by calculating the packing fraction from the residence time distribution.

Mean Squared Displacement
The mean squared displacement of the grains over time may be expressed as an integral equation in terms of P(r,t ): (r(t )Ҁr(0)) 2 P(r,t )dr (8) or in the discretised form where i is the i th interval of displacement, n r is the number of displacement bins and ∆r is the DPD bin width. A typical result for the mean squared displacement is shown in Figure 5. Here we can observe a number of important features. The evolution of the mean squared displacement is clearly partitioned into a number of regions. At the shortest times (∆t 20 ms), we can see the quadratic toe characteristic of ballistic motion [1,8,20]. At longer times, the system becomes diffusive and the mean squared displacement is linear with time. At the longest times, however, we see that the gradient of the mean squared displacement reduces. This is the effect of the caging of the particles in the cell. The walls confine the particles in the cell limiting the maximum mean squared displacement that is possible. The behaviour of the grain in the x and y directions is somewhat similar. Nevertheless the density does not vary significantly in the radial direction, and the DPD in the x direction remains broadly Gaussian at short times. As the boundary regions are reached, clipping of the density function occurs, and this is ref lected in the reduction of the gradient of the mean squared displacement. In the y direction, the presence of gravity has a profound effect, but nonetheless . This indicates that the long-time steady state behaviour of the system is broadly independent of the start location and suggests that an assumption of ergodicity is likely to be appropriate.

Fig. 4
Comparison of the packing fraction profile determined from the displacement probability density (start height ҃22.5 mm) and one produced from the time residence distribution of the particles. The two curves should in principle be equivalent, though the time residence distribution gives a more accurate picture of the solids distribution.
the saturation effect is similar; once the base is reached or the particles reach the apex of their free trajectories at the top of the cell, the DPD and the mean squared displacement are constrained.

Granular Temperature
The granular temperature E O , or mean kinetic fluctuation energy, is defined in the usual manner, (10) where m is the mass of the grain and v X , v Y and v Z are the f luctuating particle velocity components about the mean, resolved in the x, y or z directions and E X , E Y , and E Z are the granular temperatures for each of the respective directions. It is well known that due to dissipation, the granular temperature in a vibro-f luidised bed is anisotropic [1,5,8,19,20,22], and that, unlike the case of elastic spheres, equipartition of energy does not hold. In most cases E Y E X ҃ E Z [20]. This is an important consideration when probing the spatial development of a system as it causes a corresponding anisotropy in the self-diffusion coefficients [1,23]. For simplicity, however, we treat D as a scalar throughout this investigation.
We measure the granular temperature by analysing the short-time behaviour of the mean squared dis-1 3 1 3 placement curves [24]. At times of the order of or less than the mean collision time, the random walk behaviour embodied in Eq. (4) no longer holds, and in the limit t→0, one can deduce from the equation of motion for a free particle that By utilising the short-time mean squared displacement one can extract the mean squared speed, and thence the granular temperature. Bulk convection currents in the granular bed are small and are neglected in this procedure [21]. This method can be applied here for extracting the granular temperature from the particle DPDs. Figure 6 shows an example of the granular temperature profile measured in this way as a function of altitude. This is quite typical; high kinetic energy is apparent at the base of the cell that is dissipated as grains collide with those positioned above them. Boundary effects are visible close to the bottom surface as the grain-base collisions lead to deviations from the expected form of the mean squared displacement versus time curve (Eq. (11)) [20].
The relationships between D, E O and packing fraction h have been investigated in two dimensions [1]. The Enskog kinetic theory equation in two dimensions was shown to be obeyed to within 10% for packing fractions up to about 0.6 and for the restitution 110 KONA No.20 (2002)  X Y coefficient ε҃0.92. In three dimensions, the corresponding Enskog relation reads [14] D҃ ( ) . (12) where n is the number density and g O (d ), the radial distribution function at contact is given by At least two alternative expressions for D have been derived in which the Enskog theory is modified to take account of inelasticity [23,25]. However, for the case of reasonably elastic collisions (ε҃0.91 in the experiments described here), both modified expressions for D agree with Eq. (12) to within about 10%. The convective term in (6) can also be derived using kinetic theory and Brownian motion methods [14]: where g is the acceleration due to gravity.

Self-dif fusion
It can be seen from Eq. (6) that the dispersion of the grains in the system is effectively controlled by the diffusion coefficient. Provided with experimental values for the DPDs at given times, there are sufficient equations to calculate D as the system is overdetermined. This method will be presented in detail elsewhere, but will be brief ly described here.
Equation (6) can be expanded into a form that is a polynomial in the operator ∂ /∂ y. The coefficients of these terms are functions of D and u, each of which can be extracted by taking the experimental data and estimating the 1 st and 2 nd derivatives of the DPDs using a finite difference approach. There are as many DPDs as there are time steps, and as a result the random errors in using this method are reduced through the long timescale over which experiments are performed. This results in an overdetermined set of equations that then allows the coefficients, and thence D, to be extracted using non-linear regression methods. Figure 7 shows the values of the self-diffusion coefficient extracted in this way. This compares reasonably with an estimate of D which can be obtained from the gradient of the mean squared displacement in Figure 5 and predictions of D from Eq. (12), which suggests that D ȁ 0.005Ҁ0.01 m 2 /s. However, this method tends to seriously underestimate D at high altitudes and Eq. (12) shows that D is expected to scale as 1/h which it clearly does not. This is most likely an artefact of the boundaries of applicability of the Smoluchowski equation. In rarefied systems, such as those near or in the Knudsen regime, the grains tend to move as if in free f light rather than as in a random walk. Therefore, at large heights, one sees a diffusion coefficient broadly independent of packing fraction and dominated by wall collisions rather than by grain-grain interactions.
One can attempt to quantify the conditions necessary to be able to extract the self-diffusion coefficient in this way as follows: Supposing that at least three grain-grain collisions must have occurred in the period it takes for a grain to return to its original height under the influence of gravity (i.e., if less than three occur then we assume that the grains are no longer in the random walk regime), then using simple kinetic theory predictions, one can estimate for a mean speed of ȁ 0.43 m.s Ҁ1 and particle diameter of 5 mm that a packing fraction of ȁ 0.045 must exist for the extraction of D from Eq. (6) to be successful, indicating that the above method is only applicable at heights of less than about 30 mm.

Binar y systems
In reality, most granular systems do not consist of mono-sized particles; but rather contain a distribution of particle sizes. Therefore, in order to understand true powder f low behaviour, one must go beyond monodispersity and study systems containing more In the following experiments two grain sizes were used, 5 mm and 4 mm. The number of 5 mm particles N 1 was 525, and the number of 4 mm particles, N 2 was 270. The displacement probability density was measured for each phase by first using a 5 mm tracer particle, and then repeating the experiment using a 4 mm tracer. This allowed the particle coordinates for each phase in the shaker system to be determined independently. Subsequently, the methods described above for calculation of the granular temperature and the DPDs (Fig. 8) were used to determine these functions for each size phase. One can see systematic deviations in the development of the DPDs as the phases move towards their steady state distributions. The packing fraction profiles are shown in Figure 9 and highlight the differences in the phase's steady state distributions of the two phases.
The granular temperature profiles for each phase are shown in Figure 10. As previously observed [17], the system does not normally relax towards a single temperature [26,27]. The balance between energy introduced into the system by the vibrating base and that lost during inter-particle collisions and wall-particle collisions determines the granular temperature in a vibrated system. Because of differences in the mass and size of the particles, there are dif- 112 KONA No.20 (2002)  ferent solids distribution profiles for each phase, resulting in significantly different collision rates and dissipation rates. This has the effect that the particles relax towards different granular temperatures, which are coupled through the inter-phase particle collisions. The differences in the granular temperature and packing fraction propagate in turn due to differences in the diffusion coefficients. The diffusion coefficients for each phase calculated by the method described in Section 7 are shown in Figure 11. Though the diffusion coefficient of a particle through particles of a different phase is usually referred to as the mutual diffusion coefficient, in this case it is a composite of both the mutual and self-diffusion coefficients, since the particle collides with other particles of both phases. The fact that the fraction of each phase varies as a function of height means that the relative dominance of the self-diffusion and mutual diffusion components also varies. As expected, the fact that each phase has its own packing fraction and granular temperature profile results in differing diffusion coefficients for the two phases.
Knowledge of the granular temperature and the packing fraction, as well as enabling the determination of the self-diffusion coefficient, is vital for investigating granular flow behaviour at the single particle level, and provides a powerful route to validating kinetic theory methods. The measurement of the particle displacement probability density provides both a qualitative insight into the behaviour of grains in the system, and a means by which the above variables may be extracted, providing much required data for the understanding of granular flow mechanisms.

Conclusions
Particles in a granular gas are often observed to follow a random walk path. By investigating the motion of a single tracer particle in a vibrated granular bed, it is possible to estimate the displacement probability density. This function provides us with information on the system, both in a dynamic way, with the evolution of the displacement density, and on the steady state variables such as granular temperature and packing fraction and on the self-diffusion coefficient when the density is sufficiently high for the stochastic approach to be appropriate. Similar investigations on binary systems were used to quantitatively compare the diffusion coefficient for particles of different diameters.