Resonance effects on the dynamics of dense granular beds: achieving optimal energy transfer in vibrated granular systems

Using a combination of experimental techniques and discrete particle method simulations, we investigate the resonant behaviour of a dense, vibrated granular system. We demonstrate that a bed of particles driven by a vibrating plate may exhibit marked differences in its internal energy dependent on the specific frequency at which it is driven, even if the energy corresponding to the oscillations driving the system is held constant and the acceleration provided by the base remains consistently significantly higher than the gravitational acceleration, g. We show that these differences in the efficiency of energy transfer to the granular system can be explained by the existence of resonances between the bed’s bulk motion and that of the oscillating plate driving the system. We systematically study the dependency of the observed resonant behaviour on the system’s main, controllable parameters and, based on the results obtained, propose a simple empirical model capable of determining, for a given system, the points in parameter space for which optimal energy transfer may be achieved.


Introduction
One of the major challenges facing researchers and industrialists concerned with the investigation of vibrated granular media is the choice of appropriate control parameters which may be used to accurately characterize the state and behaviour of a granular system. Indeed, without a knowledge of the correct control parameters, it is effectively impossible to reliably predict the dynamical behaviours of a granular system subjected to vibrational excitation. Numerous, widely-used industrial processes are reliant on the vibration of granular beds [1][2][3][4], including processes relevant to highly lucrative (e.g. mining and pharmaceutical) industries, as well as those with applications of considerable environmental impact, such as the recycling of hazardous electronic waste [5]. An improved predictive knowledge of the relevant bed dynamics has the potential to significantly increase the efficiency and reduce the cost associated with these processes through improvements in the design and operation of the apparatus involved.
The energy with which a continuously, harmonically vibrated granular bed is excited is determined by the frequency, f, and peak amplitude, A, of the oscillations applied to the system. However, considerable effort has been made within the field of granular dynamics to determine a single, simplified parameter which may univocally describe the steady state of a dynamic granular system. There exist three main parameters which are typically used to this end: the peak vibration velocity, ω = v A, the dimensionless peak acceleration, Γ = ω A g 2 [6], and the 'shaking strength parameter', = = ω S v gd A gd 2 22 [7-9], a dimensionless measure of the peak kinetic energy associated with the driving vibrations. Here, ω π = f 2 , where f is the oscillatory frequency, while d and g denote the diameter of a particle and the gravitational acceleration at the Earth's surface respectively. It is widely suggested in the literature that many key properties possessed by a granular system-for example its mean potential and kinetic energies, its granular temperature [10], and its typical packing density-may be expected to follow simple, monotonic relationships with one or more of these parameters [11][12][13][14]. However, it has recently ). It is suggested that these striking results may in fact be related to the existence of resonances within the system, making the work of [20,21] even more relevant to the current study, as will become apparent later in this paper.
In the current work, we demonstrate that the lack of a simple dependence on v S , or Γ of a harmonically vibrated system's internal energy, as observed in the shallow, strongly-driven beds explored in [15], persists also in deeper, more dissipative and more weakly-fluidized beds. Indeed, it is demonstrated that such systems are in fact found to exhibit an even more complex relationship between the vibrational parameters with which the system is driven and its resultant dynamical properties. We attempt to provide a coherent explanation of the physical mechanism responsible for this newly observed behaviour. Specifically, we posit that the efficiency of energy transfer from a vibrating container to a bed of particles housed within depends sensitively on the existence of resonances between the motion of the system's oscillating base and that of the bed itself. Based on this proposed mechanism, we present a simple theoretical framework capable of predicting, for a given, quasione-dimensional (1D) system, the points in phase space for which one may achieve maximal energy transfer from the vibrating base to the granular system driven thereby.
The existence of resonant behaviour within granular systems has long been known; in 1978, Pieranski [22] demonstrated experimentally the presence of Feigenbaum-like period-doubling bifurcations [23] for the simple case of a single particle 'jumping' on a vibrating surface. The period-doubling behaviours of a single ball bouncing on a plate were further investigated, both experimentally and through simulations, by Tufillaro et al [24], who determined, for a given system, the regions of f-A parameter space for which different periodic and chaotic dynamics may be expected to occur. In this paper, it was additionally suggested that other driven, dissipative systems-such as granular beds-may exhibit behaviours similar to those shown for the singleparticle case explored; in particular, it was noted that one may expect the phase at which period-doubling occurs in more complex systems to depend only on a single parameter-the dissipation rate-as is reported to be the case in the simple model presented.
The single-particle system also formed the focus of work by Mehta and Luck [25]. In this instance, however, base-particle interactions were modelled as perfectly inelastic. Despite the simplicity of the system investigated, a rich phenomenology was nonetheless observed, with the existence of a period-doubling sequence whichunlike for the case of elastic particles, as investigated by Pieranski [22]-becomes abruptly interrupted as the system enters a 'locking region', wherein the inelastic particle will only collide with the moving base during the 'absorbing', as opposed to 'transmitting', phases of its sinusoidal motion, clearly resulting in a periodic freeflight pattern.
Work by Melo et al [26] extended the analysis of period-doubling behaviours to the multi-particle case, exploring the manner in which these behaviours, alongside the presence of a standing-wave instability [27,28], influence pattern formation in a granular bed. They show that the flight times of the beds explored can be closely approximated using the model of a single, perfectly inelastic sphere due to the rapid, near-total dissipation of energy arising from the great number of internal particle collisions within the granulate.
Comparison was once more drawn between the behaviour of a vibrated granular bed and that of a single inelastic ball in the 1997 study of van Doorn and Behringer [29], this time modelling the motion of a thin, avalanching surface layer atop a deep, weakly excited, granular bed exhibiting convective flow. It is found that the centre-of-mass trajectory for this thin surface layer as it detaches from and rejoins the bulk of the bed shows reasonable agreement with predictions based on the inelastic bouncing ball model (IBBM).
A subsequent article by Pastor et al [30], however, calls into question the accuracy with which the IBBM may predict the dynamics of a granular system comprising multiple particles, focussing specifically on the temporal bifurcations corresponding to the free flight, τ f , of the bed under investigation. It is shown that although agreement between the model and experimental results is strong for small accelerations (Γ ≲ 2.7), the model fails for larger driving strengths, systematically overestimating free-flight times and inaccurately predicting the positions of bifurcation points in τ f -Γ parameter space. Several explanations for the observed deviations are proposed, including the presence of interstitial air effects, wall friction and, perhaps most prominently, the fact that the model does not account for the finite duration of collisions between bed and plate.
Of particular relevance to our current investigation is the work of Sánchez et al [31] which shows that, under the right conditions (namely a large particle number, N, and a high particle density, η), the resonant response of a granular damper is effectively independent of the material properties possessed by the individual particles of which it is formed. It is hypothesized that this universal response occurs due to the inelastic collapse [32,33] exhibited by a granular system under the aforementioned bed conditions, meaning that the bed will possess effectively zero elasticity, irrespective of the specific restitution (ε) and frictional (μ) coefficients of the particles involved (at least so long as μ > 0 and/or ε < 1). Further discussion of this study, and its relation to the findings presented here, may be found in section 2.3.
Further work by Sánchez et al [34] investigates the variation in the dynamical behaviours of a particle damper as the frequency of its oscillation is altered. Their study details a transition, with increasing oscillatory frequency, between several differing dynamical states: a periodic motion in which the particle bed simply follows the motion of the oscillating system without detaching from the base; a 'transition region' where the bed switches between periodic and quasi-periodic motions; a region of chaotic motion; an 'optimal damping' region where the bed impacts both the lower and upper boundaries of the system; and finally, a second chaotic region, whose behaviour is, interestingly, distinct from the first. The maximization of energy dissipation in the 'optimal damping' region corresponding to the system's primary resonant frequency is attributed to the fact that in this state the bed is highly compacted-inducing a high collision rate and hence rapid collisional dissipationcombined with the fact that, at the points of collision between bed and base or bed and ceiling, the two interacting objects possess opposing velocity vectors, therefore increasing relative speed and hence further increasing dissipation [35].
Despite the considerable body of research investigating resonant behaviours in vibrated granular systems, this current work represents, to the best of the authors' knowledge, the first systematic study of the effect of these resonances on the internal energy possessed by a continuously vibrated bed of particles, and the efficiency with which this transfer may occur. and d 16 being used to allow the analysis of both fully three-dimensional (3D) beds as well as quasi-1D 'columns' of particles. For each of these geometries, the height of the container in which the particles are housed, L z , is equal to 50 particle diameters, a value adequately large to minimize the probability of particle collisions with the system's upper boundary in all systems explored. A simple, schematic representation of the systems described above may be seen in figure 1.

System details
Since the imaging technique used to acquire information from the experimental system utilizes highlypenetrating γ-radiation as opposed to electromagnetic radiation in the visible spectrum, the choice was made to construct the entire experimental cell from stainless steel. The presence of a conductive base and sidewalls prevents the build-up of static charge within the bed, meaning that triboelectric effects on the system may be safely neglected [36], particularly when we consider also the considerable size of the particles used. This substantial particle size also means that effects due to interstitial air within the system may also be ignored [37]. The use of a steel-walled system is also beneficial in that the high coefficient of restitution for interactions between steel particles and flat steel surfaces [38] will act to frustrate unwanted sidewall effects, such as walldriven thermal convection [39] and the heat-sink-like behaviour observed when highly dissipative walls are introduced to a system [40].
Sinusoidal vibrations are applied to the systems under investigation using an LDS V721 electrodynamic shaker. In order to assess the generality of our observations, in addition to the variation of N, we also explore a wide range of driving frequencies and amplitudes, specifically ∈ f (5, 80) Hz and ∈ A (0.6, 8.0) mm, thus providing a range of peak driving energies and accelerations ∈ S (0.9, 11.5) and Γ ∈ (1.0, 23.9), respectively.

Data acquisition-positron emission particle tracking (PEPT)
Experimental data are acquired from the system described above using PEPT. PEPT is a non-invasive technique which uses a dual-headed γ-camera to track and record the motion of a single, radioactively labelled particle within a bed of other, similar particles. In order to produce a suitable radioactive tracer particle for our current system, a single steel beadphysically identical to all others used in experiment-is irradiated with a proton beam from the Birmingham Figure 1. Schematic diagrams (not to scale) illustrating a quasi-1D, 'columnar' system (left) and a fully 3D system of width = L 16 x y , particle diameters (centre). Shown also (right) is a simplified illustration of the manner in which reconstructed γ-ray paths may be used to triangulate the position of a radioactive tracer particle within a given system. Although in reality the tracer particles used in experiment are physically identical to all other particles within the bed, in the above image, for the sake of clarity, the tracer is highlighted in red. cyclotron, converting 56 Fe atoms within the particle to β + -emitting 55 Co. The positrons emitted by these cobalt isotopes will rapidly annihilate with electrons in the dense medium of the tracer, producing pairs of 511 keV γrays whose trajectories are separated by 180°. If multiple such pairs of γ-rays are detected by the opposing heads of the γ-camera, the paths taken by the photons can be reconstructed, allowing the position of the tracer within the bed to be triangulated, as illustrated in figure 1. For a stationary particle, this triangulation can be performed to an arbitrary accuracy, while for a moving particle, a rapid series of these location events typically allows the motion of a tracer particle to be recorded with a spatial resolution on the millimetre scale, and a temporal resolution of the order of milliseconds [41]. The use of high energy electromagnetic radiation allows particle motion to be tracked even deep within the interior of large, 3D beds and, as is the case here, within highly dense and opaque systems. For systems in a non-equilibrium steady state (NESS), the motion of the single tracer particle used can be utilized to extract various important parameters pertaining to the dynamics of the system as a whole [42]. The quantities which may be obtained using PEPT include, but are by no means limited to, packing densities [42], granular temperatures [43] and velocities [44], as well as the spatial distributions thereof, mean squared displacement and diffusive behaviours [45] and even segregation in binary [46] and polydisperse systems [47]. In the current paper, however, we are predominantly concerned with our systems' time-averaged vertical centre of mass positions, h, which provide an accurate measure of the mean gravitational potential energy (E P ) possessed by a granular bed. It is perhaps worth noting that for vertically vibrated systems under the influence of gravity, such as those discussed here, the spatial distribution of energy and particle density within a given bed is likely to exhibit significant inhomogeneity [48][49][50]. Nonetheless, the time-averaged value h should still provide a reasonable measure of the typical potential energy possessed by a system and, most importantly, its variation with the various system parameters explored. For a 3D system, the mean centre of mass position of a system can be extracted from PEPT data simply by determining the long-time average of the tracer particle's vertical (z) position within the system as it moves through the vibrated bed during a run. Due to the principle of ergodicity, such a value is directly equivalent to the typical centre of mass position for the system as a whole. In order to obtain a value of h-and hence E P -for a given, 3D system, data are acquired over a period = t 4800 s run from a bed of fixed N vibrated at a constant, fixed frequency, f, and amplitude, A. Each data point shown in figures 3(a), (b) corresponds to such a data set. The steady-state nature of the system in all instances is ensured as follows: firstly, each system is vibrated for a period of 500 s before data acquisition begins, in order to allow the system to achieve a NESS. Once data have been acquired, each 4800 s trace is divided into a series of overlapping 300 s segments; individual values of h are then acquired for each segment, and consistency between all values obtained is ensured. This process is performed for all fully 3D systems investigated.
For a 1D column of N particles, meanwhile, where the tracer is not free to explore the entire system, N individual runs, each of duration 900 s, are conducted, with the tracer particle taking a different position within the system for each run. The average centre of mass height can then be determined through an appropriate averaging of all N data sets acquired. In this situation, where data are individually recorded for each particle position within the system and these positions are not interchangeable, the above-described assumption of ergodicity is, clearly, not required.

Discrete particle simulations-MercuryDPM
The use of discrete particle method (DPM) simulations [51] in addition to experimental work allows a larger parameter space to be explored, and a more thorough exploration of this parameter space, than would be feasible using experiments alone. In the current study, the use of simulations allows us not only to closely reproduce the behaviours of the experimental system under investigation, but to perform in-depth parameter studies capable of elucidating the underlying trends necessary to the understanding and prediction of the phenomena observed.
Simulations are acquired using the MercuryDPM software package [52][53][54][55] developed at the University of Twente. The simulations use a frictional spring-dashpot model with a linear elastic and a linear dissipative contribution for both the normal and tangential forces acting between colliding particles. The model used allows the coefficient of restitution, ε, to be directly defined by the user. This parameters, in combination with a similarly user-defined contact time, t c , may then be used to determine the relevant normal spring constant, k n , and damping coefficient, γ n , as: and γ γ = t n . With a known k n t , and γ n t , , the spring-dashpot normal force can be calculated as: ij n n ij n ij n ij n and the corresponding tangential force as: , is the relative normal or tangential velocity component and δ ij t is the elastic tangential displacement (see [56]). A static yield criterion is also applied such that the magnitude of δ ij t is truncated as necessary to satisfy , with μ the relevant frictional coefficient. In order to recreate as closely as possible the experimental system described in section 2.1, the particle diameter, d, particle density, ρ, system size, L x y z , , , particle number, N, driving frequency, f, and amplitude, A, are all implemented as their precise experimental values. When direct comparison between experiment and simulation is required, experimentally determined values of the normal coefficient of restitution, ε, are also implemented [38]. However, one of the major advantages of the inclusion of DPM simulations in our study is the ability to additionally explore the effect of variations in ε without altering any other particle properties-a near impossible feat in practical experiments; once the simulations have been 'calibrated' (i.e. once agreement between experiment and simulation has been demonstrated) using the known, experimental ε values, it is then possible to investigate the specific influence of ε by changing this particular parameter whilst holding all other variables fixed at the values known to successfully reproduce the experimental situation. The frictional coefficient is taken as μ = 0.1 [38], a value previously shown to provide accurate results in similar systems [57], and the contact time, t c , is implemented as 10 −6 s. The appropriateness of our chosen t c value was confirmed by performing repeated runs using values of t c a factor of 10 above and below this value and ensuring consistency in the resultant data.
Since, unlike for results acquired using PEPT, simulational data provide information for all particles within the system, the bed's vertical centre of mass position can, for a given time step, be easily calculated from the known spatial positions of all particles within the system. The time-averaged value of h is then taken simply as the mean value across all timesteps.
Examples of the strong agreement produced between experiment and simulation may be seen in figure 4. For the case of the 20 particle system, the close agreement observed is perhaps unsurprising as the large particle density and correspondingly high collision rate mean that the bed may be expected to behave as a single, perfectly inelastic body, resulting in relatively simple dynamics which are not dependent on the detailed interaction parameters of particles [31]. Nonetheless, the continued close correspondence observed between experiment and simulation for the (lower density) 15 particle case, where the system is unlikely to experience inelastic collapse, strongly suggests that the contact models used, and the relevant parameters chosen, are indeed adequate to emulate the experimental situation. Further support for the validity of the simulations used is provided by the consistent accuracy with which this model has been shown to capture the behaviour of similar experimental systems in previous studies, where the strongly driven, low density beds explored exist far from the limit of inelastic collapse [57][58][59][60][61].
The use of simulations additionally facilitates the confirmation of the fact that the systems explored may indeed be expected to exhibit ergodicity. Such a check is performed for a given system by comparing the density and temperature fields produced through an ensemble average of all particles within a system over a period of time Δt to those calculated from the motion of a single particle within an identical system over a period Δ N t · . Comparison was also drawn between the calculated vertical centres of mass for both cases, as this parameter forms a major area of interest for the present study. Checks such as these were conducted for a range of system variables spanning the relevant parameter space explored, and using various Δt values. The motion of a single particle was found to accurately capture the system's dynamics and physical state for all systems tested, lending credence to the validity of our PEPT results.
For the interested reader, further details of the implemented contact law may be found in our reference [56] or, alternatively, in the supplemental material of our previous work [58]. An overview of contact laws in general may be found in [62].

Results and discussion
We begin by discussing the general form of the relationship between the system's driving parameters and the amount of energy successfully transferred from the vibrating base to the granular bed itself. We initially quantify this energy in terms of the gravitational potential energy, E P , possessed by the granulate, as this quantity can be easily and, more importantly, accurately extracted both from simulational and experimental data sets. The trends discussed here are, however, qualitatively similar to those observed in the total energy of the system, E T , as may be seen from figure 2. This image clearly demonstrates the strong similarities between the forms of E P and E T as well as, indeed, the kinetic energy, E K -all three quantities follow a similar series of local maxima and minima, the positions of these extrema lying approximately coincident for the three parameters. The agreement is observed to break down only in the low-f (low-Γ) limit, as may well be expected from our previous work [15]. It is worth noting that the agreement exemplified in figure 2 is in fact general to all results presented.
In order to test the efficiency of energy transfer from base to bed, a system of fixed N, N L , ρ, d, ε and L x y z , , , is driven with a fixed, constant energy input S. A series of individual data sets are then taken maintaining the given value of S, but varying the combination of frequency (f) and amplitude (A) used to produce this specific value. The average energy attained by the granular system for each of these data sets can then be plotted as a function of f, A or Γ, allowing the dependence of energy transfer on these parameters to be assessed. Exemplary images may be seen in figure 3.
In the case of relatively dilute systems-i.e. those possessing relatively small bed heights (N L ) and comprising relatively elastic particles-one simply observes a single maximum (see figure 3(a)) in the energy transferred to the vibrated granulate, the existence of which is explained in detail in [15]. For deeper and/or less elastic systems, however, one observes instead a series of local maxima (see figure 3(b)), seemingly 'superimposed' on the general form expected for the lower-density case. As is clear from the images shown, these local extrema are observable both in 1D and 3D systems.
In order to allow a simpler, more reliable analysis, we wish to isolate, as far as possible, the 'new' behaviour observable in the denser systems. Thus, the current study focuses largely on the column geometry, which will act to suppress the amplitude-dependence of the energy transfer previously discussed in [15] by ensuring that only a single particle may contact the base, thus limiting the variation in the mean particle-base collision rate. The efficacy of this manner of isolating the newly observed local maxima can be demonstrated through a comparison of panels (b) and (c) of figure 3 which show, respectively, a 3D system comprising 10 layers of 256 particles and a 10 particle, 1D column of identical particles. Figure 4 shows images corresponding to a columnar system of 15 particles, as well as a similarly excited column of 20 particles, the latter showing data for a greater range of driving frequencies, f. Immediately notable from this second image is the approximately uniform spacing of the local maxima, suggestive of a resonant behaviour [22, 24-26, 29-31, 34] within the system. Support for this possibility may be seen in figure 5, which shows the fast Fourier transforms (FFTs) of the bed's centre of mass motion (i.e. the transformation of the bed's centre of mass trajectory over time) for data sets corresponding to each of the local maxima.
The power spectra shown in figure 5 all exhibit a series of distinct, delta-like peaks at frequencies ( ) closely corresponding to integer multiples of the driving frequency (f) associated with the first peak in the E-f plot of figure 4. Although peaks in the FFT of the bed's centre-of-mass motion appear also for driving frequencies which do not match the system's resonant frequency, as may be expected [34], these peaks are typically smaller in amplitude, as illustrated in figure 6. The increased amplitude observed for f = 13.4 Hz and its multiples, relative to other driving frequencies, demonstrates a greater transfer of energy to the system at these values, thus providing an indication that we are indeed observing resonant behaviour. Specifically, these results strongly imply that the system depicted possesses a characteristic resonant frequency of ∼13.4 Hz, explaining the marked increase in the energy transferred to the bed when the driving oscillations approximately match this frequency or one of its harmonics 4 .
We now attempt to provide a cogent explanation of the origin of the apparent resonant behaviour. The sharp peaks observable in the FFTs of the system's centre-of-mass motion demonstrate the system to be in the 'bouncing bed' state [9], wherein the bulk of the particles in the bed undergo coherent motion, behaving similarly to a single body repeatedly impacting upon the vibrating base plate in between periods of free flight. We propose that the observed peaks in the system's total internal energy, E T , correspond to the cases in which the period of the base's vibrational motion is congruent with the typical free flight time of the bed. This hypothesis would also explain the absence of these peaks in the case of shallower, less dissipative and hence more dilute systems where particle motion is fluid-like and randomized. Despite the simplicity of this proposed mechanism, the prediction of the frequency values for which these resonances will occur is not a trivial matter, as the free flight time of the bed is sensitively dependent on a variety of factors, most notably the number of particles (or particle layers) within the system, the material properties (in particular the restitution coefficient, ε) of the particles and the strength with which the system is driven, S. The influence of these various parameters on the system's resonant behaviour is clearly illustrated in figure 7.
The free flight time of a system of particles whose dynamics match those described above may be simply determined by measuring, over a large number of oscillation cycles, the average period between the detachment of the bottom-most particle from the system's base and its subsequent recollision. However, such a determination of the free-flight time, τ f , ignores the possibility of period-doubling bifurcations. In order to avoid this potential complication, we analyse the distribution of free flight times from which our average t f value is calculated, ensuring that said distribution is unimodal-i.e. free of any additional peaks indicative of perioddoubling behaviour.
An initial visual inspection of the data shown in figure 7 shows a qualitative agreement with the basic expectations of our above rationale: as the energy provided to the system is increased, one would expect the bed to be 'launched' from the base at a greater velocity, resulting in an increased flight time and hence a reduced  resonant frequency, as is indeed shown in image (a). As particles become less elastic, however, one would expect a decrease in flight time and a correspondingly increased resonant frequency, in agreement with the general trend shown in image (c). Similarly, as the particle number, N, is increased, the accordingly greater collision rate between particles will act to increase dissipation within the system, once again leading to a reduction in flight time and an augmented resonant frequency (see image (b)).
Having provided compelling evidence of the physical mechanism underlying the observed variation in systems' internal energy with driving frequency, we now attempt to construct a simple framework capable of predicting the points in parameter space for which one may expect resonances to occur, i.e. where one might achieve optimal energy transfer from the driving wall to the granulate itself. We first consider the limiting case ε → 0, i.e. the inelastic limit; in this case, any energy input to the granular system would be entirely and immediately dissipated, meaning that we would simply expect the system's behaviour to reduce to that of a perfectly inelastic body, a situation which has been explored extensively in the existing literature [22, 24-26, 29, 30, 63]. This allows us to establish an appropriate minimum free flight time, τ f 0 , which can, using simple kinematics, be determined from the system's peak base velocity, , the large number of particle contacts-and hence collisions-will rapidly dissipate any energy input, meaning that the system may once again be expected to approximate the behaviour of a single, perfectly inelastic mass. Finally, it must also hold true that as the energy input, and hence velocity, tend toward zero, the free flight time approaches zero also. Although the assumption that τ f 0 will vary linearly with the peak base velocity, v, is somewhat over-simplistic (it has been suggested in previous studies [63][64][65] that the velocity at detaching, v d , is, for the case of a perfectly inelastic ball, a more appropriate parameter) this direct proportionality is nonetheless well supported by our results, as demonstrated in figure 8.
With the above restrictions in mind, we find that all data sets acquired are well described by the functional form: where α is a constant whose precise value may be expected to depend on the various system properties not explicitly discussed here, such as the particle-wall coefficient of restitution, as well as inter-particle and particlewall friction. It is also perhaps worth noting that an exponent with an − N ( 1)dependence-as opposed to a simple N dependence-is implemented as it is the number of particle contacts, as opposed to the number of particles, which will determine the energy lost through collisional dissipation. It should be noted that equation (5) calculates τ f using the peak base velocity, v, as opposed to v d , as v, like N, is a well-defined, controllable system parameter, whereas v d must either be directly measured (complicating the use of the formula) or determined theoretically, thus introducing additional assumptions regarding the system's state.
The phenomenological formula presented above is found to demonstrate a pleasing agreement with our simulational results over a range of N, S (v) and ε values, as illustrated in figures 9 and 10. In all cases shown, the theoretical fits share a constant α = 12.8. Using this fixed α value, theoretical and simulated results are found to agree to within 10% for all cases investigated, and to within 5% in the majority (≳90%) of cases. Furthering the above discussion relating to the linear trend between τ f 0 and v employed in equation (5), it is interesting to note the continued agreement across a range of driving parameters (see figure 10), even at relatively high S, Γ and f values where one might expect our simplified assumption τ ∝ v f to break down. Using the above formulation, it is now trivial to predict the frequencies, f r , at which resonances, and hence local maxima in energy transfer, may occur from the relation = τ f r n n f , (where n is an integer). Although, at present, the α term in (5) remains an unknown fitting parameter, in practice α may be simply determined for a given system through the acquisition of a few preliminary data sets, and then subsequently used to predict f r values for any given combination of S, N and ε within the same system-a potentially powerful capability both in research and industry. Of course, the formula can only be expected to hold for systems in the bouncing bed state-for beds which cannot be considered to move with a single flight time, i.e. higher-energy (large S, small N, high ε) systems in which particle motion is less coherent, clearly the resonances discussed in this paper cannot occur. Moreover, the resonant behaviour discussed within this manuscript is also observed to break down in the high dissipation (large N, low ε) limit, irrespective of driving strength, S. Again, this makes sense as in the inelastic limit, where all energy provided to the bed is instantaneously dissipated, one will, naturally, not see variations in the base-bed energy transfer characteristics. For the not inconsiderable ranges of parameter space for which resonance effects are observed, however, we find strong agreement between our simulated results and the predictions of (5) in all instances. In other words, our findings strongly suggest that for any system exhibiting the above-discussed resonances, the f r values corresponding to these resonances can be accurately determined using (5). This lends a great deal of credence to the validity of this empirical formula, thus supporting its applicability as a means of predicting the resonances of a given system-and hence determining the optimal manner in which to drive a given system in order to produce maximally efficient energy transfer.
Although the current framework is only directly applicable the somewhat restrictive case of a quasi-1D, vertical column of particles, it is hoped that with further research the findings presented here may be extended to more 'realistic' 3D systems. Nonetheless, even in its current form, the phenomenological formula presented in (5) may prove useful in certain applications, for instance in selective damping [66,67], where a simple column of particles-or the use of multiple individual columns-may be sufficient to produce the desired effects. Through this application alone, our results may prove highly useful in a variety of industries ranging from the manufacture of industrial machinery and power tools [68][69][70] to the construction of skyscrapers [71] and potentially even the design of future space vehicles [72,73].
Setting aside the specific, quantitative model presented , it should additionally be noted that despite the simple geometry investigated, the results of our study provide a valuable insight into the dynamics of dense, vibrated systems of any dimensionality-it is clear from our data that a similar resonant behaviour is found in 3D beds, albeit complicated by the presence of additional phenomena simultaneously influencing the bed dynamics.

Conclusions
Through the use of discrete particle simulations, verified by experimental data acquired using PEPT, we have investigated the manner in which the efficiency of energy transferral to a dense, vibrated granular bed may be influenced by the frequency of oscillations used to drive the system.
We show that even when the energy of the driving vibrations is held perfectly constant, the amount of energy successfully transferred from a vibrating base to a granular system can vary significantly dependent on the specific form of the driving used.
We demonstrate the existence of resonant behaviour in both 1D and 3D granular beds and, for the 1D case, provide a first systematic study detailing the dependence of the resonances observed on the system's key controllable variables.
Finally, we propose an empirical theoretical framework capable of predicting the resonant frequencies of a granular column based on known, controllable parameters, namely the size of the column, the elastic properties of the particles and the strength with which the system is excited.
This work also provides great scope for future research, for instance an extension of our investigation to the situation of higher-frequency driving, whereupon it may be possible to determine the point in parameter space for which a simple, linear relation between the bed flight time and the peak driving velocity of the system can no longer be assumed and a phase-dependent correction becomes necessary. We also hope that, with further research, a theoretical framework similar to that provided here may also be formulated for the case of 3D systems.
The findings of this study have a range of potential applications and implications, for instance in industry, where the maximization of energy transfer is often desirable, or geology, where a knowledge of resonant behaviour in particulate systems may prove highly valuable. Our results additionally demonstrate that the states and behaviours of a granular system cannot be sufficiently characterized by a single driving parameter, as is frequently assumed in past and current research. This work therefore highlights the necessity of a more rigorous approach to the systematic study of the dynamics of granular systems, and potentially a reassessment of the generality and applicability of existing studies and theories wherein a single energy-input control parameter is assumed.