Low-frequency oscillations in narrow vibrated granular systems.

Low-frequency oscillations in narrow vibrated granular systems.

N Rivas, S Luding and A R Thornton Multi Scale Mechanics (MSM), MESA+, CTW, University of Twente, PO Box 217, 7500 AE Enschede, The Netherlands. N.A.Rivas@utwente.nl

We present simulations and a theoretical treatment of vertically vibrated granular media. The systems considered are confined in narrow quasi-two-dimensional and quasi-one-dimensional (column) geometries, where the vertical extension of the container is much larger than both horizontal lengths. The additional geometric constraint present in the column setup frustrates the convection state that is normally observed in wider geometries. This makes it possible to study collective oscillations of the grains with a characteristic frequency that is much lower than the frequency of energy injection. The frequency and amplitude of these oscillations are studied as a function of the energy input parameters and the size of the container. We observe that, in the quasi-two-dimensional setup, low-frequency oscillations are present even in the convective regime. This suggests that they may play a significant role in the transition from a density inverted state to convection. Two models are also presented; the first one, based on Cauchy’s equations, is able to predict with high accuracy the frequency of the particles’ collective motion. This first principles model requires a single input parameter, i.e, the centre of mass of the system. The model shows that a sufficient condition for the existence of the low-frequency mode is an inverted density profile with distinct low and high density regions, a condition that may apply to other systems too. The second, simpler model just assumes an harmonic oscillator like behaviour and, using thermodynamic arguments, is also able to reproduce the observed frequencies with high accuracy.

45.70.Mg, 45.70.Vn

1 Introduction

Vibrated beds of granular materials present a wide range of different behaviours: phase separation [1, 2], Faraday-like pattern formation instabilities [3, 4], heap formation and convection [5, 6], segregation [7, 8], clustering [9] and periodic cluster expansions [10], as well as many others. These systems generally present a remarkable collection of distinct nonequilibrium inhomogeneous stable states for relatively small changes in the energy injection parameters. Hence, they are specially suited for the study of nonequilibrium phase transitions, as well as non-linear phenomena in general. Careful analysis of the microscopic mechanics behind the different transitions improves the comprehension of the complex dynamics present in driven granular systems. This gives further insight into when, and until what point, granular media behave like classical gases, fluids or solids, or whether they require an altogether different theoretical approach.

As can be seen by the aforementioned studies, the geometry of the system plays a fundamental role in determining the phenomena. Just by reducing the effective dimensionality of the system it becomes possible to observe behaviour not easily identifiable in fully three dimensional systems. The natural approach of study is then confining the grains to quasi-two-dimensional systems, where also particle-tracking methods become possible. Our study is inspired by one specific quasi-two-dimensional geometry that presents several distinct states in the energy injection parameter space: a vertical narrow box. That is, we focus on a vertically vibrated Hele-Shaw cell with the walls parallel to gravity, inside which the grains are located. The first reported classification of the different states present in this geometry was realised by Thomas et al. [11], in what would now be considered the low energy injection limit. Research then focused on the wave-like dynamics of the granular bed, and its variations with the frequency and the amplitude of oscillation [12, 13]. It was with simulations that the energy input was considerably increased, and the existence of a density inverted state was first reported [14]. This state, named Leidenfrost after the analogous water-over-vapour phenomena [15], was then experimentally studied in depth by Eshuis et al. [16, 17], as well as the buoyancy driven convection regime that is observed for higher energy inputs.

Figure 1: Snapshots of three vertical narrow box systems with the same number of filling layers , with N the total number of particles, the (dimensional) diameter of the spherical particles, and and the (dimensional) width and depth of the container, respectively; and energy injection parameters, but different widths . From left to right, , , and . The rightmost corresponds to the column geometry. Particles are coloured according to their kinetic energy.

In the following simulational study the dimensionality of the vertical, narrow box is progressively reduced until both the width and depth are only five particle diameters wide, rendering the system effectively one-dimensional (see Figure 1). We refer to this last setup as the granular column. The first direct consequence of this further confinement is the frustration of the horizontally inhomogeneous states present in the wider geometries. Particularly, the suppression of convection makes it possible to directly observe the grains collectively oscillating at a much lower frequency than the energy injection frequency. Appropriately, we call these oscillations low-frequency oscillations (LFOs). Effective frequencies and amplitudes are defined and studied in the container length and energy injection parameter space. We then argue that LFOs are an essential feature of the dynamics of the narrow vibrated geometry, but it is only in the quasi-one-dimensional column setup that they can be easily isolated from the other collective grain movement of convection. Simulational measurements confirm this, as well as a continuum description of the system, which captures the correct frequency response for high energy inputs. The frequency behaviour is actually analogous to a forced harmonic oscillator, and is obtained mainly by considering a vibrated media with a high density region suspended over a low density one. This density inverted situation is indeed present, to different extents, in both the Leidenfrost and the convective regimes.

2 Simulations

Simulations are performed using an event-driven molecular dynamics algorithm [18]. In this approach, particles move freely under the effect of gravity until an event take place, namely, a collision with another particle or a wall. The motion of the particles in between successive events does not have to be simulated: if their trajectory equation is known, time can be advanced in variable steps. This makes event-driven simulations considerably faster than usual soft-particle simulations, where time is advanced at constant steps, independent of particle interactions. However, the need of having an analytical expression for the particle motion is a strong condition that limits the possible interaction between particles. In the following, we consider the most common approach: perfect hard-spheres, which imply binary collisions and no overlap or long-range forces between particles. Collisions are then modelled by normal and tangential restitution coefficients, , and also static and dynamic friction coefficients, [19]. In order to avoid inelastic collapse we use the TC model [20], with a constant value , with the (dimensional) diameter of the spheres, and the (dimensional) gravity. (In the following, quantities without a tilde are dimensionless). That is, collisions between two entities are considered elastic if they occur more frequent than gravity timescale units. These values are known to reproduce complex behaviour observed in similar vibrated setups using stainless-steel spheres of  mm to mm in diameter [10, 21].

The setup consists of an infinitely tall container of width and depth inside which the grains can move. The boundaries of the container are considered solid, and have the same collision parameters as particle-particle collisions. The whole box (both the bottom and the side walls) is vertically vibrated in a bi-parabolic, quasi-sinusoidal way with a given frequency and amplitude . The use of a quadratic instead of a sine function gives a considerable speed advantage in simulations, as the prediction of collision times with the moving walls becomes substantially faster. Test simulations were performed using a sine function for exemplary cases, and no significant differences were observed [22]. Furthermore, we considerably varied the collision parameters and found essential phenomena to be robust, although friction was observed to play a significant role in the dynamics.

We now introduce dimensionless variables, which will be used for the rest of the text. The depth of the box is fixed, , and the horizontal width is varied in the region. is always taken so that the number of filling layers , which implies that varies in the range. Three different oscillation amplitudes are considered, . This allows us to compare with previous results, obtained for , as also to extrapolate to lower amplitudes, where the vibrating bottom wall can be considered as a spatially fixed source of energy (i.e. a temperature boundary condition).

The dimensionless gravity timescale is given by , with . Correspondingly, the dimensionless oscillation frequency is scaled as . Nevertheless, it is almost always more meaningful to measure time in periods of box oscillations, , and thus we use . In order to compare simulations with different energy injection parameters the dimensionless shaking strength is used, . Finally, the mass scale is set to unity by taking as the mass of one particle, .

Simulations are generally run for , unless otherwise stated. Particles are initially arranged in a strongly perturbed, low density crystalline state. We confirmed that this initial configuration has no influence on the steady dynamics by running a few simulations using the end state of the previous simulation as the initial configuration. Contrary to the experiments realised in [17], where the frequency of shaking is continuously increased, the energy injection parameters are kept fixed during any given simulation.

2.1 Phase Space

In order to validate our simulations, and explore further previous research, we first focus on the case, where the comparison with previous experiments and soft-particle simulations undertaken by Eshuis et al. [17] is straightforward. Event-driven simulations are able to reproduce all previously observed states, as can be seen in the phase diagram in the space presented in Figure 2a. Furthermore, a quantitative comparison is possible by looking at the transition points in the case, where the experiments were realised. There is excellent quantitative agreement, within %, for the bouncing bed-undulations and the Leidenfrost-convection transition points, but a % error in the undulations-Leidenfrost one. The deviations could in part be explained by the nature of the transitions, as they are not sharp and are seen to present wide ranges of metastability. This makes it harder to define a precise transition point value, and motivates the use of transition regions, which we show in gray.

As can be seen from Figure 2a, for , and as is increased (keeping ), the system goes through a sequence of different non-equilibrium stable states: bouncing bed, bursts, undulations, Leidenfrost, convection and gaseous (, not shown). Some of these states disappear or appear as and are varied, but their relative order remains. In the following, we will focus on the Leidenfrost and convective states, where LFOs take place. For a description of the other states, we refer the reader to [11, 17].

The density inverted Leidenfrost state owes its name to the analogous liquid-over-vapour phenomena, where a thin layer of vapour over a hot surface significantly slows the evaporation of the droplet above it, by keeping it floating over the hot surface [23]. Figure 2b shows the packing fraction and the granular temperature as a function of , for different amplitudes and frequencies of oscillation, all in the Leidenfrost state. The granular temperature is defined as twice the fluctuating kinetic energy per degree of freedom: , where is the position of particle , its velocity, and the average velocity field. Indeed, a low temperature, high density region is suspended over a low density, high temperature one. Notice that the difference in density between the solid and gaseous regions is greater for higher (red vs. black), but lower for higher (blue vs. red): these features will be relevant in our model discussion for the validity regions of a density profile approximation.

Figure 2: (a) Phase diagram of the vertically vibrated system in the dimensionless container width () and shaking strength () space, for fixed box oscillation amplitude . The equivalent box oscillation frequency is shown on the right axis. All previously reported states are seen: bouncing bed (b.b, yellow), bursts (green), undulations (purple), Leidenfrost (blue) and convection (red). Transition regions are shown in gray, and are defined by the regions of bistability of every pair of states. The borders between different numbers of convective rolls ( and ) is also delimited (dashed lines). (b) Packing fraction (solid) and granular temperature (dashed) vertical profiles, for the exemplary amplitudes and frequencies shown in the inset, and . All these systems are in the Leidenfrost state.

When is further increased, the density of the solid region is seen to progressively decrease, leading to a buoyancy driven convective state (see Figure 1). Horizontal homogeneity is lost, leading to low density regions where particles go up and circulate around high density regions, where particles agglomerate and move mainly in the horizontal directions, towards the low density regions. The number of convection rolls () diminishes with increasing , until the energy input is so high that particle motion is essentially uncorrelated and the system enters the gaseous state (, data not shown).

We now turn our attention to the lower amplitude regions. Figure 3 shows a phase diagram again in the parameter space, for different shaking amplitudes . As observed previously [17], and confirmed here for a wider range of parameters, the dimensionless shaking strength is a better parameter than the dimensionless acceleration, for the characterisation of the Leidenfrost-convection transition. On the other hand, the transition points of bouncing bed-Leidenfrost (or undulations-Leidenfrost for ) vary significantly with , but stay within when compared in . In general terms, the most significant influence of reducing is the disappearance of the bursts and undulations states; the large amplitude of the box oscillation plays a dominant role in the dynamics of these states.

We briefly remark that simulations were done until , and no new states were observed, except for the coexistence of convection and Leidenfrost states for . The possibility of this coexistence provides new insight into the nature of the Leidenfrost-convection transition, and suggests further research.

Figure 3: Phase diagram of the vertically vibrated narrow box in the shaking strength () and container width () space, for different oscillation amplitudes, as shown in the legend.

If the length is reduced further, below the limit, the frequency needed to trigger convection progressively increases, until at (a value slightly dependent on , see Figure 3) no convection was observed even for . For , undulations and bursts are also frustrated by the small size of the container. It is in this geometry that it becomes possible to observe the Leidenfrost state for higher , where low-frequency oscillations (LFOs) can be directly observed and eventually, as is increased, dominate the collective dynamics of the system.

2.2 Low-Frequency Oscillations (LFOs)

Finally, we reach the column limit, where . In order to study LFOs the evolution of the vertical centre of mass of the particles is considered, . Figure 4a shows for fixed and several different . For comparison, non-stroboscopical and stroboscopical are shown for the and cases: the distinct high and low frequencies become immediately recognisable. The amplitude of the oscillations is seen to increase from the to the order, and present an appreciable regularity in time. While at both oscillations are comparable in amplitude, and thus very hard to identify from direct observation, at they have become clearly differentiable. Although LFOs are seen to be fairly chaotic (recall that there are only 300 particles in the column geometry, hence fluctuations play a leading role), we characterise them by a constant amplitude and a single frequency , as an initial first order description.

First, let us focus on the frequency of the LFOs, , which is clearly recognisable from the power spectra of , presented in Figure 4b. The spectra are obtained by taking the discrete fast Fourier transform of over after an initial transient of , with a sampling rate of . An average is then taken over simulations with identical parameters but different initial conditions; although the shape and peaks are already recognisable from single simulations, the ensemble averaging reduces the noise considerably. The time window, the sampling rate and the transient time were varied and no significant differences were observed.

Figure 4: (a) Centre of mass evolution, , for and different dimensionless shaking strengths , as a function of time in gravity timescale units . The light colour data are taken with sub-period resolution, while dark colour data are taken every oscillation cycle at the point of maximum wall amplitude. (b) Fast Fourier transform of the centre of mass of the particles, , for and several different . The arrow indicates the direction of increasing . Different amplitudes, not shown, present the same qualitative behaviour.

All spectra present two main features: the expected delta-like peak at and its harmonics, and a broad peak one to two orders of magnitude lower, corresponding to the LFOs. The LFO frequency, , is defined as the frequency of the maximum of this broad peak. After observing the different spectra it becomes evident that depends on the energy injection parameters. Figure 5a shows for different and , remarkably scaling all LFO data. Notice that decreases as increases, i.e., the collective grain movement becomes slower as the shaking gets faster. The decay is faster than inverse linear, and can be fitted by a power with a error (not shown). Let us also notice that the length of the container makes no discernible difference, as long as the system stays in the Leidenfrost state; the decreased data in the case are due to the Leidenfrost-convection transition. The collapse of the different amplitude curves is very good for and , while for data slightly deviates. We interpret this decrease as the influence of the undulations state in the Leidenfrost regime; notice that for and the system is almost at the boundary between both states (see Figure 2a).

In order to quantify the relevance of the LFOs, we define the relative intensity of the peak, , as the normalised distance from the low frequencies asymptotic limit to the maximum of the broad peak. Figure 5b shows for different and . Although the dependency is not straightforward, it can be seen that LFOs become increasingly distinguishable from other movements until , after which there is a decline, except for the highest amplitude case. Already at oscillations should be discernible in the spectra as a peak twice as big as the low-frequency asymptotic limit. is seen to have a pronounced effect on ; higher amplitudes of oscillation lead to more pronounced LFO peaks.

Finally, we define the amplitude of the LFOs, , as the standard deviation of : . Data is considered only after , to disregard transient states. Figure 5c shows increasing in an almost linear way. The curves coincide, within their error, for and , while for all other cases is consistently smaller. Nevertheless, makes all curves comparable, further confirming its relevance for this system.

Figure 5: (a) LFO frequencies, , as a function of , for different container lengths , and shaking amplitudes , as given in the inset. (b) Intensity of , , defined as the height from the assymptotic low-frequencies value of the spectra to the broad peak, for the same data as (a). (c) Amplitude of the LFOs, defined as the standard deviation of , as a function of , for the same data as (a).

2.3 LFO’s in convective state

We now consider in detail the peculiar change of behaviour of and for in the case. This is a sign of the Leidenfrost-convection transition, still present at this container length (see Figure 2a). During convection, becomes a less relevant quantity, as there is no longer horizontal homogeneity. Nevertheless it is still possible to identify LFOs, even if the oscillations are entangled with the convective flow. The presence of LFOs in the convective regime should not be surprising if one notices that it also presents the essential feature of the Leidenfrost state: a high density, low temperature region suspended over a low density, highly agitated one, although there is an additional low density, highly convective zone above. Our model, derived in Section 3 below, suggests that when density inversion is present, LFOs exist. Figure 6 presents several different fields and snapshots that show that, indeed, density inversion is present in the convective regime, in addition to the horizontal inhomogeneity. All data is taken from the same simulation, and fields are time-averaged over after an initial transient of , with data taken every . The average velocity field, Figure 6a, clearly shows the presence of convective flow, with a small downwards band and a wider upwards region. Particles agglomerate at the bottom of the downwards flux side, as can be seen from the average density field (Figure  6b), and the two snapshots (Figures  6d and 6e). This happens when downwards and upwards particles collide, leading to a high granular temperature region (Figure  6c). Note, then, that both sides correspond to low density, high temperature regions sustaining high density, lower temperature ones, although the density and temperature profiles vary considerably from left to right. The profile is more similar to the Leidenfrost case in the upwards flow region (left in the shown figures), as in the downwards flow region the high density area presents a comparable, although lower temperature to the low density region below.

Figure 6: (a) Averaged velocity field of an system in the convective state, for and (). The colour of the arrows corresponds to the average speed, increasing from blue, green, yellow, until red. (b) Average density field of the system in (a). Colour scale from blue (low densities) to red (high densities). (c) Average granular temperature field, as defined in main text. (d, e) Two snapshots of the system taken at the minimum (d) and maximum (e) of a low-frequency oscillation. Colour corresponds to the particles kinetic energy.

2.4 Summary

Having possible experimental realisations in mind, the general picture is that LFOs are easier to observe for higher amplitude and frequencies of oscillation of the box, while keeping small; it is at these configurations that LFOs have the highest amplitudes and better defined frequencies, as quantified by and , respectively. Let us now remember that at this limit we also observed the most clear phase separation in the Leidenfrost state, with distinct low and high density regions. In our model, presented next, the separation of the phases and the confinement of the system to a one-dimensional geometry implies the existence of LFOs, and the frequency is essentially determined by the ratio of the low and high densities.

3 Continuum model

After observing the collective movement of the particles in the column geometry, an oscillator-like description naturally comes to mind. The two coexisting frequencies observed in the spectra suggest a forced oscillator model, with clearly defined forcing and response frequencies. In the following we derive such frequency behaviour from a continuum description of the granular media. We begin by considering Cauchy’s equations for mass and momentum conservation:


where corresponds to the material density, is the velocity vector, the stress tensor and the gravitational acceleration in the downwards direction, . Furthermore, the material derivative is defined as . We consider the same scaling as in simulations, with length scales in units of particle diameters , time units given by gravity , as also , taken as the mass density of a single particle, , with .

As has been observed in simulations, the dynamics of the system in the column limit is effectively one-dimensional. This immediately suggest the consideration of , and . Substituting in (1) yields


Furthermore, expanding (2), and using (3), one reaches a one-dimensional momentum conservation equation


3.1 Two phases approximation

In order to solve (4) it would be necessary to know both the density and the velocity profiles, and . Our approach consists in eliminating the -dependence from (4) by integrating in the vertical direction, and taking a first order approximation of the density profile , and average values for the vertical velocity profile . We begin by integrating (4) in the vertical direction


with the bottom boundary, , and top boundary, , dependent on time, due to the movement of the bottom wall and the free surface at the top.

The approximation of consists in dividing the system in two separate, constant density regions, inspired by the measured Leidenfrost state density profile. Let us remember that this approximation becomes increasingly better as increases and decreases, as shown in Figure 2b. Consequently, a low density region is defined where , for ; and a high density region where , for , with the position of the interface between the two regions. Figure 7 shows a schematic representation of this approximation, and the origin of its motivation. It then follows that the first integral in (5) can be expanded as


Analogously, the third integral in (5) becomes


with the height of the gaseous region, and the height of the solid region.

Notice that the second integral in (5), corresponding to the stress term, is a perfect integral, and thus only the stress boundary conditions are needed for its evaluation. We assume the stress through the system to be continuous in , and thus it is not necessary to evaluate at the interface position . Thus, from (5) we finally obtain:


3.2 Boundary conditions

It now becomes necessary to specify the boundary conditions. The shaking of the container implies that , with and the amplitude and frequency of energy injection in the model. At the top, s(t), we consider a free surface, and thus the kinematic boundary conditions are given by


Furthermore, the stress at the bottom and top of the granular media are needed. The free surface at the top is straightforward: . At the bottom, on the other hand, we divide the stress contribution in two: mean () and fluctuating () terms, where the mean term is straightforward: , with the total mass of the system, ; and the area of the base of the container, .

For the fluctuating part of the stress, , we first consider the force applied to the granular medium by the moving bottom:


with the mass being pushed by the bottom wall. In order to obtain , let us consider a moving platform of surface area pushing an ideal, incompressible gas of density , in analogy to the moving box and the low density region observed in our system. Accordingly, the mass pushed by the box in time is given by . Notice that this is valid for high , where gravity effects on the dynamics of the particles can be ignored. Integrating, we directly get that . Substituting in (11):


Notice that we have naturally obtained as the amplitude of the force applied by the oscillating bottom, further suggesting that the shaking strength is the relevant parameter for the system in the high limit. It then follows, from (11), that:


Finally, substituting the stress boundary values in (8), we obtain:


3.3 Height averaging

The remaining two integrals in (14) involve the velocity profile, , which varies in the vertical direction. In order to solve these integrals we height-average, that is, for a given quantity , we consider its average value


Notice that, from the first integral in (14), would correspond to . Thus, before applying (15), we express the integral as a total time derivative. Considering that the boundaries are time dependent, it becomes necessary to use Leibniz integration rule, and thus the first integral in (14) can be expressed as


Analogously, the second integral in (14) becomes, after using (10),


Substituting (16) and (17) in (14), and using (15), we finally obtain:

Figure 7: From left to right, snapshots from simulations showing an LFO period, at phases , , , and ; the corresponding time averaged density profile, a representation of the two phases approximation made for the continuum equations, and finally a schematic representation of the model. The dashed line shows the position of the centre of mass, , which in the model corresponds to the position of the interface between the two phases, , which also corresponds to the position of the mass of a forced harmonic oscillator.

Based on the behaviour observed in simulations, we now assume that the high density region is incompressible. This implies that , as also that the velocity of the continuum media at the interface position is equivalent to the velocity of the interface, and hence to the velocity of the surface, that is, . Thus (3.3) becomes


Furthermore, we now use the fact that . Thus, substituting and dividing by , we obtain


3.4 First order approximations

It now becomes relevant to consider the relative importance of each of these terms, in the region of phase space where simulations show that LFOs are present, that is, for . First, we consider that , a condition that holds better for and low , as shown in Figure 2b. On the other hand, , as can be seen from Figure 4. Furthermore, we measure from simulations that and , meaning that the dynamics of the LFOs are considerably lower than the typical velocity of grain diameters per gravity timescale, as can also be deduced by the previously obtained frequencies . Let us also notice that , again, from Figure 2b. Finally, from simulations we obtain that , and . Taking into account all these considerations, it becomes straightforward to see that the first term in (3.3) is , the second term is at most , the third is then , and the fourth term is . Moreover, all terms on the right side are . Thus, disregarding small terms in (3.3), after dividing by , we obtain


where we have defined the mass of the solid region per unit base area , , and the equivalent of the gaseous region, . Equation (21) corresponds to a forced harmonic oscillator equation of the form:


with natural frequency


amplitude of forcing , and constant .

3.5 Model and simulations comparison

We have shown that, considering Cauchy’s equations for continuum media, and making assumptions in concordance to the observed granular Leidenfrost state, the system becomes equivalent to a simple forced harmonic oscillator, expressed by (22). In this case, is the displacement of the centre of mass around the equilibrium position at , the natural frequency of the system, and and the amplitude and frequency of the forcing. The analogy of the forcing with the granular column is straightforward: and would be equivalent to and , respectively. Furthermore, we choose to correspond to , in order to directly compare with previous measurements.

Notice that the natural frequency does not explicitly depend on the forcing frequency , as can be seen in (23). The implicit dependence comes from the variation of and with , as observed in simulations, where, for fixed , increases with , giving the correct inverse proportionality of with . Therefore, in order to obtain a frequency from the model, only and need to be specified, which we measure from simulations.

Figure 8: Low-frequency oscillation frequencies , as a function of the dimensionless shaking strength , for different box oscillation amplitudes: (a), (b) and (c). All systems have , . Simulation (black) corresponds to frequencies obtained from fast-Fourier transform of the simulation data, while continuum and thermodynamic/kinetic theory data points (blue and red, respectively) are obtained from models presented in sections 3 and 4, respectively, using data acquired from simulations.

Both quantities can be obtained from , the density in the granular column as a function of height. In order to obtain an accurate average, we consider , which makes all profiles directly comparable. This is analogous, in the model, to centering the profiles at the interface between the two distinct regions. It is then straightforward to compute as the average value of the density for . On the other hand, we take as the total mass for , taking care not to count particles that are in free flight above the solid region, as they do not have influence on the oscillator dynamics. This implies that although the center of our profiles is , is not equal to .

The comparison between the frequencies obtained in simulations and from the model is presented in Figures 8a-c. For low amplitudes, , and high frequencies, , the agreement between the frequencies is within the error bars. For lower , or higher amplitudes, the assumption of two distinct phases, as also the approximation of , become less justified, resulting in the model consistently underpredicting the frequencies, with more than disagreement at the point of the bouncing bed-Leidenfrost transition for . We believe that the prediction could be improved by considering more complex density profiles, as also by including terms of lower orders, although this exceeds the scope of our work. In general terms, the resulting one-dimensional model turns out to be a remarkable well approximation for high and low , showing that this many-particle, out-of-equilibrium system actually behaves as a regular forced harmonic oscillator when confined in a column, in the corresponding energy injection region.

4 Thermodynamic model

Remarkably, it is possible to obtain another accurate expression for using a completely different approach, considering basic concepts from thermodynamics. Assuming a spring-like behaviour, the natural frequency of our medium is given by , with the stiffness constant of the spring-like medium, and the mass sustained by the spring. We also know that , with the bulk modulus, the area of the spring, and it’s height at rest. Assuming an adiabatic ideal gas, it is possible to relate the bulk modulus with the pressure, , with the adiabatic index. Notice that the ideal gas approximation is being used only for the gaseous region of the Leidenfrost regime, where densities are low and no significant correlation of the particles is observed. We then obtain:

The pressure is taken as the force caused by the solid mass the spring sustains, , divided by the area of the container: . Thus, we finally reach:


This significantly simple expression is remarkably accurate when compared with simulation measurements. Figures 8a-c also show for this model, taking to be the same as in the previous section, . The adiabatic index is considered as , the theoretical value for an ideal monoatomic gas. The agreement is again within error bars for high frequencies, and deviates considerably for lower frequencies, except for the case, where low frequencies are also captured. Relating the two obtained LFOs frequencies, equations (23) and (24), one obtains ; the interpretation of this result remains a challenge.

5 Conclusions

A vertically vibrated bed of grains presents low-frequency oscillations (LFOs) due to the decoupling of the driving frequency and the dynamics of a high density region suspended by a lower density one. The relevance of these oscillations increases as the distinction between the two densities increases, that is, proportional to the frequency and inversely to the amplitude of oscillation of the system container. The LFO frequencies are inversely proportional to the driving frequency, and follow a common power law for a range of amplitudes. The amplitude of the oscillations, on the other hand, increases in an almost linear way with the frequency.

Event-driven simulations give an overall excellent qualitative and quantitative agreement with experiments and soft-particle simulations done in wider systems, although they show discrepancies in some critical transition values. We remark that the hard-sphere approximation can be meaningful even in systems with very high-density regions, as present in the Leidenfrost state. The considerable speed advantage makes it extremely useful, and sometimes the only means to systematically study high dimensional parameter spaces.

Starting from Cauchy’s equations for conservation of mass and momentum, integrating in the vertical direction and assuming two distinct low and high constant density regions, it is possible to reproduce the frequency behaviour observed in simulations. That is, a forced harmonic oscillator, with the natural frequency proportional to the ratio of the densities. This simple model is able to predict the natural LFO frequency for high excitation frequencies, where in fact the two phases are well separated. The non-linear terms, discarded in our analysis, should provide the necessary corrections for lower frequencies, as well as the consideration of a more realistic density profile. A second approach, using thermodynamic arguments, also gives a remarkably accurate expression for the frequencies, although in this case just a simple mass-spring system behaviour was assumed. The quantitative agreement of both models is nevertheless remarkable, taking into account the low number of particles involved, and the presence of very high and low density regions.

Further insight could be gained by appropriately coarse graining the granular medium in order to obtain stress fields, which would directly relate both models. A point of interest, not studied here, is how well do kinetic theory predictions hold in such a system, taking into account the reduced container size, the small number of particles and the presence of considerably different densities. Current work is being done on verifying the consistency of macroscopic fields obtained by theoretical arguments and coarse-grained simulational data.

We suggest that LFOs, here shown to be ubiquitous to vertically vibrated density inverted systems, could play a fundamental role in the Leidenfrost-convection transition. More specifically, LFOs could be the primary source of density fluctuations observed before convection is triggered, when one region of the system oscillates at a different phase than another. Understanding this will need further simulation and experimental research.

6 Acknowledgements

The authors wish to acknowledge Thomas Weinhart for a careful revision and important suggestions regarding Section 3. This work was financially supported by the NWO-STW VICI grant 10828.



  • [1] J. S. Olafsen and J. S. Urbach. Clustering, order, and collapse in a driven granular monolayer. Physical Review Letters, 81(20):4369, 1998.
  • [2] P. Melby, F. Vega Reyes, A. Prevost, R. Robertson, P. Kumar, D. A. Egolf, and J. S. Urbach. The dynamics of thin vibrated granular layers. Journal of Physics: Condensed Matter, 17(24):S2689–S2704, 2005.
  • [3] F. Melo, P. Umbanhowar, and H. L. Swinney. Transition to parametric wave patterns in a vertically oscillated granular layer. Physical Review Letters, 72(1):172, 1994.
  • [4] S. Luding, E. Clément, J. Rajchenbach, and J. Duran. Simulations of pattern formation in vibrated granular media. Europhysics Letters (EPL), 36(4):247–252, 1996.
  • [5] S. G. K. Tennakoon and R. P. Behringer. Vertical and horizontal vibration of granular materials: Coulomb friction and a novel switching state. Physical Review Letters, 81(4):794, 1998.
  • [6] M. Medved. Connections between response modes in a horizontally driven granular material. Physical Review E, 65(2):021305, 2002.
  • [7] K. Ahmad and I. J. Smalley. Observation of particle segregation in vibrated granular systems. Powder Technology, 8(1–2):69–75, 1973.
  • [8] A. Kudrolli. Size separation in vibrated granular matter. Reports on Progress in Physics, 67(3):209–247, 2004.
  • [9] S. Luding and H. J. Herrmann. Cluster-growth in freely cooling granular media. Chaos: An Interdisciplinary Journal of Nonlinear Science, 9(3):673, 1999.
  • [10] N. Rivas, S. Ponce, B. Gallet, D. Risso, R. Soto, P. Cordero, and N. Mujica. Sudden chain energy transfer events in vibrated granular media. Physical Review Letters, 106(8):088001, 2011.
  • [11] B. Thomas, M. O. Mason, Y. A. Liu, and A. M. Squires. Identifying states in shallow vibrated beds. Powder Technology, 57(4):267–280, 1989.
  • [12] S. Douady, S. Fauve, and C. Laroche. Subharmonic instabilities and defects in a granular layer under vertical vibrations. Europhysics Letters (EPL), 8(7):621–627, 1989.
  • [13] C. R Wassgren, C. E Brennen, and M. L Hunt. Vertical vibration of a deep bed of granular material in a container. Journal of Applied Mechanics, 1996.
  • [14] B. Meerson, T. Pöschel, and Y. Bromberg. Close-packed floating clusters: Granular hydrodynamics beyond the freezing point? Physical Review Letters, 91(2):024301, 2003.
  • [15] J. G. Leidenfrost. De aquae communis nonnullis qualitatibus tractatus. 1756.
  • [16] P. Eshuis, K. van der Weele, D. van der Meer, and D. Lohse. Granular leidenfrost effect: Experiment and theory of floating particle clusters. Physical Review Letters, 95(25), 2005.
  • [17] P. Eshuis, K. van der Weele, D. van der Meer, R. Bos, and D. Lohse. Phase diagram of vertically shaken granular matter. Physics of Fluids, 19(12):123301, 2007.
  • [18] B. D. Lubachevsky. How to simulate billiards and similar systems. Journal of Computational Physics, 94(2):255–283, 1991.
  • [19] S. Luding. Granular materials under vibration: Simulations of rotating spheres. Physical Review E, 52(4):4442–4457, 1995.
  • [20] Stefan Luding and Sean McNamara. How to handle the inelastic collapse of a dissipative hard-sphere gas with the TC model. Granular Matter, 1(3):113–128, 1998.
  • [21] M. G. Clerc, P. Cordero, J. Dunstan, K. Huff, N. Mujica, D. Risso, and G. Varas. Liquid-solid-like transition in quasi-one-dimensional driven granular media. Nat Phys, 4(3):249–254, 2008.
  • [22] Sean McNamara and Stefan Luding. Energy flows in vibrated granular media. Physical Review E, 1998.
  • [23] J. G. Leidenfrost. On the fixation of water in diverse fire. International Journal of Heat and Mass Transfer, 9(11):1153–1166, 1966.
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