Abstract
We present comparison of two closure models for dissipative flow of granular gas. Initial validation of the models is achieved using MD simulations for a vibrated cell with elastic side walls. With the dissipation at the side walls, convective rolling in the cell is observed. We show that dense inelastic granular gas model exhibits deviations from near elastic model even at low dissipation and dilute conditions. The flow physics shows that the strength of convective roll can be related with the differential dissipations on the side walls. The discrepancy between the two models is significant as we reevaluate the scope of near elastic model in the presence of dissipative boundary conditions.
1. Introduction
Dissipative gaseous flows are often countered in nature as well as in industries such as chemical, pharmaceutical, and power generation [1]. Understanding them from first principle is important not only due to their application but also due to their diverse nature and complex physics [2, 3]. The nature of these flows changes fundamentally as the inelasticity increases especially with the added effect of higher densities. A number of mathematical/hydrodynamic models have been proposed in last few decades to describe the physics of these flows [4–8]. However, validation of these models against experimental data is hardly found in literature apart from much widely used near elastic modelling approach (for instance, see [9–11]).
The investigation of the exact form of constitutive relationships for energy and momentum transport fascinates many, for instance [12], but at the same time poses a great deal of difficulty while simulating them. The use of different hydrodynamics descriptions is still an open issue due to the nonlinearity and complexities of the models [10, 13, 14]. There is still a bridge to be made between the theoretical framework and a real-time tool for use in practical industrial scale granular systems. Many industrial processes such as drying, cooling, agglomeration, granulation, and coating involve transport of mass, momentum, and energy of the granular media. Conventional methods used for fluids involve balances of these quantities, but for granular flows there are some fundamental differences. In order to understand granular flows a common method is to explore a simplified system that allows us to deduce the key physics, including, for example, shear cells, rotating drum and fluidized cells. A number of experimental methods such as NMR [15] and high speed photography [16], particle level molecular dynamics simulations [17] as well as bulk hydrodynamic models [18, 19] have explored systems that use these geometries to test and analyze granular physics.
Vibration is very helpful for fluidizing materials that have a high irregularity of shapes and sizes. It is especially helpful with those materials which fluidize poorly due to a broad particle size distribution. In simple dry vibrofluidized beds vibration is the only mode which fluidises the bed. In such a system, one observes that the behaviour of the fluidized material is strongly dependent on both the strength and type of shaking. Typically, the particles are set into motion via interaction with the boundary of the system, usually in a vertical position. When the boundary acceleration exceeds that due to gravity the particles are able to move ballistically, setting the system into a fluidized state. The motion of particles in such a state resembles the classical picture of a molecular gas and a number of successful models have built upon this analogy over the past three decades [4, 6]. As the principal difference between a classical and a granular fluid is the inelasticity of particle-particle and particle wall collisions [20, 21], a rational description requires considerations of factors such as the anisotropy, roughness of particles (which couples translational and rotational degrees of freedom), soft-particle interactions (which mean that particles deform during collisions) [18] alongside the inherent transient nature of shaken beds. In addition, the coefficient of restitution is a function of collisional velocity, and realistic particles tend to consist of varying sizes and shapes (poly disperse) [18].
It would be very difficult, however, to develop a hydrodynamic model with constitutive relations explicitly dependent on all of the above factors. With the exceptions of some attempts to approximately include friction and the use of a collisional model with a velocity-dependent restitution coefficient [18], most of these factors have not been incorporated to date. Besides this most of the work has been limited to dilute monodisperse systems utilizing constitutive relationships applicable to near elastic systems [9, 19, 22]. It is important, therefore, to test whether the available models capture the granular physics at moderate-to-dense flows with higher levels of inelasticity.
In this research work, we use monodisperse vibrated beds as a tool to investigate the effects of different types of constitutive models for granular physics. For comparison and further investigation molecular dynamics simulations are also carried out using the vibrofluidized dry granular bed. We compare results of two different models in the limit of near elasticity with dilute loading conditions. The models are solved using a finite-difference-based solver. At first the validation of dilute near elastic model for granular gas is carried out using MD simulations in a three-dimensional vibrated bed having depth of one particle diameter with elastic boundaries. For simplicity we assumed that the out of plane walls of the cell are frictionless and elastic for all simulations carried out in this paper. The effectively 2D bed contains spherical single-size grains with sides and base walls. Afterward the dense inelastic gas model is validated in comparison with the results for near elastic dissipative granular gas model for slightly dissipative gas with elastic side walls.
In the second stage we test the models with dissipative boundary conditions (side walls) and compare the granular temperature and packing fractions predicted using near elastic as well as dense inelastic model. We show that even at low inelasticity/dissipation and dilute state the near elastic model failed to capture the detailed physics predicted by dense inelastic gas models. Afterwards we show that dissipative side wall causes the formation of convective rolling in the cell. It is also shown that the strength of the roll and its location can be linked with the side wall dissipations.
2. Molecular Dynamics Simulations
Here an event-driven MD simulation is used which models particles as inelastic hard spheres where collisions are instantaneous and binary. To generate the initial configuration particles are inserted sequentially and randomly in the cylinder without overlap. A preliminary run was performed typically for 106 total number of collisions in order to allow the system to reach a stationary nonequilibrium state. To generate hydrodynamic scale trends using the molecular dynamics code, a virtual grid is generated to find the localized variations of granular temperature and packing fractions in the domain. The locations of the grains' centres in their local grid cells are noted, and associated properties are stored. The results are averaged over the predefined period of time and/or number of collisions, thus generating temperature and packing fraction readouts at all required instances at all places within the domain of cells. For steady state the averaging is carried out for the total duration of the simulation. To produce meaningful results, the simulation is run for 106 collisions per particle.
3. Mathematical Models
3.1. Hydrodynamic Model
The kinetic energy of the dissipative flows is generally achieved through balance of energy input with energy sink. The dissipative loss of energy in the bulk means that constitutive relationships are dependent on dissipation rate along with granular temperature and density. The steady-state balance of momentum equation describing a slightly compressible flow in a viscous fluid is given below along with the mass balance [6]:
where is the density, is the velocity field, is the pressure and , are the dynamic and dilatational viscosities respectively. Body forces are included in the momentum balance, , which can be due to, for example, gravity.
For the dissipative flows additional term corresponding to energy sink appears. The resultant energy equation in nonconservative form is given by [6]
where is the granular temperature [23], is the heat flux, is the density, and is a sink or source term.
3.2. Constitutive Relationships
A number of researchers have proposed transport coefficients to achieve closure of (1)–(3) [4, 7, 23]. For pressure , here it is calculated using the following equation of state [24]:
where is given by ,we use the approximation suggested by Carnahan and Starling [25] and Torquato [26] for as follows:
where is the packing fraction, and is the mass of a single particle of dissipative gas with particle diameter .
We use three closure models for comparison. The “near elastic model” uses constitutive relationships suitable for dissipative flows having coefficient of restitution very close to elastic limits. The relationships used here are proposed by Jenkins and Richman [6]. While the second form of constitutive relationship relies on the dense inelastic flow suitable for highly inelastic moderately dense flows. We use relationships as proposed by Garzó and Dufty [27]. The third model uses constitutive relationships developed using simplified approach for dilute dissipative flows. The constitutive relationships are proposed for dilute inelastic gas by Brey et al. and also include the packing fraction gradient term in heat conduction law [28]. The presence of additional term, additional constitutive relationship, has been argued in the past, but it has been shown as a part of the closure.
The expressions for transport coefficients of viscosity, thermal conductivity, and dissipation are only introduced here, and readers are referred to original research papers for complete expressions. The viscosity coefficients are given by
where and are elastic limits of viscosities [27]. The heat flux relationship has different forms for near elastic model and inelastic models. For near elastic model, it is related to thermal conductivity through conventional conduction law whereas for both dilute and dense inelastic gas models additional term appears in the relationship such that
where is the thermal conductivity, and is an additional coefficient that appears in inelastic model [27, 28]. The constitutive relationships for both are given as
where is the dilute limit of [27, 28]. The sink/source term is given by
where the first term on the right hand side is the energy dissipation rate due to particle-particle losses, and second term corresponds to work and viscous heating.
The relationships (6)–(10) are functions of packing fraction, granular temperature and coefficient of inelasticity. Here the lengthy expressions are omitted in the main manuscript for the sake of coherence and ease of readers. Readers are referred to the reference papers for complete derivation of transport coefficients.
3.3. System Description
We take a three-dimensional system with lengths of sides and in horizontal and vertical directions, respectively. The thickness of the cell is unit particle diameter, . The top of the cell in the vertical direction is kept open (see Figure 1). The dimensional parameters and their definitions are given in Table 1. The out-of-plane walls of the cell are assumed frictionless and elastic. In this effective 2D box, no mass can enter or leave the boundaries; however, energy can enter or leave. In order to maintain a gaseous state continual input of energy is required. The cell is initially filled with gas of uniform density under the influence of gravity () in the vertical direction. The cell contains spherical particles of equal diameter (diameter, = 2 mm), each having mass of -6 kg. The cell has energy intake from the lower boundary which is dissipated through bulk dissipation and through side walls. The forms of energy and momentum boundary conditions are described in the following section.
3.4. Boundary Conditions
The rectangular domain has four boundaries (Figure 1) which can be labelled as bottom, top, right, and left sides. In order to obtain complete description of the system thermal, and momentum boundary conditions are required at each of the four edges of the system.
3.4.1. Bottom Side
The heat flux at the vibrating boundary is related using the relationship developed by Jenkins and Louge [29] as follows:
where is the coefficient of restitution between particle and base. The momentum boundary condition is set to be a no-slip boundary as in the case of viscous fluids. The negative sign is removed to allow inflow of energy from the base.
3.4.2. Top Side
The boundary condition consisted of a zero heat flux with no-slip momentum boundary condition.
3.4.3. Left Side
The left-side wall is modelled as a dissipative no-slip boundary. The heat flux due to the inelastic particle side wall collisions is then given by
where is the particle side wall coefficient of restitution.
3.4.4. Right Side
The right-side wall is also modelled as a dissipative no-slip boundary. The heat flux due to the inelastic particle side wall collisions is then given by
where is the particle side wall coefficient of restitution.
Solution Methodology
An adaptive first-order time stepping with a second-order spatial derivative integration finite-difference scheme is used to integrate (1)–(3). Off-center difference scheme is used at the boundaries to ensure second-order accuracy at the boundaries as well. Adaptive time step is based on a Courant condition combined with maximum and minimum time step limits. Three different grid sizes (200 × 200, 250 × 250, and 270 × 270) were tried to obtain the grid independence. 250 × 250 point - grid was found to generate satisfactory accuracy. The boundary layers near the walls were also resolved using the same grid.
4. Results and Discussion
4.1. Validation of Near Elastic Model
The near elastic model of Jenkins and Richman [6] has been widely used and validated against MD simulations and experimental results for simplified geometries at low levels of dissipation. See for example, [19, 30, 31]. First we compare the results of the near elastic model of Jenkins with MD simulations. The comparison is carried out for the simple two-dimensional cell with vibrating base as shown in Figure 1. The rms base velocity is set at 3 in the cell which contains a total of 256 grains of 2 mm diameter each. The side walls of the (see Figure 1) are set as elastic walls while the top end is open. The coefficient of elasticity at the vibrating base is also set to unity whereas the coefficient of restitution of grain-grain interaction is fixed at 0.98. The comparison of packing fraction and granular temperature is shown in Figures 2(a) and 2(b), respectively. The comparison shows a good agreement between the MD simulations and hydrodynamic model. The results show that the near elastic model gives good predictions especially for the case of near elastic dilute cell loading.
(a)
(b)
One can observe a slight increase in granular temperature for MD simulations results in the region of > 20. The near elastic model does not show such increase in temperature. Despite the simplicity of the near elastic models it has been under extensive criticism especially in the presence of better mathematical models [27]. It has been argued that additional terms are not negligible even in simplified geometries.
4.2. Comparison of Dense Inelastic and Near Elastic Models
In the next step we compare the findings of near elastic model with that of dense inelastic model. Garzó and Dufty showed that dense inelastic model reduces to near elastic model for dilute and slightly inelastic case except for the influence of additional term [27]. Therefore, here for the sake of highlighting the difference between the two we take a case of further higher base velocity of for comparison. This causes the increase of granular bed height and reduction in peak packing fraction as observed in Figure 3. The validation in previous section at higher peak packing fraction and lower height of fluidization suggests that the near elastic model provides reasonably good predictions at higher base velocities. Therefore, a comparison with near elastic model is shown here while the results of MD simulations are not shown in Figures 3 and 4. The comparison between the packing fraction variations along y* is good between the two models as seen in Figure 3. Similarly the comparison between the granular temperatures predicted by both models show good agreement for . For , the dense inelastic model shows considerable deviation from the predictions of near elastic model. The increase in the granular temperature is due to inclusion of additional term as noted previously. The upturn is consistent with the earlier findings and is often referred as thermal inversion and near elastic models shows no such phenomenon as inclusion of such additional term is neglected. Evidences from MD simulations (e.g., see [17]), and NMR experimental results [32] suggest that the granular temperature upturn or thermal inversion is observed in dry vibrated granular beds.
4.3. Effect of Dissipative Side Walls
In the case of one-dimensional near elastic and dilute conditions both models show good agreement except away when due to the temperature upturn. The dissipation of side walls is now introduced to create gradients in the horizontal direction also. Different coefficients of restitution are used at the left and right hand sides of the two-dimensional section (Figure 1). is taken as 0.62 at the left side wall while is set as 0.79 for the right side wall. The resultant temperature variation is shown in Figures 5(a) and 5(c) for near elastic model and dense inelastic model. It can be noticed that the horizontal variations are influenced by the increase of temperature towards the top side of the cell. Near elastic model shows less horizontal variations compared to others. Also the spread of high-temperature zone near the base is wider and higher in the near elastic model.
(a) Near elastic
(b) Near elastic with correction
(c) Dense inelastic
Influence of side-wall dissipation is also seen in the normalized velocity vectors (Figure 6) showing a rolling pattern centred towards the left side of the cell. This is expected as previously rolling pattern has been linked with the strength of side-wall dissipation in vibrated beds [30]. Also the shifting of the location of roll towards the left side of the cell is expected due to increased dissipation on left side wall compared to right hand side wall. For both models no apparent change in rolling pattern is seen.
(a) Near elastic
(b) Near elastic with correction
(c) Dense inelastic
The fundamental difference between the near elastic and dense inelastic model is the influence of additional term in the heat flux relationship. Brey et al. proposed such an inclusion for dilute case [28]. As the granular gas used for simulations is dilute, we compare the results for the constitutive relationships proposed by Brey et al. suitable for dilute inelastic loading in Figures 5(b) and 6(b). It can be noticed that variation in granular temperature predicted by dense inelastic model is similar to Brey et al.. A comparison of granular temperature at the side walls is also presented in Figures 7(a) and 7(b) for right and left hand sides, respectively. It can be noticed that modifications introduced by Brey et al. are significant as it matches reasonably with the dense inelastic model at dilute loading conditions. On the other hand Jenkins’ near elastic model predicts significant deviations compared to both models. It can be concluded that influence of thermal inversion term is not only limited to unidirectional variation of granular temperature rather this term affects other spatial variations as well.
(a)
(b)
Despite extensive mathematical formulations and evidences from MD simulations and experimental results, use of near elastic models is still favoured due to their simplicity and ease of use in most of the commercial dissipative/granular gas software (see, e.g. Fluent etc.). At higher level of inelasticity and higher cell loadings the deviations are expected to grow larger and such a scenario shall limit the use of near elastic models. In the next section, we show the behaviour of both models under the influence of different side-wall dissipation combinations to understand the limitation of near elastic models.
4.4. Convective Rolling Strength and Side-Wall Dissipation
The influence of the side-wall dissipation on the bulk physics is considerable. The generation of convective roll in the domain and its strength has been linked with the side-wall dissipation [30]. It is anticipated the higher the dissipation at the side walls, the stronger will be the roll strength. The nondimensional angular velocity () [30] of the roll is shown in Figure 8 for the net differential coefficient of dissipation at the walls. The net difference is obtained by setting one wall as elastic and reducing the coefficient of restitution at the other side. The net differential coefficient of restitution () is varied from 0 to . It can be seen from Figure 8 that both models show no rolling with elastic walls. As the relative side wall dissipation is increased, we observe deviations between the near elastic model and dense inelastic model. This shows that influence of the boundary conditions results in change in flow physics even for the case of dilute loading and low grain-grain dissipation. It can be observed that the dense inelastic model shows reasonable match with the MD simulations results. While near elastic model shows lower strengths of roll. The result of Figure 8 suggests that the choice of constitutive relationships for simulating granular flow is not simply based on the level of dissipation between grains. The dissipative nature of boundary may influence the flow physics even in simplified geometry.
5. Conclusion
The description of granular flows using constitutive relationships of different mathematical foundation is an issue of debate. The bulk physics is reasonably well predicted by the hydrodynamic formulations used here. However, it is shown here that without the inclusion of additional coefficients (thermal inversion term especially), the complete description of the dissipative flow is not achieved even at low dissipation and dilute loading. The results show that the boundary conditions are also influential despite low intergranular dissipation even at dilute loading. The findings here are important as we evaluate the scope of the near elastic models commonly used in commercial codes. This research work does, however, set the basis for describing the granular flow using appropriate closure models, and in the future we anticipate that effects of further inelasticity and enhanced density could be studied.