## Abstract

We report the results of nuclear magnetic resonance imaging experiments on vertically vibrated granular beds of mustard grains. A novel spin-echo velocity profiling technique was developed that allows granular temperature, mean velocity and packing fraction distributions within the three-dimensional cell to be measured as a function of both vertical position and vibration phase. Bimodal velocity distributions were observed at certain portions of the vibration cycle, and in general the ability to acquire time-resolved data demonstrated the significant distortions to the velocity distributions and the systematic errors in calculated temperature distributions that may arise with time-averaged measurements. The experimental behaviour was compared with predictions from a time-varying one-dimensional hydrodynamic model using the experimental parameters as input to the code. In both cases, damping of longitudinal sound waves was linked to significant volume heating effects, which contrasts with the usual heat transport mechanism (i.e. diffusion from the boundaries) currently assumed in most steady-state models. This leads to a new explanation for the counterintuitive upturn in granular temperature in vibrofluidized granular beds, based on amplification and damping of sound waves in the high-altitude region.

## 1. Introduction

Many materials in the chemical, food and pharmaceutical industries are processed, transported and stored in granular form. Over the years, an extensive engineering literature has been built up which is capable of describing the constitutive behaviour of such granular materials over certain operating ranges. However, until the past 15–20 years, the approach was largely empirical and developed at the macroscopic level. More recently, the emphasis has shifted towards a more fundamental theoretical and experimental understanding built up from the microscopic (single grain) level. For the simplest state—so-called granular gases—there is now a theoretical framework relating variables such as packing fraction, granular temperature, shear rate and heat flux for idealized monodisperse spherical particles as well as for some binary systems (Jenkins & Richman 1985*b*; Jenkins & Mancini 1989). ‘Granular temperature’ is a measure of the velocity fluctuations of the grains and is defined as the average kinetic energy of the grains by analogy to the usual thermodynamic temperature, which is the average kinetic energy of the molecules. The thermodynamic temperature is normally irrelevant in these problems because the mass of a grain is so much greater than the mass of a molecule. The dissipative nature of the collisions between grains means a granular gas needs continual input of energy to sustain its motion and is therefore far from equilibrium. Nevertheless, the dynamic behaviour of a fluidized granular material shows some striking similarities with the behaviour of a fluid in thermal equilibrium, e.g. in its microscopic structure (Warr & Hansen 1996), self-diffusion properties (Wildman *et al.* 1999) and convection behaviour (Wildman *et al.* 2001*a*). In spite of the progress, however, there is still active debate about issues as fundamental as the correct form of the equation to describe heat flux in such dissipative materials (Barrat *et al.* 2005).

One of the main geometries used over the years for investigating granular gases in the absence of shear has been the vibrofluidized bed. This consists of a cell that undergoes high-frequency agitation of sufficient amplitude to fluidize the granular medium contained within it. A range of experimental techniques has been applied to such cells, including high-speed photography (Warr *et al.* 1995) and—in three dimensions—positron emission particle tracking (PEPT; Wildman *et al.* 2001*b*) and nuclear magnetic resonance (NMR; Caprihan *et al.* 1997; Yang & Candela 2000; Yang *et al.* 2002; Huan *et al.* 2004). Although all three techniques can, in principle, provide data that are resolved in time through the vibration cycle, much greater emphasis has normally been placed in such studies on the time-averaged behaviour rather than on the transient effects. On the other hand, molecular dynamics simulations and hydrodynamic modelling by Bougie *et al.* (2002) have demonstrated the probable importance of dynamic effects. One of the main reasons for the experiments described in this paper was therefore to investigate the time-varying axial velocity distributions in a three-dimensional vibrofluidized granular bed that has largely been ignored in past experimental literature.

A second motivation was to explore new mechanisms for the counterintuitive upturn in granular temperature at high altitudes observed both experimentally (Warr *et al.* 1995; Wildman *et al.* 2001*b*; Huan *et al.* 2004) and numerically (Helal *et al.* 1997; Brey *et al.* 2001; Ramirez & Soto 2003) in some vibrofluidized granular beds. The phenomenon is counterintuitive because, for a cell with smooth vertical sidewalls, the energy needed to sustain the fluidized state is injected almost entirely by the base of the cell. One would therefore expect a monotonically decreasing temperature with distance from the base, the temperature gradient in the vertical direction at any horizontal slice being sufficient for the energy flux through the slice to balance the rate of dissipation within the volume of material lying above it. One explanation for this phenomenon is that, for dissipative granular materials, the usual Fourier law may need to be modified to include a term in the packing fraction gradient, as well as the temperature gradient (Brey *et al.* 1998; Sela & Goldhirsch 1998; Soto *et al.* 1999; Barrat *et al.* 2005), i.e.(1.1)where *z* is height; *T* is granular temperature; *η* is packing fraction; *J* is energy flux; *κ* is thermal conductivity; and *μ* is an equivalent transport coefficient for the packing fraction gradient term. This term arises theoretically through the expansion of the hydrodynamic equations to Navier–Stokes order. *μ* is proportional to (1−*e*), where *e* is the coefficient of restitution, and so disappears in the case of an elastic fluid. The presence of the extra term leads to the prediction of a finite constant temperature gradient at high altitudes in a vibrofluidized bed (Brey *et al.* 2001; Ramirez & Soto 2003).

An alternative viewpoint however (J. T. Jenkins 2005, personal communication) is that, since the expansion is done to first order, terms that are quadratic in the spatial gradients, the dissipation and their products may be neglected. As a consequence, terms in the heat flux that are a product of (1−*e*) and ∇*η* should, for consistency, be ignored. From this viewpoint, the zero-flux condition at high altitudes implies a zero temperature gradient there, and therefore—for a stationary system—provides no explanation for the upturn in temperature. The correct heat flux law for granular materials is therefore still an open question. It should also be pointed out that an alternative mechanism based on convection has been proposed (Wildman *et al.* 2001*b*), though the magnitude of this effect would appear to be too small to explain the phenomenon. Dynamic effects have however been largely ignored, despite the fact that all experiments on steady-state vibrofluidized beds involve continual input of kinetic energy and are inherently dynamic. It is therefore of interest to investigate whether alternative dynamic mechanisms could be responsible for the presence of the upturn.

The structure of the paper is as follows. First, in §2, an experimental method based on NMR for providing the necessary phase- and height-resolved data within a three-dimensional system is described. This is followed in §3 by an outline of a numerical solution to the hydrodynamic equations for a dissipative granular gas, using experimentally measured parameters as input to the code. Section 4 describes the main experimental and numerical results, which are then compared and discussed in §5. Finally, in §6, we propose what we believe to be a novel dynamic mechanism for the granular temperature upturn observed in time-averaged experiments on vibrofluidized granular beds.

## 2. Experimental

### (a) Granular system

Most NMR acquisitions generate signal from a liquid phase such as water or oil within the sample. For this reason, naturally occurring grains such as poppy or mustard seeds (Yang *et al.* 2002) have been used in previous granular materials studies. For the present work, we used mustard seeds on account of the relatively large diameter and narrow size range (*d*=2.04±0.23 mm), and relatively small deviations from spherical geometry (aspect ratio=1.17±0.09). The mean grain mass was 5.66 mg.

A cylindrical cell containing the grains was machined from a permanently anti-static acetal co-polymer to reduce electrostatic charging effects, with a glass disc insert (thickness 1 mm) at the bottom. The internal radius of *R*=9 mm was constrained by the size of the NMR spectrometer bore, and the height was sufficient to prevent collisions with the top wall. The cell was mounted on top of a glass rod within the vertical bore; the rod was vibrated vertically by means of a camshaft and separate aluminium horizontal drive shaft attached to a DC electric motor. The lengths of the two driveshafts (1.0 m horizontal and 1.0 m vertical) were imposed by the need to avoid placing the DC motor in the large stray magnetic field of the superconducting magnet, thus restricting the maximum practical drive frequency to approximately 40 Hz, somewhat lower than the value of 50 Hz typically used in experimental studies using electromagnetic shakers (Warr *et al.* 1995; Wildman *et al.* 2001*b*). The cell was evacuated to reduce air drag on the grains. Coefficients of restitution were measured using high-speed photography as follows: mustard–plastic, 0.60±0.01; mustard–glass, 0.58±0.02; mustard–mustard 0.68±0.02. For the measurement of the grain–grain coefficient of restitution, each of the two seeds was suspended by a pair of thin Kevlar filaments to allow repeatable ‘head-on’ collisions with negligible post-impact rotation. The values quoted represent the mean, and standard deviation in the mean, resulting from 10 impacts (mustard–plastic), eight impacts (mustard–glass) and 20 impacts (mustard–mustard).

### (b) NMR acquisition

A novel spin-echo velocity profiling technique allows velocity distributions within the three-dimensional cell to be measured as a function of both vertical position (*z*) and vibration phase (characterized by time *t* after the trigger pulse). The pulse sequence developed for this study is described by Mantle *et al.* (in press). We chose here the vertical velocity component (*v*_{z}) since the granular temperature is normally highest in this direction, although other components could be measured equally easily at the expense of additional acquisition time. The measured signal integrates through the thickness of the cell; in effect, we traded spatial resolution in the horizontal directions (*x* or *y*) for the time and *v*_{z} resolution that are the primary focus of this paper. Previous studies on two- and three-dimensional vertically vibrated beds using high-speed photography and PEPT (Warr *et al.* 1995; Wildman *et al.* 2001*a*) have demonstrated that the main gradients in temperature and packing fraction occur along the *z*-direction, with relatively minor perturbations in the horizontal directions.

All NMR experiments were carried out on a Bruker Biospin DMX 300 spectrometer operating at a ^{1}H frequency of 300.13 MHz. Spatial resolution was achieved using a three-orthogonal axis gradient system capable of producing a maximum gradient strength of 1 T m^{−1}. The field of view in the *z*-direction was 50.0 mm and the number of data points acquired was 128, thereby giving an axial pixel resolution of *δz*=390 μm.

A 25 mm ^{1}H birdcage resonator was used to excite and detect the magnetization from the mustard seeds and the ^{1}H 90° pulse length was 32 μs. The NMR pulse sequence of Mantle *et al.* (in press) has two distinct advantages over a single-spin or stimulated-echo velocity profile sequence (used by, for example, Huan *et al.* (2004)): (i) it refocuses magnetization dephasing due to constant motion in a linear background gradient (Pope & Yao 1993; Fukushima 1999) and (ii) the readout, or spatially encoding, gradients are also compensated for velocity artefacts. Ramped gradients were used throughout to minimize extra spin dephasing from magnetic fields created by eddy currents due to gradient switching. The gradient ramp-up and ramp-down times were 100 μs. The total echo time, TE, was 2.46 ms. Velocity encoding was achieved by using 64 equal gradient increments from −0.6 to +0.6 T m^{−1}. The length of the velocity encoding gradient, *δ*, was 580 μs. The delay between velocity encoding gradient pulses, *Δ*, was 1.40 ms. The NMR signal excitation was triggered at a fixed phase of the sample vibration. Twelve increments in the vibration phase were used to sample one complete vibration cycle. The phase increments were calculated by dividing the inverse of the vibration frequency by 12. Eight scans, at a recycle time of 365 ms, were averaged to obtain a sufficient signal-to-noise ratio. Hence, the total experimental time for a single vibration frequency was approximately 38 min. Following acquisition, the raw data were zero filled to 256 data points in the velocity encode dimension, and then a two-dimensional Fourier transform was applied to each of the 12 datasets to give spatially encoded, vibration phase-dependent, velocity profiles, denoted here by *S*(*v*_{z}, *z*, *t*). These parameters enable determination of velocities within a range −0.95 to +0.95 m s^{−1} and with a resolution of 0.007 m s^{−1}. A unique feature of the data presented here is that the ^{1}H background signal resulting from the acetal co-polymer sample holder, which is seen as a vertical line in the Fourier transformed data, provides a useful reference for the phase of the vibration cycle and an internal reference for the velocity.

## 3. Numerical

The experimental system was simulated using continuum equations derived by Richman and Jenkins (Jenkins & Richman 1985*b*,*a*). The time-varying three-dimensional equations were reduced to one spatial dimension by assuming homogeneity in the horizontal directions. A modified no-particle-flux velocity boundary condition was used to allow the granular layer to smoothly fly off the bottom plate (Shattuck in press). At the bottom plate, a temperature boundary condition for a smooth boundary (Jenkins & Louge 1997) was also used, with a no-flux condition at the top. An extra dissipation term proportional to the collision rate was added to account for the collisional losses with the vertical walls of the chamber. Without this term, the mean temperature was found to be higher than that of the experiment. However, a time-averaged temperature upturn is seen regardless of the presence or absence of this term. The equations were solved using an adaptive first-order time step with a second-order spatial derivative integration scheme. Some numerical viscosity was introduced to smooth the shocks during the initial transient, but this was subsequently turned off and the simulations were run until a periodic steady state was achieved. The system was driven by sinusoidally vibrating the entire system, with the experimentally measured amplitude and frequency, as well as cell dimensions and grain parameters, as the input to the code. More details may be found in appendix A.

## 4. Results

Experiments were carried out over a range of vibration amplitudes and frequencies, and for different numbers of grains (*N*_{g}) in the cell. Most of the results presented in the paper come from a single dataset with the following parameters: amplitude *A*_{0}=1.84 mm, frequency *f*=38.2 Hz and *N*_{g}=110, which corresponded to approximately two layers in the condensed phase.

The raw data for the distribution *S*(*v*_{z}, *z*, *t*) are shown in figure 1 as a series of frames, with an interframe time of 2.18 ms, and where within each frame the horizontal and vertical axes represent *v*_{z} and *z*, respectively. Qualitatively similar behaviour was found in the numerical simulation results. Frame 1 is the first after the trigger pulse. Within each frame, the central ‘blob’ is the signal from the mustard grains. Some signal from the cell wall and base is also apparent: since the entire cell is moving upwards at the same velocity, this shows up as a vertical line whose horizontal position varies from frame to frame. The increase in signal strength near the bottom of this line occurred at the interface between the glass insert and the base of the plastic cell (marked on frame 1 of figure 1) and provided a convenient marker for the origin of the *z*-axis, defined as the mean location of the top (i.e. impacting) surface of the glass plate. The vibration amplitude *A*_{0} derived from the wall velocity was found to be consistent to within 0.6% of the value determined from the vertical position of the glass base. The amplitude of the first harmonic in the base position was approximately 11% of the fundamental, with the amplitude of each of the higher harmonics under 4%. This relatively high spectral contamination of the fundamental sine wave is likely to be due to mechanical resonances in the long drive train and represents one of the main differences between the experimental situation and the numerical model.

Integration of *S* along the *v*_{z} direction provides a signal proportional to particle density as a function of *z* and *t* which is shown in figure 2*a*. The equivalent output from the numerical code, using the experimentally determined coefficients of restitution, is shown in figure 2*b*. When plotting figure 2*b*, it was necessary to estimate the experimental phase offset. If we approximate the position of the top glass surface as(4.1)a cycle-averaged phase angle *ϕ* at frame 1 (*t*=0) was estimated as 1.38 rad from the real and imaginary parts of the Fourier transform of the measured *z*_{b} values. However, a value of 1.05 rad gave a closer match to the *z*_{b} values during the part of the cycle at which the majority of impacts occurred and was used to plot the image in figure 2*b* as well as the white lines on figure 2*a*,*b*.

Both the experimental and numerical results presented in figure 2 are qualitatively consistent and show the bed detached from the base and essentially in free fall for the first few frames, before impacting on the base from frame 8. The numerical results show somewhat less expansion of the bed than the experimental ones, suggesting that dissipation has been overestimated with this particular choice of coefficients of restitution. The main difference between the two is that the strong signal peak adjacent to the base in frames 9–12 is largely missing in the experimental dataset. This was borne out quantitatively by a drop of 25% in the *u*_{z}-integrated signal for frames 9–11 compared with the average of frames 1–5. The ‘missing’ signal corresponds to regions with intense accelerations and is most probably caused by cycle-to-cycle variations (e.g. subharmonic behaviour of the bed) or a slight drift in the vibration frequency during the course of the experiment. However, this effect is a significant one only for low *z* values on frames 9–11.

### (a) Granular temperature distributions

In order to extract quantitative data from *S*(*v*_{z}, *z*, *t*), it is necessary to separate out the signal due to the sidewall from that due to the mustard seeds, and furthermore to take account of possible bimodal velocity distributions. As stated in §3, a ^{1}H background signal is produced from the acetal co-polymer sample holder resulting in a signal from the cell sidewalls regardless of the presence or absence of any grains.

An example in which the bimodal behaviour is particularly marked is shown in figure 3 from frame 9 of a dataset with the following experimental parameters: *f*=31.8 Hz, *A*_{0}=1.68 mm, *N*_{g}=55 and *z*=4.5 mm. In general, the grains leaving the base have a mean velocity and granular temperature that differ significantly from those approaching the base, and therefore a sum of two Maxwellian distributions is the obvious choice of fitting function, *S*_{f}(*v*_{z}, *z*, *t*). We also use a narrow Gaussian function to approximate the sidewall signal, leading to the following general form of *S*_{f}:(4.2)where the index *j*=0, 1, 2 represents, respectively, the sidewall signal and the primary and secondary granular peaks, and *C*_{3} takes account of uniform background noise. Equation (4.2) was fitted to the measured *S*(*v*_{z}, *z*, *t*) on a row-by-row basis so that *C*_{j} (peak amplitude), *σ*_{j} (peak width) and (velocity offset) are in general functions of *z* and *t*. The number of free parameters was minimized first by estimating the values of *V*_{0} and *σ*_{0} from data at the top of the cell (where no grains were present) and second by including the secondary peak only when justified by the data. Only 6% of the rows containing useful signal required both peaks, and these occurred exclusively in the last four frames of the cycle; for the remaining 94%, therefore, only five parameters were free: *C*_{0}, *C*_{1}, *σ*_{1}, and *C*_{3}.

Figure 4 shows the *v*_{z} probability density function (PDF) at a fixed height of *z*=6.8 mm for each of the 12 frames. The initial starting phase was assigned by the location of the chopper slot on the drive shaft, which was essentially arbitrary; the frame ordering for this and remaining figures has been shifted cyclically for clarity of presentation, with the original frame numbers at the right of the plot to allow cross-referencing with figures 1 and 2. Figure 4 is a horizontal cross-section through the 12 frames of figure 1, after subtraction of the wall peak and the background level, and normalization of its integral to unity. The best-fit Gaussians are also plotted for comparison: considering the rapid cyclical changes in temperature and mean velocity, these are remarkably good fits to the measured distributions. The corresponding time-averaged PDF at the same height is shown in figure 5*b*, together with its best-fit Gaussian. Unlike the time-resolved plots, the time-averaged plot shows significant deviations from Maxwellian form. Furthermore, the width of the best-fit curve implies a non-dimensional granular temperature of 0.972, which is 23% higher than the average temperature calculated from the 12 frames of figure 4. This effect is more pronounced at lower altitudes (figure 5*a*, *z*=2.91 mm) and less pronounced at higher altitudes (figure 5*c*, *z*=10.7 mm).

The *z*-component of temperature quoted above was calculated for the general (bimodal) case as(4.3)The corresponding expression for mean velocity is(4.4)Equations (4.3) and (4.4) reduce to the simple forms and for the case of a unimodal distribution (*C*_{2}=0). Equivalent non-dimensional forms for the variables are defined as follows:(4.5a)(4.5b)(4.5c)and(4.5d)where *g* is the acceleration due to gravity.

### (b) Packing fraction distributions

The time-varying non-dimensional temperature and mean velocity profiles are shown in figures 6*a* and 7*a*, respectively. Integration of the fits over the velocity direction yields estimates for the packing fraction profile(4.6)which is shown in figure 8*a*. In each of these figures, the digit at the right-hand end of each plot identifies the frame number defined in figure 1, which corresponds also to the column index in the image of figure 2*a*. Corresponding results from the simulations are presented in figures 6*b*, 7*b* and 8*b*. The digit at the right-hand end of each plot in these cases identifies the equivalent simulation frame number, which also corresponds to the column index in figure 2*b*, i.e. frame 1 is defined by *t*=0 in equation (4.1) and the phase offset value *ϕ*=1.05 rad. Differences between the experiment and the model during the impact part of the cycle are sufficiently large to cause a time lag in the propagating wave of approximately 10–20% of the vibration period between the numerical and experimental results. For this reason, experimental frame 3 is displayed adjacent to simulation frame 1, and so on for the following frames, to aid the reader during the subsequent discussion section. The results of averaging the time-varying granular temperature distributions from figure 6 over a full cycle are shown in figure 9*a*.

### (c) Discussion of experimental errors

Several factors can contribute to errors in the experimental results and we give here a discussion of some of the main ones. Among the systematic errors, the ‘missing’ signal corresponding to regions with intense accelerations, already discussed in the paragraph preceding §4*a*, is probably the most significant. The other main systematic error source is the fact that the grain velocities are in effect calculated from a finite difference formula with a time-step equal to the delay between velocity-encoding gradient pulses, *Δ* (1.40 ms in these experiments). This is a good approximation when the Enskog mean free time between collisions, *τ*_{E}, is much larger than *Δ*, but underestimates the true velocity as *τ*_{E} approaches *Δ*. An order of magnitude estimate of *τ*_{E} can be made using the values *η*=0.2, and , giving *τ*_{E}=30 ms. We therefore conclude that this is not a significant error source in these experiments except possibly at the highest particle densities that may occur during the compressional phase of the wave.

The largest contribution to the random errors is non-ideal scan-to-scan reproducibility, and this can be seen as a ‘speckled’ band surrounding the signal peak on each of the frames of figure 1. The velocity distributions are fitted to each row independently of the rest, and the noise-induced errors in granular temperature and mean velocity can therefore be estimated from the spatial fluctuations in the curves of figure 6*a* and 7*a*, respectively. It is clear from these figures that the random errors are relatively small. For reasons of clarity, error bars were therefore not added to the experimental results.

## 5. Discussion

The most important general observation is the highly transient nature of both the granular temperature (figures 4 and 6) and mean velocity profiles (figures 4 and 7). A less significant time variation in the packing fraction distribution (figure 8) is also apparent.

Figure 7 shows the presence of a relatively short intense compressional pulse (region of negative gradient on the 〈*v*_{z}^{*}〉 versus *z* curve) travelling in the positive *z-*direction as a result of the collision with the base earlier in the cycle, followed by a longer, weaker dilatational phase. The speed of sound in a fluidized granular medium of spheres is in the dilute limit (Bougie *et al.* 2002)(5.1)The experimental non-dimensional sound speed may be estimated from figure 7 as 1.1, based on the fact that, between frames 1 and 8, the compressional wave travelled approximately four particle diameters. This is in good agreement with the value *c*^{*}=1.15, obtained by substituting a time-averaged temperature value of 0.8 from figure 9 into equation (5.1).

The compressional sound wave causes local heating as it passes, and hence broadens the velocity distribution. For example, frames 11 and 12 in figure 4 have the broadest velocity distributions of the 12 and are therefore the hottest phase in the cycle at this particular altitude (*z*^{*}=3.34). Figure 7 shows that this is also the point in the cycle where the compressional pulse is passing through the same location. Comparison of figures 6 and 7 shows a local heating wave (a local maximum in the granular temperature) travelling at the same speed and in phase with the compressional sound wave in figure 7. This effect was also apparent in the numerical studies of Bougie *et al.* (2002) on vibrated beds with significantly larger numbers of grains (4–15 diameters deep in the condensed phase as opposed to approximately two layers in this study). However, since the focus of that paper was on the formation and propagation of shock waves, the local temperature enhancement was not analysed in detail.

The heat equation for a three-dimensional fluidized granular medium having a non-zero velocity component only in the *z*-direction, denoted *u*, reads as follows:(5.2)where *ρ* is mass density; *p* is pressure; *λ* is the bulk viscosity; *μ*_{s} is the shear viscosity; and *γ* and *γ*_{w} are, respectively, the dissipation rates due to grain–grain and grain–wall collisions (Jenkins 1999; Bougie *et al.* 2002). Out of the six terms on the right-hand side of equation (5.2), the last two are sink terms since *γ* and *γ*_{w} are always positive and can therefore be immediately dismissed as the source of the temperature rise. The remaining four terms can, in principle, all contribute in a positive sense to the ∂*T*/∂*t* term.

The first term (−*u*∂*T*/∂*z*) represents advection of temperature due to the mean velocity of the grains. The second term (−*p*∂*u*/∂*z*) represents internal energy changes due to adiabatic compression or expansion and will be positive or negative depending on the sign of ∂*u*/∂*z*. The third term, which represents viscous dissipation, is always positive giving a non-zero average over a vibration cycle. The region of highest damping of the wave—and hence the highest heating rate of the granular gas—is the region of the graphs in figure 7 with the greatest slope. Finally, the fourth term represents redistribution of heat through diffusion.

In contrast to the usual static modelling procedures (Helal *et al.* 1997; Yang *et al.* 2002; Martin *et al.* 2005), in which heat is assumed to diffuse into the gas from the base of the cell, a very different and potentially equally important mechanism for energy transport in granular gases may therefore be volume heating within the gas, due to viscous damping of travelling acoustic waves. The relative importance of the different terms is shown in figure 10 after time-averaging over one complete vibration cycle. The terms are plotted as a function of *z*_{s}, i.e. the *z* coordinate measured in the frame of the moving shaker, because this is the frame in which the hydrodynamic equations are solved. The sum of all six terms was found to be very close to zero in this frame, confirming that a steady state has been achieved. In the high-altitude region, the third term (viscous damping) is seen to dominate over all other terms. It is matched by significant negative contributions from the fourth term, which represents diffusion of this heat back into the bulk of the vibrated bed, and from the sixth term, representing grain–wall dissipation. We consider this point further in §6.

The velocity of the centre of mass of the bed was calculated for each of the 12 frames. Over the first six frames, the variation of this velocity with time was close to linear with a gradient of −9.68 m s^{−2}, indicating that over this period the bed is essentially in free fall. Thus, for approximately half the cycle (approx. 11 ms), one has in effect a ‘microgravity’ experiment. Although of shorter duration than traditional microgravity experiments in aircraft following parabolic trajectories (Louge *et al.* 2001; Leconte *et al.* 2006), the ability to probe the interior of a three-dimensional bed may allow a ground-based facility of this type to offer a useful alternative for investigations of short-timescale phenomena in freely cooling granular gases.

The time-averaged experimental temperature distribution in figure 9 is comparable to previous time-averaged results from techniques such as PEPT and high-speed photography (Warr *et al.* 1995; Wildman *et al.* 2001*b*): a rapid decrease in temperature as one leaves the base, down to a plateau region, with a hint of a temperature upturn at higher altitudes. In the simulation results, a more marked upturn in temperature is visible over the range *z*^{*}=5–8. We consider further the issue of the temperature upturn in §6. In the remainder of this section, we discuss the observation that significant differences are apparent in figures 6–9 between the experimental and numerical results for *z*^{*} values above approximately 5–6.

As pointed out by Martin *et al*. (2005), at sufficiently high altitude the packing fraction drops to a level at which the hydrodynamic description is no longer applicable. Two possible critical packing fractions were proposed to predict the onset of the Knudsen regime. The first, *η*_{v}, occurs when the mean time between collisions is large enough for gravity to completely arrest the mean particles' vertical motion; and the second, *η*_{r}, occurs when the mean time between collisions exceeds the mean time required to cross a horizontal distance equal to the radius of the vessel. These relationships are(5.3)and(5.4)Inserting values of *T*^{*}=0.8 (taken from figure 9), and the cell geometry parameters, gives the values *η*_{v}=0.092 and *η*_{r}=0.017. The corresponding heights, estimated from the experimental time-averaged packing fraction distributions, are *z*^{*}=6.2 and 7.9, respectively. The divergence between experimental and numerical results above *z*^{*}≈6 can therefore be interpreted as being due to the non-applicability of the hydrodynamic equations in this region.

## 6. New mechanism for temperature upturn

The dramatic upturn in time-averaged temperature predicted by the model (figure 9*a*) raises interesting questions. The fact that the model is one-dimensional means that convection cannot be responsible for the upturn. Furthermore, no term proportional to concentration gradient appears in the heat flux expression, and therefore this oft-cited reason for the upturn is also not applicable in this case. One possible explanation is that the result is an artefact arising from the use of constitutive equations outside their range of applicability (i.e. a coefficient of restitution insufficiently close to 1). To test this possibility, we repeated the simulation for the close-to-elastic case *e*=0.94 (with the same base velocity and a cell height of 12*d*). The resulting time-averaged temperature distribution shown in figure 9*b* still demonstrates a significant temperature upturn. We discuss here in semi-quantitative terms the alternative physical mechanism, outlined in §5, which is based on the transport of kinetic energy through the transmission of acoustic waves.

For simplicity, only variations along the *z-*axis will be considered. We also restrict our attention to the dilute upper region of the cell, within which the time-averaged mass density is well approximated by a Boltzmann distribution for reasonably elastic collisions (Wildman *et al.* 2000)(6.1)where is a constant. The final approximation is that, within this region, the temperature variations are sufficiently small that the wave speed *c* (equation (5.1)) can be regarded as independent of *z*.

The admittance of the cell, *Y*(*z*), defined by(6.2)where *A* is the cross-sectional area, then takes the form(6.3)where *α*=*g*/*T*_{z} and *Y*_{0}=*Y*(0). The propagation of longitudinal waves in such a system is mathematically equivalent to the case of the loudspeaker horn with exponentially increasing area, for which solutions to the wave equation are well known (Lighthill 1978). The excess pressure, *p*_{e}, and fluid velocity, *u*, vary with distance *z* along the horn as follows for a sinusoidal wave of angular frequency *ω*(6.4)and(6.5)where *p*_{0}=*p*_{e}(0). Equations (6.4) and (6.5) show that sound waves can propagate into the ‘isothermal atmosphere’ region of the granular gas only if(6.6)Provided equation (6.6) is satisfied, *p*_{e} decreases as *p*_{0} exp(−*αz*/2), and *u* *increases* as *u*_{0} exp(*αz*/2), where *u*_{0}=*u*(0). In practice, the disturbance propagating into the high-altitude region of the cell comprises a wide range of frequencies, some of which are reflected since they do not satisfy equation (6.6), and the rest of which have frequency-dependent phase velocities. The situation is further complicated by the existence of significant nonlinear effects (Bougie *et al.* 2002). Nevertheless, in both experimental and numerical cases, amplification of the fluid velocity with increasing height occurs as predicted by this analysis within the region governed by equation (6.1). This is demonstrated graphically in figure 11, in which the time-varying mean velocity distributions from figure 7 (equivalent to *u*^{*}, the non-dimensional form of the fluid velocity *u*) have been differentiated numerically and then time-averaged to calculate the cycle-averaged term . We focused on the derivative rather than on *u*^{*} itself, since by equation (5.2) it is this term that dictates the rate of damping of the longitudinal sound wave. In each case, a centre difference formula was used, with a step length of *δz* and 2*δz* for the numerical and experimental data, respectively. The overestimation of damping by the simulations suggests that the experimental upturn will be less extreme than that predicted by the simulations, in agreement with the results shown in figure 9*a*.

In the low-density limit, the term is dominated by the shear viscosity, the limiting form of which is (Jenkins 1999)(6.7)and is therefore independent of *z* in the isothermal approximation. This fact, combined with the height variation of in figure 11, shows that the source term from equation (5.2) increases with height, above *z*^{*} values of approximately 5, a result that is consistent with the ‘Term 3’ curve in figure 10.

The dissipation rate of kinetic energy in the low-density limit is controlled by the grain–wall dissipation rate, *γ*_{w}, which scales as *ρT*^{3/2} (appendix (A 7)). If *T* is indeed constant, then, as *z* increases, *γ*_{w} decreases due to the exponentially decaying density *ρ*. Furthermore, a constant *T* would imply that the fourth term (diffusion of heat) is zero. However, in the steady state, the source and sink terms must balance when averaged over a complete cycle. This can be achieved, therefore, only if the average granular temperature increases with *z*, either (i) to increase *γ*_{w} or (ii) to set up a positive temperature gradient, thereby allowing a heat flux back into the bed.

In conclusion, we have shown that the term in the heat flux expression (equation (1.1)) is not necessary to explain a temperature inversion in a vibrofluidized bed, if the flow is time dependent. Equally, however, we cannot rule out the term as playing a contributory or indeed sole role to the inversion in certain cases. Previously published simulations (Soto *et al.* 1999; Brey *et al.* 2001; Ramirez & Soto 2003) have demonstrated density-driven heat currents in time-independent situations where there is no possibility of sound or shock-wave heating. Although the discussion in this section has related primarily to a phenomenon specific to vibrofluidized granular beds, a more general conclusion can be drawn, namely that the possibility of volume heating should be considered when modelling time-varying, rapid granular flows containing regions, such as free surfaces, with strong spatial variations in particle density.

## 7. Conclusions

A novel NMR technique has been developed for measuring the time-resolved one-dimensional granular temperature, fluid velocity and packing fraction distribution within a vibrated granular bed. The velocity profile sequence offers the advantage over a single-spin or stimulated-echo velocity profile sequence of refocusing magnetization dephasing due to constant motion in a linear background gradient. Using this technique, we observed strong cyclic variations in the granular temperature and fluid velocity distributions. A one-dimensional hydrodynamic model reproduced in qualitative terms many of the key features observed experimentally. A localized heating wave travelling through the bulk of the material was associated with viscous damping of a longitudinal sound wave and provided the basis for an alternative explanation of the temperature upturn sometimes observed in vibrofluidized granular beds.

## Acknowledgments

We gratefully acknowledge useful discussions with Jim Jenkins and Javier Brey during the programme on Granular Physics, Kavli Institute for Theoretical Physics, Santa Barbara (March–June 2005). This research was supported in part by the Engineering and Physical Sciences Research Council under grant numbers GR/T25040/01 and GR/R75694/01 and by the National Science Foundation under grant numbers PHY99-07949 and DMR-0134837. M.D.M., A.J.S. and L.F.G. wish to thank the EPRSC for the provision of the NMR spectrometer. J.M.H. is also grateful to the Royal Society and Wolfson Foundation for a Royal Society–Wolfson Research Merit Award.

## Footnotes

- Received October 30, 2006.
- Accepted May 23, 2007.

- © 2007 The Royal Society