Abstract
Granular flows require sustained input of energy for fluidization. A level of fluidization depends on the amount of heat flux provided to the flow. In general, the dissipation of the grains upon interaction balances the heat inputs and the resultant flow patterns can be described using hydrodynamic models. However, with the increase in packing fraction, the heat fluxes prediction of the cell increases. Here, a comparison is made for the proposed theoretical models against the MD simulations data. It is observed that the variation of packing fraction in the granular cell influences the heat flux at the base. For the elastic grain-base interaction, the predictions vary appreciably compared to MD simulations, suggesting the need to accurately model the velocity distribution of grains for averaging.
1. Introduction
Granular flows are of particular interest due to their vast and diverse applications as well as their intricate nature and dynamics. One of the major aspects of these flows is the dissipative nature of granular material. A granular material requires continual input of energy for a sustained flow. The amount of heat flux injected at the boundaries strongly influences the flow physics of granular materials especially in the case of vibrated beds. While simulating these flows, the behaviour of granular material at or due to the solid surface is an integral part of the solution for the entire flow field. Considerable efforts have been made to understand the influence of solid surface on granular flow physics, both experimentally and analytically. The experimental work of Savage and Sayed [1] shows the effect of rough walls. Profound effects of boundaries on the granular flow properties have been reported [2–5], which indicate that the drive surfaces are equally important as the granular material itself in determining the results of the test.
Boundary conditions are usually modelled from the first principles of momentum and energy balances at the surface [7]. Thus, by satisfying the expressions for the balance of momentum and energy, the granular boundary effect can be studied. Considerable work can be found on modelling the wall as rough particles which can be treated as granular materials principally, and at high speeds, the effects of shear can be neglected due to slip even leading to instantaneous collision [3–5]. In general, the energy supplied through base vibration E b is dissipated through particle-particle Dpp and particle-side wall collisions Dpw. Mathematically,
The input of energy E b of the granular cell is through the change in momentum ΔP a particle experiences in encounters with the base. The moving base interacts with the incoming/falling particle causing a change in particle energy ΔE as
where V(tδ) is the base velocity at time tδ. Through various methods of averaging, conditions have been obtained that apply to two-dimensional systems of identical disks or three-dimensional systems of identical spheres that interact with boundaries [5]. By assuming that the interaction is completely elastic and instantaneous, averaging over the time, total energy input for sinusoidal wave form can be calculated using [7, 8]:
where U is the velocity of the incoming particle, V is the RMS base velocity, N is the number of particles, g is acceleration due to gravity, and m is the mass of single particle. Depending on the unknown function
While for Maxwellian particle velocity distribution, [6] a slightly different relationship for symmetric base profile such as sinusoidal is proposed:
Both (4) and (5) coincide with (3) with different scale factors depending on the choice of
It can be concluded that the continuum level description of dissipative granular flows has led to the formation of granular flow models. However, with the increase of number density and/or dissipation, the models tend to fail. Thus, it is necessary to evaluate continuum order description for various flow regimes in order to establish their validity.
2. Vibrating Base-Heat Flux Boundary Condition
The analytical foundation of the heat flux boundary conditions proposed and used in modelling granular beds has largely been developed on the argument of instantaneous binary collisions. Conditions for a steady state of heat flux and momentum balance at a vibrating boundary have successfully been developed [3–6] and implemented in hydrodynamic models [10, 11]. One particularly effective model used in a number of simulations was developed by Richman [4] for bumpy walls and simplifies rather neatly the limiting case of a flat smooth vibrating wall. The results from the Richman method should be directly comparable with Warr's findings [3] in the limits of high frequency excitation. In this latter work, Warr and Huntley showed detailed calculations of the energy transfer from a vibrating boundary to a granular gas. From the mean change in velocity squared of the particles following a collision, they numerically evaluated the energy transfer integral, I e for a given base velocity magnitude. Keeping in view (3)–(5) and by incorporating momentum balance in the heat balance, I e can be related to the heat flux through [11] in nondimensional form:
where J b * is the nondimensional heat flux at the base, η o is the packing fraction, and T* is the non-dimensional granular temperature at the base (see Table 1 for details). G is the pair distribution function at the contact encountering the extended volume effect.
Relationship between the dimensional and dimensionless variables.
Calculations shown by Richman [4] quantified the heat and momentum flux for the rapid granular flows of identical smooth spheres that interact with bumpy boundaries through inelastic collisions. The boundaries were considered to be nonstationary as well and could translate randomly with given specific mean velocities along with deviations about the mean with specified fluctuation velocities. Richman produced expressions for mean energy transfer in all the three orthogonal directions moving with a mean square velocity of 3V2 (V2 for each of the direction), with V being the mean velocity of the base. The relationship includes addition and removal terms scaling with the square of base velocity and granular temperature. The bases of these calculations were upon the Maxwellian velocity distribution functions that relate and describe the velocities of both the flow particles and the boundaries. The rate of exchange was also calculated for the case of linear momentum and kinetic energy between the granular particle and wall particle for a nearly flat surface.
Later, Warr and Huntley showed detailed calculation for a triangular waveform inputting energy [3]. A detailed analysis considering all possible combinations of impacts with possible speeds and positions of a particle upon a vibrating base is provided by Warr and Huntley. The energy was calculated using the mean change in velocity squared of the particles due to collision, and it was evaluated numerically through an energy integral I e for a given base velocity magnitude. For these calculations, the requirement is on the vibration period which should be much shorter than the mean time between successive collisions made by a given particle and the base. This analysis showed how the energy integral varies with the speed of the incident particles and the velocity of the base while taking into account the coefficient of restitution between the particle and the base. Remarkably, it is shown that high incoming speeds of particles do not necessarily guarantee high energy inputs for small base velocities. In fact, it is also possible for incident particles to lose energy on the average during an impact. This allowed direct comparison to the work of Richman for the mean square velocity of base (3V2) and allowed a calculation of heat flux in a one-dimensional system with flat boundaries.
Kumaran also proposed energy input expressions for vibrated beds [6]. It gives an energy input proportional to the mean square speed of the bed. Although it is easy to implement, it does not allow for the possibility of energy extraction from the bed. The nondimensional form of I e * for the Kumaran boundary condition can be calculated using [6]:
While for Richman, boundary condition I e * is given by [3]
In general, the average energy input is proportional to square of nondimensional RMS base velocity V*2. As the base velocity increases, the predictions of the three heat flux models (Richman, Kumaran, and Warr) start to grow. Figure 1 shows the calculation of the energy integral for each of the formulations with e = 0.9. We see that at high temperatures (high base velocity), Warr's model is comparable with Richman's calculations, but that at low base temperatures (low base velocity), there are substantial discrepancies. One may notice that predictions using Kumaran are lower than Richman for the given granular temperature range. On the other hand, Warr's model approaches the predictions of Richman's model with the increase of base temperature. At low temperatures, the discrepancy between Warr's model and Richman's model is significant. The discrepancy between Kumaran's and Richman's approach can be related to the additional energy extraction/dissipation term at the base.

In vibrated beds, low granular temperature is often seen along with high packing fractions. The higher densities of granular material especially at higher levels of dissipations may influence the total heat flux at the vibrating boundary. The effect of dissipation during grain-base interaction is phenomenal on the amount of heat flux transfer. However, keeping the effect of grain-base dissipation apart, the effect of high densities is related to the pair distribution function, often modelled using [10, 11]:
The shielding effect due to higher packing fractions results in the likelihood of each grain to interact with base without the influence of its neighbouring. At a given granular temperature, one may observe the heat flux variation with increasing packing fraction in Figure 2 for elastic grain-base interaction. It can be observed that at lower packing fractions and lower temperatures (base velocities), the predictions of Kumaran's model match well with Richman's model. However, as the driving RMS base velocity is increased, the discrepancy between the two models starts to appear especially at higher packing fractions. At higher base velocities (higher temperatures), Kumaran's model predicts fairly high heat flux compared to the model of Richman. The observed discrepancy between the models may lead to significance differences in the flow physics predictions.

3. Steady State Hydrodynamic Model
For a vibrated granular cell, one-dimensional steady state variation in packing fraction or density and granular temperature can be described using
where ρ is the density, g is the gravity, P is the pressure, K is the thermal conductivity, and Dpp is the total dissipation of the granular media. For dilute limits K and Dpp, reduce to
The equation of state used for the case is given as
The equation set (10) requires three conditions to generate the solution. The first one is that the number of particles in the bed is fixed

Comparison of packing fraction (η) and granular temperature (T*) variation in 1-dimensional cell for N = 100 with MD simulations.
4. Simulation Geometry
The geometry used for the simulation is a vertically vibrated bed with elastic side walls and base. The vibrated bed is filled with smooth hard spheres under vertical vibration in the presence of gravity. The top surface of the cell is kept open. In absence of side wall dissipation, variations in hydrodynamic parameters of the granular flow occur in the vertical direction mainly. The coefficient of restitution between the grains is 0.9. The cell contains 100 grains of uniform diameter d = 2 mm under the influence of gravity g. The results for the variation of packing fraction and granulation temperature are compared with event-driven MD simulations.
4.1. MD Simulations
Hard sphere MD simulations are performed to validate the findings of hydrodynamic simulations. These simulations used event-driven scheme with hard sphere model. The ED scheme considers all the collisions as binary and instantaneous only. The base-particle coefficient of restitution was kept at unity to minimize any effects due to the inelasticity of the base. The initial configuration of the cell uses a random distribution of grains in the simulation cell, which is then run for a total of 109 number of collisions to avoid the possibility of biasness with initial placement in the cell. Evaluation of granular temperature and packing fraction and other related parameters is carried out through ensemble averaging over 109 collisions. The hydrodynamic scale trends are generated by mapping the location of individual grains on one-dimensional grid and storing their associated velocities.
For high base velocity, low density, and nearly elastic granular material, the density can be approximated as
By equating the energy flux at the bottom plate with the rate of dissipation of energy per unit area of the base, one may find that the base temperature is given by
It can be noticed that the base temperature is directly dependent on base velocity (temperature) and the packing fraction. With the increase of packing fraction (number of grains in cell), the heat flux variation in the granular cell also varies. Along with heat flux increase at the base, the dissipation also scales (
For dilute packing fraction of the cell, the model described using (10) with equation of state can be solved for the granular bed. Figure 3 shows the comparison of packing fraction and nondimensional granular temperature (T*). The comparison is carried out against the predictions of MD simulations. It can be observed from the variation of packing fraction (Figure 3) that peak packing fraction lies away from z* = 0 (location of base). Naturally, the highest temperature is observed at the base followed by a decreasing profile owing to dissipation. Compared to MD simulations, the comparison is reasonable with prediction of slightly higher peak packing fraction compared to MD simulations. Interestingly, the peak value of packing fraction lies near the lowest value of granular temperature (∼ z* = 10). The region 0 < z* > 10 has increasing packing fraction profile and decreasing granular temperature while moving away from the base. This pattern can be visualized as a block of granular material lying in the region of ∼ z* = 10, where dissipation rate is at its peak. Below the block of granular material, the grains have rapidly decreasing granular temperature as they move upward.
From the aspect of energy flux imparted through the base, such variation in the packing fraction profile can possibly lead to various types of grain-base collisions. Although the volumetric influence is catered using (9), still the chance for a successive double collision with base (sudden recollision) of a single particle may increase. In such a case, the grains located near the vibrating base experience push back from top material block, immediately after getting kicked from the vibrating surface. This could lead to an extended collision of the grain with the base. An extended collision is the one when the falling grain meets with the base also moving downward. In such a scenario, the resultant direction of motion of the grain will not change. In fact, the grain experiences a secondary collision when the base follows the upward part of the next cycle.
Figure 4 depicts the possible scenarios of grain-base interactions. Figure 4(a) represents the case of extended rapid recollision. The likelihood of this type of collision increases with higher packing fraction compared to case (b) and case (a) in Figure 4. At further higher packing fraction, the peak of the “maximum density” region increases with a shift in its location (see hydrodynamic predictions in Figure 5), only to increase the chances of sudden extended recollision. At the same time, one may notice that the packing fraction near the base also increases slightly.

The possible collision scenarios for grain-base interaction.

Moving packing fraction “hump” with the increase of granular cell packing fraction.
Interestingly, the two regions on either sides of the peak packing fraction are fairly similar, that is, the initial steep increase and the following decrease of packing fraction (see Figures 3 and 5). However, the difference in corresponding variation of granular temperature in both regions is phenomenal. The top layer lying above the “block” region of grains has nearly constant granular temperature, though MD simulations suggest slight increase in temperature. The increase of granular temperature away from the vibrating base has been previously observed [10]. While the bottom layer has a sharply decreasing granular temperature as the grains move away from the base, the nearly constant granular temperature indicates that the grains in top layer hardly change their energies; that is, there are less collisions between grains at that height, while the rapidly decreasing temperature suggests higher rates of dissipation in the region.
The variation of heat flux throughout the granular bed is also shown in Figure 6. As the packing fraction is increased, the heat flux at the base increases. At the same time, the increased dissipation reduces the fluidization height, resulting in shifting the “maximum density point” towards the base. In such a scenario, one may notice that the calculation of averaged heat flux in vibrated cell relies on the correct description of velocity distribution in the cell. In Figure 7, a comparison of heat flux predicted using MD simulations is drawn against the Kumaran and Richman models for elastic grain-base collision. One may notice that the predictions of the Kumaran as well as the Richman models match reasonably well with the results of MD simulations for low peak packing fractions. However, as the peak packing fraction increases, the deviation between the MD simulations and the Kumaran model becomes apparent. On the other hand, the predictions of heat flux using the Richman model are slightly lower even at low densities. With the increase in the packing fraction, results from MD simulation and the Richman model match somewhat better than Kumaran. However, the discrepancy is still present between the predictions of MD simulations and theoretical models. The results suggest that the theoretical models need improved predictions of velocity distributions of the grains near the base. That may provide better estimation of heat flux at different packing fraction conditions.

Variation of heat flux in the cell with the increase of granular cell packing fraction.

Variation and comparison of heat flux at the base with the increase of granular cell peak packing fraction at fixed base velocity of V b * = 1.74.
5. Conclusion
Granular flow in one-dimensional vibrated cell is modelled using steady state hydrodynamic model. The variation of packing fraction and granular temperature in the cell has dependence on the heat flux imparted to the dissipative media. In the limit of low inelasticity and packing fraction, estimation of heat flux at the vibrating base is predictable. However, with the increase in the granular cell packing fraction, a region of high density granular material is observed away from the base. The region acts like a top shield for the grains moving upward after receiving momentum from the base. With increase in dissipation in the region, the likelihood of grains escaping the block is very low. In such a scenario, the models for the velocity distribution at the base may no longer remain valid and shall be accommodated in the theoretical models. A comparison with MD simulations also shows appreciable deviation from the theoretical predictions. Experimental validation of the presented results can help in revealing the complex flow physics.
