Steady \& unsteady fluidised granular flows down slopes
D. E. Jessop, A. J. Hogg, M. A. Gilbertson, C. Schoof

TL;DR
This paper develops a model for fluidised granular flows down slopes, incorporating kinetic theory, and validates it with experiments on steady and unsteady flows, revealing key dynamics and simplifications.
Contribution
It introduces a kinetic theory-based framework for modeling fluidised granular flows down slopes, including asymptotic reductions and experimental validation.
Findings
Model accurately predicts steady flow velocity profiles.
Unsteady flow front positions match theoretical similarity solutions.
Sidewall resistance significantly affects unsteady flow evolution.
Abstract
Fluidisation is the process by which the weight of a bed of particles is supported by a gas flow passing through it from below. When fluidised materials flow down an incline, the dynamics of the motion differ from their non-fluidised counterparts because the granular agitation is no longer required to support the weight of the flowing layer. Instead, the weight is borne by the imposed gas flow and this leads to a greatly increased flow mobility. In this paper, a framework is developed to model this two phase motion by incorporating a kinetic theory description for the particulate stresses generated by the flow. In addition to calculating numerical solutions for fully developed flows, it is shown that for sufficiently thick flows there is often a local balance between the production and dissipation of the granular temperature. This phenomenon permits an asymptotic reduction of the full…
| kg/ms | |||||
|---|---|---|---|---|---|
| kg/m3 | |||||
| kg/m3 | |||||
| m | |||||
|
|||||
| cm3/s |
| /[∘] | /[cm3/s] | Experimental measurements | Estimate | ||
|---|---|---|---|---|---|
| /[cm] | /[cm] | /[cm] | /[cm] | ||
| 10 | 33.1 | 0.56 | 0.74 | 1.42 | 0.68 |
| 33.8 | 0.61 | 0.96 | 1.00 | 0.77 | |
| 53.5 | 0.76 | 1.35 | 1.35 | 0.93 | |
| 15 | 11.2 | 0.52 | 0.69 | 0.94 | 0.43 |
| 38.2 | 0.73 | 1.07 | 1.19 | 0.71 | |
| 42.5 | 0.77 | 1.09 | 1.25 | 0.74 | |
| /[∘] | /[cm3/s] | ||
|---|---|---|---|
| 10 | |||
| 15 | |||
| /[∘] | /[cm3/s] | /[cm/s] | /[1/s] |
|
||||
|---|---|---|---|---|---|---|---|---|
| 10 | 33.8 | 7.1 | 113 | 1.68 | 0.22 | 0.022 | ||
| 53.5 | 6.6 | 147 | 1.20 | 0.30 | 0.014 | |||
| 15 | 11.2 | 0 | 61 | - | - | - | ||
| 38.2 | 4.2 | 134 | 0.83 | 0.37 | 0.020 | |||
| 42.5 | 6.7 | 133 | 1.34 | 0.23 | 0.019 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Steady & unsteady fluidised granular flows down slopes
D. E. Jessop\aff1, 2\corresp
A. J. Hogg\aff3
M. A. Gilbertson\aff4
C. Schoof\aff2
\aff1Laboratoire Magmas et Volcans, Université Clermont-Auvergne-CNRS-IRD, OPGC, Clermont-Ferrand, France \aff2Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, Vancouver, Canada \aff3School of Mathematics, University of Bristol, University Walk, Bristol, BS8 1TW, United Kingdom \aff4Department of Mechanical Engineering, University of Bristol, University Walk, Bristol, BS8 1TR, United Kingdom,
Abstract
Fluidisation is the process by which the weight of a bed of particles is supported by a gas flow passing through it from below. When fluidised materials flow down an incline, the dynamics of the motion differ from their non-fluidised counterparts because the granular agitation is no longer required to support the weight of the flowing layer. Instead, the weight is borne by the imposed gas flow and this leads to a greatly increased flow mobility. In this paper, a framework is developed to model this two phase motion by incorporating a kinetic theory description for the particulate stresses generated by the flow. In addition to calculating numerical solutions for fully developed flows, it is shown that for sufficiently thick flows there is often a local balance between the production and dissipation of the granular temperature. This phenomenon permits an asymptotic reduction of the full governing equations and the identification of a simple state in which the volume fraction of the flow is uniform. The results of the model are compared with new experimental measurements of the internal velocity profiles of steady granular flows down slopes. The distance covered with time by unsteady granular flows down slopes and along horizontal surfaces and their shapes are also measured and compared with theoretical predictions developed for flows that are thin relative to their streamwise extent. For the horizontal flows, it was found that resistance from the sidewalls was required in addition to basal resistance to capture accurately the unsteady evolution of the front position and the depth of the current and for situations in which side-wall drag dominates, similarity solutions are found for the experimentally-measured motion.
1 Introduction
Particles are often transported in the form of dense currents under the influence of gravity. Their bulk flow rate is greatly enhanced if part or all of their weight is supported by a gas flow through them. When particles are poured onto a slope that is less than their angle of repose, they are held stationary by the action of contact friction and merely flow down the surface of the pile in a thin layer as more grains are successively added. When they are poured onto a slope that is steeper than the angle of repose, a thin, dense current forms in which the particles move in bulk down the slope (see, for example, Ishida et al., 1980; GDR MiDi, 2004). If a gas is passed vertically through the particles then the drag it exerts on the particles bears some of their net weight and hence the frictional forces decrease. Consequently, the effective angle of repose of the particles decreases as does the minimum slope angle at which bulk flow flow takes place (Nott & Jackson, 1992). When the gas flow is sufficiently large for the entire weight of the particles to be supported (i.e. the particles are fluidised), then bulk frictional forces are insignificant and very mobile currents form, even on horizontal surfaces. The influence of a fluidising gas flow through particles on their mobility is exploited widely in industrial settings where it is necessary to transport bulk materials either to move them from one place to another using air slides (which can be several kilometres long), or to keep horizontal surfaces clear of particles in pieces of processing equipment such as circulating fluidised beds (Savage & Oger, 2013). There are also features in many particulate environmental flows in which there is significant upward gas flow and this enhances their speed and range (e.g. Druitt, 1998; Roche et al., 2004).
There have been extensive studies of the flow of particles down a slope and of some of the effects of fluidisation. A common approach to mathematically modelling these motions is based on a continuum description that couples expressions of mass conservation with expressions of the balance of momentum within each phase (e.g. Nott & Jackson, 1992). Under this approach, the fluidised material is treated as two inter-penetrating phases that interact with each other. The models do not resolve the motion of individual particles, but rather the evolution of averaged, bulk properties, which depend upon the net effect of direct interactions between particles within the current and between the particles and their surroundings. The duration of contacts between the constituent particles has important consequences: if the contacts are sustained then they are likely to be frictional in nature; if they are instantaneous then they are collisional in nature (e.g. Campbell, 2006). The stresses induced by instantaneous collisions between pairs of particles (i.e. in dilute and rapid granular flows) can be evaluated through the use of granular kinetic theories (Jenkins & Savage, 1983), in which a key dependent variable is the granular temperature, , a measure of the variance of the instantaneous velocity field. Hydrodynamic equations of motion have then been derived for granular materials that are much like those for dense gases except there is substantial energy dissipation through inelastic collisions (see Jenkins & Savage, 1983; Lun et al., 1984; Haff, 1983, for example). It is possible, of course, for there to be collisions over a range of durations and these ideas do not translate to dense and slowly-shearing flows where contacts are prolonged and thus, in part, frictional. The action of the interstitial fluid is a further factor that needs to be considered when modelling granular flows. For example, in contrast to the original studies of granular kinetic theories, Koch & Sangani (1999) proposed that interaction with the fluid could generate agitation within the flows and that fluctuating viscous forces could be the generators of particle temperature.
There have been relatively fewer studies that report granular flows that are aerated or fluidised. An early approach was to treat the fluidised particles as a non-Newtonian fluid of power-law rheology, sometimes with a yield stress (Botterill & Bessant, 1973; Botterill & Abdul-Halim, 1979; Ishida et al., 1980; Savage & Oger, 2013). This approach can be made to work well in specific practical situations (Singh et al., 1978), but is entirely reliant on empirical methods to determine the effective rheology in each circumstance since such approaches do not capture the fundamental dynamics of the particle motion.
A more fundamental approach is to model the evolution of averaged properties of the inter-penetrating phases. Ogawa et al. (1980) modelled steady, one-dimensional fully-fluidised currents down slopes. They derived constitutive relations based on the collisions between a particle and its neighbours, which were represented by an imaginary spherical shell surrounding it. This resulted in a balance between collisional stresses and gravitational forces. Nott & Jackson (1992) coupled a kinetic theory for collisional grain flows with a Coulomb-like model for frictional effects to predict the bulk mass flow rate of aerated grains down an inclined channel. The experiments and model featured gas flow rates up to the minimum required to fully fluidise the particles. Their mathematical model of friction in the flows followed Johnson & Jackson (1987) and Johnson et al. (1990) and assumed that the frictional component (dominant at high particle volume fractions, ) was simply added to the collisional component (dominant at low ). They pursued a similar approach to the interaction term between the gas and the particles adding together a contribution based on the Ergun equation (dominant at high ) and one from the Richardson and Zaki equation (dominant at low ). No contribution was included from slip between the two phases in the direction of the slope. Oger & Savage (2013) took a similar approach (although with some different closures of the models), again retaining a frictional term, and solved the resulting equations using the MFIX numerical code to study the dynamics of granular motion within air slides, computing the steady, fully-developed velocity and granular temperature fields for flows within a channel of rectangular cross-section. Finally, Eames & Gilbertson (2000) reported the unsteady flow of fluidised materials along horizontal surfaces. For their system, they showed that collisional stresses would be small compared with those associated with fluid drag and so when fully fluidised, the force balance set hydrostatic pressure gradient against fluid drag terms. We will show below how our work differs from their modelling framework and yet is able to reproduce features of their experimental results.
Key to furthering our understanding of the dynamics of fluidised flows is direct and detailed experimental evidence against which theoretical models can be validated. However, there are few measurements of fluidised granular currents, especially down slopes. Previous experimental studies have presented bulk properties such as total flow depth and mass flow rate (e.g. Nott & Jackson, 1992; Eames & Gilbertson, 2000). Some measurements of local properties such as velocity have been made though this has often been with instruments such as optical probes or turbine elements (e.g. Botterill & Bessant, 1973, 1976; Ishida et al., 1980; Nott & Jackson, 1992). Whilst providing important information, the disadvantages to these techniques are that they lack spatial resolution, are intrusive (especially in fluidised particles Rowe & Masson, 1981 and offer only point measurements i.e. traverses are necessary to build velocity profiles and they are therefore only suited to steady flows. More recently, Particle Image Velocimetry (PIV) has been applied to fluidised systems such as static beds (Bokkers et al., 2004) and dam-break experiments over horizontal surfaces of initially fluidised, fine natural volcanic ash (Girolami et al., 2010). PIV has the advantage of offering high spatial resolution, and allows instantaneous velocity fields to be calculated. The experiments of Girolami et al. (2010) had a short-lived phase of quasi-constant flow following the initial release of material; however the grains were not continuously fluidised along the apparatus. This meant that even though the materials were highly expanded initially because of the very small particle size, they decelerated rapidly due to the loss of mobility associated with compaction in the terminal flow phase. As such, they are not representative of fully-fluidised flows.
The aim of the present work is to understand better the dynamics of fluidised granular flows by providing further experimental evidence and proposing a new unsteady model of these flows that fully takes into account the interaction between the particles and the fluid and incorporates collisional stresses. Both of these processes play a crucial role in the dynamics of fluidised granular flows in which the gas flow bears most of the weight of the particulate layer and the particle interaction contribute significantly to the shear stresses developed by the flow. This implies that the dynamics are different from ‘dry’ granular flows in which the role of the interstitial fluid is negligible (e.g. Lun et al., 1984; Forterre & Pouliquen, 2008; Woodhouse et al., 2010).
In this work, experimental measurements were made of granular currents over a range of slope inclinations and conditions and the experimental arrangement is described in § 2. The measurements were made in an apparatus that confined the flow between two walls, which enabled the overall size and shape of currents to be measured over time. In addition, PIV was used to measure the velocity profiles of the particles within the currents, enabling their overall behaviour to be linked to their rheology. § 3 develops the general continuum model and the equations of motion for the flowing state. This builds upon the ‘two-fluid’ approach in which the gas and grains are treated as two inter-penetrating phases (Jackson, 2000). Fully developed flows are tackled in §4 and compared with experimental observations. The continuum model in this section is analysed in the regime for which the properties of the flowing layer vary only with distance from the underlying boundary and the solutions are computed numerically and asymptotically in a regime where the flow thickness far exceeds the diameter of an individual grain. Unsteady and transient effects found in flows along inclined channels are investigated in §5 and a new model developed in the ‘lubrication’ regime where the downslope length-scale is much large than that perpendicular to the slope. Flows along horizontal surfaces differ their counterparts along inclines (§6) and measurements of their inherently unsteady motion are reproduced well by a new self-similar solution to the flow model in the lubrication regime. Finally our findings are summarised and discussed in §7. We also include two appendices. In the first we analyse the consequences of an extended kinetic theory, following the constitutive laws of Jenkins (2007). In the second the effects of the side walls are analysed in the regime that the flow depth is much less than the channel width.
2 Experimental approach
2.1 Experimental setup
The experimental arrangement conformed to that shown in Figure 1. The apparatus was a long, narrow channel (1 cm 100 cm, 50 cm in height) which could be inclined to some angle, , to the horizontal. The bottom of the channel was a porous plastic distributor material (Vyon ‘D’) through which dry air was passed from a windbox below at a speed , but for which the pressure drop over it was much larger than that through the granular flow. This ensured that the gas flow was evenly distributed i.e. the gas flow entering the apparatus was uniform and perpendicular to the distributor plate so that at the base of the granular layer, the gas velocity is where is the particle volume fraction evaluated at the base of the flow. The particles were constrained between vertical parallel walls, so that the motion is effectively two-dimensional and the motion of particles within the current could be seen. The front wall was was made from a glass sheet allowing the flows to be viewed and the other sides were made of aluminium plate. The rear plate was painted black to increase the contrast between it and the white particles. Particles entered the apparatus at one end (the uppermost when inclined) through a funnel giving a constant volumetric flux, , which could be changed between experiments by changing the aperture of the funnel (Nedderman, 1992). The flow rate was the flow rate of the current based on the bulk volume when the particles were at rest; so, the volume flow rate per unit width of particles where is the particle volume fraction of a static bed of particles, and the distance between the front and back of the flows is cm. It could be controlled by using funnels of different sizes, each of which could then be associated with a bulk flow rate, ; however, this is a nominal flow rate as the actual flow rate could vary from occasion to occasion. The apparatus had closed ends; so, to avoid ‘backing-up’ when running experiments with a non-horizontal slope, particles were removed from the downslope end using a vacuum cleaner. This had no measurable effect on the height profiles obtained but allowed experiments to be run for longer. No removal of particles was necessary for the slower-moving horizontal flows. The value of was accurate to within .
The material used for all the experiments was approximately spherical, glass beads (Potters Ballotini) with particle diameters in the range 250–425 m and a mean diameter m. We measured the particle volume fraction of densely packed, static material (i.e. a maximum) as , which is close to the maximum value of for random, close-packed, mono-sized spheres (Jaeger & Nagel, 1992), and the (unfluidised) bulk density was g/cm3. The powder corresponds to a class B powder according to Geldart (1973), so no bubble-free expansion when fluidised is expected within a ‘static’ bed. The minimum fluidisation speed, , was found by independent experiment where the gas flow rate through a static bed of material was gradually increased and the resulting pressure drop through the bed measured (Davidson & Harrison, 1963). We found that the entire weight of the bed was supported when cm/s.
2.2 Shape of currents and front position extraction
The flows, viewed from the side, were recorded using a digital video camera. Calibration was performed using an image of a block of known dimensions placed in the apparatus once the camera was set up in position for a given experiment. Still images from the recorded experiments were analysed by transforming the RGB images to grey scale. These were then turned into binary images through thresholding. Though the threshold value was calculated automatically, the contrast between black back wall and white particles meant that the resulting binary images were robust and consistent. The upper and lower surfaces of the outline of a current were defined as the first and last white pixels when descending a column of the binary image. 95% confidence intervals for height measurements are cm. The front position was taken as the point where the top surface met the bottom one.
2.3 Velocity measurements
PIV was used to make measurements of the velocity fields of the flows using a high-speed video camera capturing at 500 frames per second, close up to a particular region of the flow. The PIV measurements required the flow to be seeded with marker particles for which we used the same-sized particles as for our other experiments but approximately one third of which were dyed black. The properties of the dyed particles (, angle of repose etc) were identical to the non-dyed particles.
Two-dimensional velocity fields were calculated by processing image pairs (two consecutive frames) from the video taken by the high-speed camera using the open-source Matlab-based DPIVSoft2010 code. The software makes an initial estimation of the velocity field on a coarse grid and then uses this to translate and deform the interrogation window in the second image in keeping with the deformation of the flow field. Errors associated with image pattern distortion, as is the case when velocity gradients are large, are greatly decreased using this method (Meunier & Leweke, 2003). Several initial iterations were run to get a good approximation for the flow field. A final run was performed with an interrogation window of 3232 pixels () and velocity vectors were calculated using a 50% overlap between adjacent windows. A median filter was then applied with a limit of 0.5 to remove spurious vectors (e.g. Adrian & Westerweel, 2010, pp. 406).
Instantaneous velocity profiles may not be representative of the flow as a whole. In particular, the bubbles of gas that could form spontaneously in the flows often disrupted the instantaneous velocity profiles. However, flows down the steeper slopes in our experiments, 10∘ and 15∘, reached a steady state very quickly, and for these flows an ensemble average of the flow velocity could be found by averaging over both many points in time and at several positions along the flow. The quality, and hence the accuracy, of time-averaged velocity fields has been shown to be greatly improved when the average instantaneous correlation function is used to calculate the velocity field (Meinhart et al., 2000). We therefore modified the PIV routines accordingly to produce a single time-averaged velocity field per experiment using an interval of twenty frames (0.04 s) between image pairs, and fifteen image pairs per experiment. This interval is larger than characteristic time for shear ( s), so the velocity fields at successive intervals are uncorrelated. Velocity profiles were then formed from the stream-wise vectors of the time-averaged velocity field lying on a depth-wise transect at points separated at intervals of 1 cm and averaged to form the ensemble average velocity profile. The resolution of PIV measurements can be expressed as (Adrian, 1991),
[TABLE]
where is the uncertainty of locating the centroid of the correlation peaks, is the magnification factor of the lens, and s is the time step between images. For our setup, , and so that cm/s.
For the steady flows it is more useful to define error based on the sum of variances of the all the profiles used to calculate the ensemble averaged standard deviation over the images given by
[TABLE]
This average standard deviation was then used to calculate 95% confidence intervals for the velocities.
3 Equations of motion
We investigate the motion of granular currents down an inclined surface when the particles are fluidised, as shown in figure 1. These flows are gravitationally-driven, but do not accelerate unboundedly; instead the principle action of particle interactions is to contribute to the shear stresses that balance the down slope acceleration and potentially lead to steady motion. We formulate a mathematical model of the two-phase motion that couples mass conservation for each phase with expressions for the balances of momentum and we show how this formulation may be applied to steady fully-developed flows that vary only with distance from the underlying boundary (§4), and to unsteady, relatively thin flows for which the acceleration perpendicular to the underlying boundary is negligible (§§5,6).
The mathematical model is built upon a continuum description of two inter-penetrating phases which interact with each other. These models do not resolve the motion of individual particles; rather, they allow the computation of the evolution of averaged properties. Such approaches have been employed often for confined, horizontal fluidised beds (e.g. Bokkers et al., 2004; Goldschmidt et al., 2004), but these studies differ from the dynamics of the flows analysed in this contribution where there is persistent shear flow down the inclined surface. The flows analysed here also differ in an essential way from non-fluidised granular motion down inclines since the support of the weight of the grains by the imposed gas flow significantly reduces resistive forces and increases mobility. Nevertheless, we find that steady flows are admissible and thus the motion must develop sufficient shear stresses to balance gravitational acceleration. Our model assumes that these stresses arise from particle interactions and are collisional and the particle fluctuations may be characterised by a granular temperature since friction as a bulk property is virtually eliminated by fluidisation and the viscous forces associated with interstitial gas flow are negligible. The granular temperature will be shown to be relatively small and thus the interactions generate only relatively weak shear stresses, but these are sufficient to balance the gravitational acceleration.
The collisional nature of the motion is justifiable in all but some small regions of the currents, for example close to the surface of the slope. The model captures only the relatively slow evolution of averaged quantities. In particular bubbles (i.e. volumes largely evacuated of particles that travel through fluidised particles) are not explicitly resolved. Bubbles are an important feature of deep, static fluidised beds as apart from strongly affecting the local instantaneous volume fraction of particles and they are the primary source of granular temperature in such a bed (Menon & Durian, 1997). There are several processes that might lead to the generation or suppression of bubbling, most notably including the dissipation of granular temperature through collisions, which is prone to clustering instabilities (Goldhirsch & Zanetti, 1993; Fullmer & Hrenya, 2017). Some studies have sought to predict the onset of bubbling in static beds through linear stability analysis (e.g. see the review by Jackson, 2000). The flows of fluidised materials analysed here are somewhat different from these stability analyses, however, due to the persistent production of granular temperature by work done by the velocity field shear against the shear stresses, a process absent in static beds; hence, by means of a scaling analysis Eames & Gilbertson (2000) showed that the contribution of these bubbles to the overall balance for granular temperature is likely to be negligible for this downslope motion. Furthermore, shallowness in the bed is thought to suppress bubbling (Botterill et al., 1972; Tsimring et al., 1999), as is shear (Botterill & Abdul-Halim, 1979; Ishida et al., 1980). We therefore assume that bubbling is likely to have a limited influence on the fluidised currents. Extensive bubbling was not observed in the currents. The photograph shown later in figure 9 is typical with no apparent bubbles. While agitation was visible at the top of the currents, bubbles sufficiently large to fill the width of the bed were hardly ever seen.
Most of the theoretical developments in this study will be for two-dimensional flows and the effects of the front and back walls of the apparatus are neglected. The use of this planar set-up allows the structure of the system to be seen and measured (as described in §2), but at the expense of it being bounded by walls not present in realistic, three-dimensional systems. Arguably, because fluidisation eliminates internal friction, a large part of the effect that the presence of these walls might also be eliminated. Here, in most of what follows, we analyse the motion in the regimes that the side-walls play a negligible role; however in §6.1, we also analyse the case when the side-walls have a dominant effect on horizontal currents and in Appendix B, we derive the extra, weak retardation on flows down slopes that arises from side wall drag when the depth of the flow is much smaller than the width.
The general equations of motion for a continuum model, known as a ‘two-fluid model’, of a gas-particle system have been developed by Jackson (2000). The conservation of mass in each phase is given by
[TABLE]
where denotes the volume fraction of solids and and the velocity field of the gas and solid phase, respectively.
Following Jackson (2000), the balance of momentum for the gas is given by
[TABLE]
and for the particles
[TABLE]
where and are the densities of the gas and solid phase respectively, is the spatially-averaged stress tensor of each phase with the superscript (, ) denoting the gas or solid phase, respectively, is the drag force exerted by the particles on the fluid due to the difference in their velocities and denotes gravitational acceleration. The material derivatives, and denote the rate of change moving with the gas and the solid phase respectively. A number of researchers, including Ergun (1952) and Jackson (2000), have suggested that , where is a drag coefficient. Virtual mass and particle shear forces are neglected.
These equations will be solved for the situation shown schematically in figure 1. The slope is inclined at angle, , to the horizontal with the underlying boundary at and the upper surface of the current at , while the -axis is aligned with the basal boundary. A mixture of solid particles and gas runs down the slope under the influence of gravity.
4 Fully developed flows
4.1 Model for fully-developed flows
First, fully developed flows are investigated, in which the dependent variables are functions only of the distance from the boundary, , and the velocity fields of the gas and solids are given by and , respectively. Conservation of mass for the solid phase is automatically satisfied by this form, but for the fluid phase we deduce that
[TABLE]
where is the fluidising gas flux per unit area normal to the boundary.
The expressions for the balance of momentum follow those proposed by Johnson & Jackson (1987) and Agrawal et al. (2001) where for the gas phase down the slope
[TABLE]
where denotes gravitational acceleration and is the gas viscosity. Perpendicular to the slope, we find
[TABLE]
where is the pressure within the fluid phase. In (8) and (9) we have assumed that the gas phase is incompressible and can be treated as Newtonian with constant viscosity, .
For the solid phase, the balance of down slope momentum is given by
[TABLE]
while normal to the slope,
[TABLE]
In (10) and (11), and are components of the solid phase stress tensor, , and at this stage we have not yet invoked any constitutive model for these stresses. Further, from (10) the driving force for the current is gravity and within this framework, currents over horizontal surfaces are inherently unsteady as they decelerate.
The downslope balance of momentum (10) differs from previous contributions. Nott & Jackson (1992) implicitly assumed that there was no relative component of velocity downslope between each of the phases and thus there was no drag force (i.e. ). Eames & Gilbertson (2000) did not consider momentum balance for the fluid phase and imposed ; thus within their model, the drag force is dominant and by assumption the shear stress associated with the solid phase is negligible. We do not invoke either of these assumptions at this stage, instead maintaining the various dynamical processes until their relative magnitudes have been fully assessed below.
Adding the normal momentum equations (9) and (11), we find that
[TABLE]
When the current is homogeneous so that particle volume fraction is constant, from (7) the vertical component of the gas velocity is also constant; so, (12) expresses the hydrostatic balance between the vertical gradient of the normal stress from both solid and fluid phases and the weight of the fluidised grains.
It is also insightful to eliminate the fluid pressure field between (9) and (11) to find that
[TABLE]
This expression reveals the fundamental dynamical role played by fluidisation. The slope-normal component of the inter-phase drag, incorporated into the model by , can balance the weight of the grains and thus it is possible for the normal stress tensor of the solid phase to be much reduced from its non-fluidised magnitude. This in turn reduces the magnitude of the solids shear stress and thus the mobility of the fluidised flows is greatly enhanced. Equation (13) is different from the classical model of a static fluidised bed because velocity gradients lead to normal stresses in the solid phases and these may contribute in a non-negligible way to the balance between weight and drag as shown below.
4.1.1 Inter-phase drag and constitutive equations
The drag on the solid phase due to the fluidising gas flow is given by , where the drag coefficient, , may be written
[TABLE]
where and are given in Table 1 (Ergun, 1952). The first term on the right-hand side of the equation represents the drag associated with viscous processes, and the second term with inertial processes. For the regime of interest in this study, the inertial effects are negligible since the Reynolds number, based on gas velocity and particle size, is sufficiently small ; however for completeness at this stage we maintain it in the model formulation. Other expressions for the drag coefficient have been used (e.g Agrawal et al., 2001; Oger & Savage, 2013) and these could replace (14) within this modelling framework.
It was noted above that particle interactions are dynamically important because of the momentum transfer arising from particle collisions (Lun et al., 1984; Garzo & Dufty, 1999). Here we examine the collisional stresses and follow Nott & Jackson (1992) and Agrawal et al. (2001) amongst others who incorporate these effects into models of fluidised and aerated flows, to write the shear and normal components of stress in terms of a granular temperature , which measures the fluctuations of velocity about the mean, the volume fraction of solids , and the coefficient of restitution , which characterises dissipation in the instantaneous collisions. While the constitutive laws invoked here have been validated in some scenarios by simulation and experimentation, there remains some uncertainty about their generality. Hence we pose the model quite generally so that the constitutive relations could be updated as required.
For fully developed flows, we write
[TABLE]
where and are dimensionless functions given in Table 1. In this study, we employ the constitutive formulae derived by Garzo & Dufty (1999) and recently used for modelling dense avalanches by Jenkins & Berzi (2010).
Granular temperature may be generated and ‘conducted’ via the flow processes and dissipated in the collisions. Following Lun et al. (1984) and Garzo & Dufty (1999) amongst others, these effects are encompassed in the following expression of energy balance within the flow
[TABLE]
where the flux of granular temperature is given by
[TABLE]
and , and are dimensionless functions given also in Table 1. In posing this balance of granular temperature we have neglected generation and dissipation of the granular temperature mediated by viscous interactions with gas (Koch & Sangani, 1999).
Following the formulation of Agrawal et al. (2001), dissipation by viscous processes is much smaller than dissipation through inelastic collisions when
[TABLE]
Furthermore, the generation of granular temperature by viscous processes is much smaller than that by granular interactions when
[TABLE]
The constitutive laws, , as well as those involved the boundary conditions (, see §4.1.3) feature the radial basis function, . Various authors have suggested forms for and we employ an expression that is close to the suggestion of Vescovi et al. (2014), who empirically fitted a function to match data from discrete element simulations. Importantly, the radial basis function diverges as the volume fraction approaches maximum packing (as established by Torquato (1995)) and following Vescovi et al. (2014) we write
[TABLE]
where the weighting function is given by
[TABLE]
Thus, when the radial basis function is given by the formula proposed by Carnahan & Starling (1969), but it exceeds this value when and diverges as maximum packing is approached. Vescovi et al. (2014) suggest that and that . While this choice ensures that and its derivative are continuous at , the second derivative is discontinuous. This is problematic for the system of differential equations that we will integrate numerically; therefore, we employ the values and , which ensure that is sufficiently smooth. Moreover, our expression (22) with these values is close to those proposed by Vescovi et al. (2014) and Torquato (1995) and appears to match the simulation data adequately.
The energetic balance encompassed in (17) assumes that the particles are sufficiently agitated so that the particle diameter is the appropriate correlation length scale over which dissipated occurs (and the rate of dissipation is then given by ). Recently, however, Jenkins (2007) has suggested that at relatively high concentrations, clusters of particles begin to form and thus the correlation length increases to () and then the rate of dissipation is given by . This extended kinetic theory has been applied to unfluidised flows of grains down inclined planes by Jenkins & Berzi (2010, 2012), where an empirical formula for , informed by comparison with simulations and experimental measurements, is proposed in terms of the dependent flow variables. In our study there is potentially the need to include this phenomenon into the modelling framework to obtain good comparison between the predicted and measured results. However, as shown in appendix A, we find that extended kinetic theory makes negligible difference to the model predictions for the fluidised flows in our regime of interest and so we do not include it in the calculations that follow.
4.1.2 Coefficient of restitution
The dynamical effects of collisions between particles and between particles and the underlying boundary are characterised in the model by three parameters: the coefficients of restitution between the particles, , and between the particles and the boundary, , and the specularity coefficient , which governs the dynamic interaction between the particles and the bottom surface (36). These parameters are relatively difficult to measure directly.
The coefficient of restitution, , plays an important role in continuum models and in Discrete Particle (or Element) Models (DPM), which endeavour to calculate the motion of large ensembles of particles and to resolve individual particle collisions. In continuum models, the difference of from unity is proportional to the rate at which the collisions dissipate energy (see the definition of in table 1), whereas in DPMs it controls the ratio of normal velocities before and after binary collisions and in these models, there are potentially additional means of energy dissipation. Often values for the coefficient of restitution are adopted without independent experimental confirmation and for DPM studies, typical values are relatively high (for example, and respectively in the studies of Goldschmidt et al. (2004) and van der Hoef et al. (2008)). These values are close to measured values of discrete collisions (see, for example, Kharaz et al., 2001). When used in kinetic theory models, commonly adopted values of are rather lower and Jenkins & Zhang (2002) suggest a means by which the the appropriate value for kinetic theories can be derived from directly measured normal and tangential coefficients of restitution and the tangential coefficient of friction. For glass spheres of 3mm diameter, the measured data of Foerster et al. (1994) corresponds to an effective coefficient of if the method of Jenkins & Zhang (2002) is employed and this is the value we employ in this study. We have no direct measurements of the appropriate coefficient of restitution for the collisions between the particles and underlying boundary; we choose , but note that its magnitude has very little influence upon the computed flow profiles apart from within thin basal boundary layers.
4.1.3 Boundary conditions
The boundary conditions for this problem follow the formulation of Johnson & Jackson (1987) and Johnson et al. (1990). At the base, there is no-slip for the fluid phase, a slip condition for the particle phase, and a condition specifying the flux of granular temperature. These are respectively given by
[TABLE]
What this means physically is that solid-phase slip is allowed and stress is transmitted in the down-slope direction by specularity i.e. the degree to which the angle of exit of a particle after collision with the base is different from the entry angle. This is mathematically represented by the specularity coefficient . Furthermore, fluctuation energy (granular temperature) is generated at the bottom surface and potentially dissipated by inelastic collisions, encompassed through a coefficient of restitution, .
At the top of the current , there are the free-surface boundary conditions that fluid shear and normal stresses vanish, given by
[TABLE]
However in addition, the flux of granular temperature vanishes and the solid phase normal and shear stresses adopt small values, representing the surface as being the location where collisional behaviour ends and instead the particles follow ballistic trajectories (see Johnson et al., 1990). Thus we enforce
[TABLE]
The boundary conditions (23)-(25) are of the same character as those employed by researchers in other flow regimes (see, for example, Jenkins & Berzi, 2010) and as for the constitutive laws, the framework for analysing these flows is robust to variations in the closures used for these conditions.
4.1.4 Non-dimensionalisation of equations
We now identify typical dimensional scales for the dependent variables and assess the magnitude of the various terms in the governing equations. It is convenient to sum the down-slope momentum equations of each phase (8) and (10) to eliminate the inter-phase drag so that
[TABLE]
In this expression, the key driving force is the down-slope gravitational acceleration and it is this term that the other terms must balance. Since the density of the gas is much smaller than that of the solid phase and the effects of gas viscosity are negligible away from boundaries in this streamwise balance, we deduce that the dominant resistance is provided by the shear stress associated with the solid phase. Coarsely scaling the variables and assuming the volume fraction and the constitutive functions of it are of order unity,
[TABLE]
Furthermore, if the granular temperature is in local equilibrium between production and dissipation (an assumption that will be tested in the numerical solutions that follow), then from (17),
[TABLE]
whence, the scaling for the velocity field is given by
[TABLE]
It is now convenient to introduce dimensionless variables, given by
[TABLE]
This set of scalings for the dependent variables of fluidised flows differs from those for flows of unfluidised, collisional granular media (Woodhouse et al., 2010). For non-fluidised flows, the granular agitation must provide sufficient normal stress to support the weight of the overlying layer. This would require the granular temperature to be of magnitude , which is considerably larger than the estimate deduced here (28) unless the motion is along relatively steep inclines . For fluidised flows, however, granular temperature is generated by collisions but the imposed gas flow through the underlying particles provides most of the normal stress to balance the weight of the flowing layer. The granular temperature, therefore, is lower and consequently the shear stresses are lower, which in turn significantly increases the mobility of these flows. Hence these fluidised flows are characterised by relatively high flow speeds and relatively weak resistance.
The model is characterised by five dimensionless groups:
[TABLE]
which represent respectively the inclination of the underlying boundary, the relative density of the gas to the solid phases, the size of the particles relative to the flow depth, the magnitude of the drag exerted by the fluidising gas flow relative to the weight of the granular layer and the reduced Stokes number, which compares particle inertia to fluid viscous effects.
It is also possible to define a particle Reynolds number,
[TABLE]
Notably, the Reynolds number defined in this way is independent of the inclination of the slope. The magnitude of the various model parameters for the experiments are set out in Table 2.
These scales may be used to show that the granular temperature dissipation and generation by viscous processes are negligible compared with direct particle interactions: (19) is satisfied when and (20) when .
We have five governing equations: mass conservation (7), down-slope fluid momentum conservation (8), the combined normal momentum equation from which the fluid pressure has been eliminated (13), the down-slope solids momentum equation (10), and the equation for the conservation of granular temperature (17). Non-dimensionalised, these become
[TABLE]
where measures the magnitude of the dimensionless relative velocity between the phases and is given by . From (23), the dimensionless boundary conditions at the base are given by
[TABLE]
while from (24), we enforce at the top surface
[TABLE]
The system of governing differential equations (31)-(35) and boundary conditions (36)-(37) form a seventh order differential boundary value problem. We use (31) to eliminate in favour of and we also evaluate in (35) by explicitly differentiating (33). The system is then integrated numerically. (For this task we employ the boundary value problem solver bvp4c in MatLab.) Example solutions for the volume fraction, , the granular temperature, and the velocity of the solid phase, are plotted in figures 2–4 for various values of the governing dimensionless parameters. We do not plot the gas velocity, , because outside of thin basal boundary layers, it is indistinguishable from the solids velocity, . This basal boundary layer exists because the while the gas phase satisfies a no-slip condition, the solid phase exhibits slip.
The general trends are that the volume fraction is approximately uniform while the granular temperature decreases with distance from the bottom boundary. Additionally there is a small slip velocity for the solid phase and the velocity shear decreases with distance from the the boundary. There are some systematic deviations from these general trends, most notably in regions close the upper and lower boundaries. Since the interactions with the basal boundary are more dissipative than interactions with the constituent grains, the granular temperature decreases within a region in the vicinity of the boundary. This boundary effect is diminished as the flow becomes thicker (i.e. as decreases, see figure 2), but is magnified where either the fluidising gas flow is increased (figure 3) and the slope is increased (figure 4). We also find that throughout the bulk of the domain, away from thin layers adjacent to the upper and lower boundaries, the production of granular temperature is in close balance with its dissipation, thus confirming the dimensional scales identified above (28).
It is of particular interest to evaluate the dimensionless flux of solids per unit width and the average concentration of particles, respectively given by
[TABLE]
and these are plotted as functions of the dimensionless parameters in figures 5, 6 and 7.
From figure 5, we observe that the dimensionless volume flux per unit width, , and the average volume fraction, , do not vary strongly with the relative particle size, . Thus, we deduce that boundary-related effects on the bulk characteristics are negligible for flows that are in excess of twenty particles thick. This result is of particular significance when unsteady shallow flows are analysed (§5).
The effects of the fluidisation velocity, , are rather more subtle (see figure 6). Increasing increases the normal support of the weight of the granular layer due to the gas flow and this lowers the concentration of the layer. However, the net volume flux, , does not vary monotonically with . Indeed, for the parameters in figure 6, it attains a maximum at a dimensionless gas flow rate and at that value of , . This reflects the trade-off between the increased mobility but lower solids fraction of more dilute currents. Finally there is also relatively complex behaviour with increasing slope angle (figure 7). For the computations in this figure, as we increase the inclination, we also adjust and , but maintain constant (see (29) and (30)). We find that as the slope increases, the normal stress developed by the particle collisions increases with increasing granular temperature and this supplements the fluidising gas flow, leading to a progressively decreasing average volume fraction. However the dimensionless volume flux exhibits a more complicated dependency because while the velocity fields increase with increasing slope, the volume fraction diminishes and eventually becomes sufficiently dilute for to be maximised at some finite value of . (For the parameters analysed in figure 7, the local maximum in the flux occurs at .)
4.1.5 Asymptotic solution
In the bulk of the flow away from the boundaries, it is possible to deduce an asymptotic solution to the governing equation for the regime and . This regime will have a widespread validity as in order to use a continuum approach, and for gas-solid flows . From (32), we note that to leading order and away from boundaries, the downslope velocities of the two phases must be equal (). Furthermore, from (35) there is a local balance between granular temperature production and dissipation such that
[TABLE]
provided the volume fraction of particles is not too small (i.e , so that the ‘conductive’ effects of the granular temperature remain negligible).
The governing equations for the normal and perpendicular momentum balances, (33) and (34), are then given by
[TABLE]
These reduced governing equations neglect shear stresses in the fluid phase, which become non-negligible as the basal boundary is approached and which allow the velocities of the two phases to differ. Also the ‘conduction’ of granular temperature is neglected because, to leading order, we find a balance between production and dissipation (39). When , the lower boundary layer corresponds to a region within which the velocity of the solid phase is small, while the upper boundary layer to a region within which the granular temperature is small. The leading order boundary conditions are then given by and the volume fraction at the base is given by , which is determined by substituting for from (39) into (36),
[TABLE]
where the constitutive functions are evaluated at . The basal volume fraction is thus a function of , , and .
Rearranging (40) and (41) and denoting , we find that
[TABLE]
where denotes differentiation with respect to . It is straightforward to integrate numerically these coupled first-order equations subject to the boundary conditions and . The solutions are plotted in figures 2, 3 and 4 and it is evident that these asymptotic solutions accurately reproduce the numerical solution of the complete system (very often in these figures, the asymptotic curves are indistinguishable from the numerical solution of the complete system).
There is also an even simpler approximate solution. The coupled system admits a homogeneous solution when
[TABLE]
where the constitutive functions are evaluated at . In this case, the temperature gradient is constant, , with . This solution is ‘attracting’ in the sense that trajectories in the phase plane approach it when
[TABLE]
which in turn demands that and that (see figure 8c). If these inequalities are not held then the reduced system evolves towards a state different from a uniform volume fraction and may not admit solutions at all. Physically, when , the dissipation of granular temperature, here encapsulated through a coefficient of restitution , is insufficient to allow for a steady balance between the weight of the flowing layer, the fluidising drag and the normal stresses generated through particle interactions (a balance expressed by (40)). When , we find that this limitation does not play for parameter values associated with the flows considered in this study and that a flow with a homogeneous volume fraction of particles provides a good representation of the more complete dynamics (see figures 2-4).
When the dimensionless fluidisation velocity and the slope are set, the average volume fraction within the current, , can be calculated using (45). Figure 8a shows the effect of on the curves of as a function of slope when the fluidisation flow is constant (). The curves in this plot are continued up to the maximum value of the slope, for which the reduced model leads to a homogeneous volume fraction and it can be seen that the slope at which this can be achieved is successively reduced as dissipation in the collisions is decreased. From figure 8b, for a given slope, , and fluidisation gas flow rate, , flows with lower coefficients of restitution lead to higher dimensionless volume fluxes per unit width. This is simply rationalised: as a high coefficient of restitution implies reduced dissipation and high granular temperatures. Consequentially there are higher stresses and greater resistance to the downslope motion.
Since the granular temperature must vanish at the surface to leading order, we find for the simple approximate solution with uniform volume fraction that the granular temperature is given by
[TABLE]
and the velocity field of the solid phase is given by
[TABLE]
where . The scaled slip velocity at the wall can be added to (48), but when , it is negligible. The velocity profile (48) is similar to the dimensionless ‘Bagnold’ velocity profile up to the factor , which controls the mobility of the flowing layer and is influenced by the fluidisation velocity. In many situations the approximate solution provides a very good representation of the solution to the complete system (see figures 2 to 4).
Also in this regime , the approximate solution yields
[TABLE]
and this is plotted in figures 5, 6 and 7, once again illustrating the utility of this asymptotic solution. The quantity thus plays a crucial role in determining the dimensionless flux, and in figure 8d, we plot its dependence on volume fraction for a range of values of . We note that is maximised for (with the precise value weakly dependent on ) and vanishes both when vanishes and when it approaches maximum packing. This variation reflects the balance between fast moving dilute flows and slow moving concentrated flows, leading to a flux maximum at intermediate values . Finally, we note that a dimensional estimate of the depth of a fluidised current may be obtained
[TABLE]
where is the dimensional flux per unit width at the source and the effects of slip at the wall has been neglected.
4.2 Experimental measurements of fully developed flows
4.2.1 Depth of currents
A typical velocity profile is shown in figure 9, superimposed on a captured image from the recording of an experiment. There is a small slip at the lower boundary, an approximately linear increase in velocity with distance from the wall until a maximum velocity is attained and then a progressive drop to zero. There appears to be a top to the current where the particle volume fraction suddenly drops and there . is greater than the height at which the maximum velocity is attained. Above particles are detected, but their velocity drops with increasing height and it has a large variance. This is consistent with there being a ballistic region into which individual particles may be projected. The height at which particle velocity drops to zero is the top of the entire current and is denoted by . The depth of the current can fluctuate a little with time (see figures 18 and 19, top). As a result, the averaging process will occasionally include points that are above the average height of the current so that the averaged velocity at these points will be necessarily lower than in the bulk of the current and the variance will be higher. The measured heights for the different currents are summarised in Table 3 compared with the height estimated from (50).
The asymptotic solution yielded an approximate formula linking height and source volume flux (50), and this is shown in figure 10. In this figure, fixed values of were used because from (45), depends on and it was not possible to find a solution for the full experimental range of for a fixed value of . The decreasing effect of on mobility as its value approaches reflects the effect it has on mobility shown in figure 8d. The theoretical formula contains no adjusted parameters and is an approximation to the more complete description, but it yields a reasonable quantitative representation of the relationship between the depth of the flowing layer, the source flux and channel inclination.
4.2.2 Velocity profiles
The ensemble averaged velocity profiles for 10∘ and 15∘ slopes measured halfway along the tank are shown in figure 11. The features described for figure 9 are reflected in each of the profiles. The structure of the currents posited by the model is similar in structure to the experimental measurements up to the point at which the velocity is a maximum (see figure 2). Above this point, the variance of the velocity measurements increases markedly as the velocity drops-off. As described above this effect could be because the PIV is simply sampling ballistic particle trajectories.
The velocity profiles measured on inclinations of , with source fluxes and are distinct from each other with the former forming a current that is deeper and much faster than the latter. It is not clear why there is such a large difference between the measured profiles. One possibility could be that the flowing layer exhibits multiple states for the same imposed flux. This behaviour is known in models of unfluidised granular flows (Woodhouse et al., 2010); however we failed to find such multiplicity of solutions in the governing equations examined in this study in the parameter regime corresponding to these experiments. The relatively fast and expanded flow with leads to small estimates of particle volume fraction with an excessive portion of the flow where (see below, section 4.2.3. Based on , ; based on , ), which seem physically unlikely, and so this experimental run is not reported further.
Figure 12 shows some instantaneous velocity profiles halfway along the tank for shallower slope angles whose motion may not be steady. The velocity profiles for the 3∘ and 5∘ slopes are similar in character to the averaged profiles for steeper slopes.
4.2.3 Particle volume fraction
The expectation from the model is that the particle volume fraction would be approximately uniform within the fluidised currents. It is not possible to analyse the degree to which is a function of position in the currents from our experimental set-up; however, it is evident that towards the top of the current above , drops sharply so that the current loses its opacity. This is consistent with the decreasing particle velocity there. Some bubbles were seen in the currents, but not many and at a small number of sites, even once the current had traversed the bottom of the container, and those that were seen were small in size.
From the measured velocity profiles, it is possible to estimate the average volume fraction, , by integrating the particle velocity profiles and dividing the measured particle flow rate by the result. The results are shown in Table 4, using as the overall depth of the current, to give measured volume fraction of the currents (). The measured values of particle volume fraction can be compared with estimated values, , which have been calculated using (45).
There can be quite good agreement between and despite several inherent uncertainties in their computation. Overall, the values of were comparable to the values of and with with the typical values of for static fluidised beds (Epstein & Young, 1962). In between and particles are present, but in practice the fall-off of velocity above is sufficiently rapid that this makes very little difference to calculations of : if is calculated on the basis of the top of the current being rather than , then its value decreases by less than 0.02, except when and , when it reduces by 0.09.
4.2.4 Slip at the wall
Specularity coefficients have not been measured for fluidised granular currents and some modellers think they should not be used at all for individual collisions (Goldschmidt et al., 2004). Their value is sufficiently badly defined that in their investigations of bubbling fluidised beds of glass particles Altantzis et al. (2015) used values between and and Li et al. (2010) from [math] to . Our computations showed that apart from within relatively narrow layers close to the boundary, the magnitude of the specularity coefficient had relatively little effect upon the flow profiles. It is, however, possible to estimate the value of from the directly measured slip velocities and gradients using the boundary condition (23) and the definition of in Table 1, so that in terms of dimensional variables
[TABLE]
The results are shown in Table 5, and it can be seen that a measured average value of is 0.28. The values for shown in Table 5 should be treated as only indicative as the velocity gradients close to the wall are shallow and the slip velocities small, so small variations in the velocity profiles can result in significant changes in the value of ; however, despite these uncertainties, the value of is reasonably consistent. The effect of on the velocity profiles is to introduce a slip velocity proportional to . Even for the relatively large values of estimated here, the magnitude of this dimensionless term is relatively small.
4.2.5 Scaling of the velocity profiles
The measured velocity profiles scaled as and are plotted in figure 13. With the exception of the lowest flow rate when , the data collapses well in the region close to the wall. The model predicts dependence of the theoretical curve on the slope angle for the flow through its influence on and hence , but this is not reflected in the experimental curves for which the scaling appears to eliminate the effect of . is also affected by the value of . Increasing causes a decrease in the predicted dimensionless velocity: the theoretical curves shift towards the left and the difference between the curves for and becomes less (see inset, figure 13).
5 Unsteady, developing flows on slopes
The mathematical model may be extended to unsteady, developing flows of fluidised currents down slopes, but it now takes a somewhat different form because stream-wise gradients can no longer be neglected. In this situation, we analyse the motion in the ‘lubrication’ regime, for which a representative streamwise length scale, , far exceeds a representative length scale perpendicular to the boundary, (). This means that accelerations perpendicular to the boundary are negligible and that to leading order, the normal stresses adopt the ‘hydrostatic’ balance given by (12). We again assume that the flows are many particles thick, , the density of the gas is negligible relative to that of the solids, , and the effects of gas viscosity are negligible (see §4.1.5). The leading order, dimensional momentum equations of each phase parallel with the incline then take the form,
[TABLE]
where the average volume fraction in the flowing layer, determined by the balance between the fluidising gas flow and the particle weight, is given by (45). It is interesting to note from (52) that now there must be a leading order difference between the downslope velocities of the two phases. Furthermore, since the flow is spatially and temporally evolving, we must include the inertia of the solid phase, which in (53) is given by the term (here denotes the material derivative). Summing these two momentum balances to eliminate the inter-phase drag and assuming further that the stresses in the solid phase are isotropic and that the current is in hydrostatic balance (12), we deduce that
[TABLE]
The granular temperature of the flow in this regime is assumed to be in local balance between its production and dissipation through collisions, and these processes dominate its advective and diffusive transport. We proceed further by adopting the appropriate dimensionless scales following the distinguished scaling identified in §4.1.4 and embodied in the dimensionless variables of (28) and (28). However, here we non-dimensionalise the depth of the flowing the current by a representative depth-scale , , and additionally
[TABLE]
Then, using the approximate form of the solution established in §4.1.5, we find that the depth-integrated, dimensionless momentum equation is given by
[TABLE]
where and
[TABLE]
In this setting, as for Kumaran (2014), measures the relative magnitude of the inertial to resistive terms. Unlike a Reynolds number for viscous fluid flows, it features only the length scales in the problem and , because both the inertial terms and the shear stresses are proportional to the square of velocity.
To proceed further we assume that the velocity field adopts similar dependence to (48) on distance from the boundary and this permits the evaluation of the integral and boundary quantities in terms of the average velocity and the depth of the layer:
[TABLE]
To complete the model, we express conservation of mass,
[TABLE]
This system is subject to the boundary condition that we impose a sustained source of particles at the origin
[TABLE]
Additionally, if the flow is supercritical then we must enforce the Froude number at the source. We impose the initial condition, and the current forms a front , such that .
Before constructing solutions, it is convenient to relate the height and streamwise length scales. We choose and thus . The lubrication regime requires that streamwise lengths far exceed the thickness of the flow; since the current is expanding in the streamwise direction, this regime is inevitably entered after sufficient time. However the adopted scaling may imply that in terms of these variables, the initial evolution may not be well captured by the lubrication assumption. We further choose the dimensional height , using (50) so that
[TABLE]
The governing equations now entail the single dimensionless parameter,
We construct travelling wave solutions for the dimensionless height and velocity fields. We write and , where is the dimensionless wave speed which is to be determined. Conservation of mass then implies that and in particular, the front speed is given by
[TABLE]
Balance of momentum leads to
[TABLE]
where a prime denotes differentiation with respect to . Distant from the front, the flowing layer carries a constant volume flux of material determined by the source conditions ( as ). Here we note that the travelling wave solutions do not satisfy the source condition precisely at , but instead it is satisfied as . For these flows, we find that the current adjusts over a short distance behind the front to a uniform depth and velocity and thus the travelling wave solution provides an accurate representation of the solution for the flow. Thus, we deduce that the position of the front and the far-field depth are given by
[TABLE]
Experimental measurements of the distance travelled by fluidised currents with time are shown in Figure 14, in which the inclination of the channel, the source volume flux and the fluidising gas velocity were varied; the measured flow speeds ranged over a factor of five. The measurements for the fully fluidised currents () after scaling are shown in Figure 15. From (62), should be proportional to and this is true even when the slope angles are small. Individually, the currents display a constant speed; however, the measured speeds can be significantly different from that expected from the model ().
The degree of data collapse for different parameters - , the nominal flow rate , and - is shown in Figure 16. For all three parameters, the scaling eliminates much of the scatter, but it is not fully eliminated. The collapse of data onto different lines with the same values of and for is excellent. It can be quite good for , especially at low . For the collapse of data is often incomplete. There will be some variation reflecting the difference in value of the true value of from the nominal value . It can be seen in Figure 16 that the effect of is small, but significant; however, it is eliminated after scaling, as shown in Figure 16b.
A systematic omission from our model is the effect of side wall drag and this could provide an additional resistance to motion, thus slowing the speed of propagation. In appendix B, we analyse the effects of the side walls when the height of the current, , is much less than the breadth of the channel, . We demonstrate that there is a weak retardation to the dimensionless speed proportional to when . We analysed the speed of the flow from the data plotted in figure 14 and found no systematic dependence on and thus there is no evidence that these relatively shallow currents were significantly slowed by side wall effects.
5.1 Time taken to establish steady uniform behaviour
The dimensionless profile of a fluidised current moving down a slope is determined from (63) and is implicitly given by
[TABLE]
where . We plot in Figure 17 the height of the travelling wave of material as a function of distance from the front for various values of the inertial parameter, , and note that the length scale over which the flow adjusts to the uniform depth, , increases with increasing . One measure of the streamwise length, , over which the flow attains its uniform depth is given by evaluating when , which when is given by
[TABLE]
At a fixed location, it is then possible to evaluate the dimensionless timescale over which the uniform depth is established, .
Figure 18 shows examples of the development of the fluidised currents at different angles of channel inclination. In figure 18, the flows on the steeper slopes relatively rapidly attain a uniform state in which the current does not vary along the apparatus, whereas those on shallower slopes and with smaller sources fluxes take much longer to approach this state. Flows along horizontal channels never approach a uniform state; instead currents adopt the shape of a wedge and do not progress at constant speed. These are are analysed in §6.
Figure 19 shows the change of height with time for five flows with the same nominal flow rate halfway along the apparatus before and after scaling. It can be seen that the scaled times at which the currents achieve a constant height (and systematically with ), as would be expected, but they are an order of magnitude larger than . For the expected values of , the currents would have to achieve their constant height very quickly, almost instantly, and for the values of corresponding to the experimental flows, from figure 17 the front of the currents would be expected to be ‘blunt-nosed’, with quite steep gradients of height at the front of the current. In fact, the front of the currents (figure 18) had a relatively shallow gradient.
6 Horizontal flows
Figure 20 shows the distance travelled by the front of a fluidised flow over horizontal surface as a function of time. It is evident that the currents do not travel at a constant speed. The shapes of the currents are shown in figure 21 and, ignoring the disturbance at the start of the currents at the point that they are poured into the system, they have an approximately triangular shape, though one with a low aspect ratio (i.e. their extent far exceeds their depth). Furthermore they flow through ‘bulk’ motion, not through the build-up of lamina arising from the constant avalanching down the current’s top surface seen for non-fluidised granular flows. They must also be scaled differently because the length and time scales introduced in (55) become singular when . To this end, we introduce the characteristic height scale, and define the following dimensionless variables
[TABLE]
The depth-integrated expression of momentum balance is then given by
[TABLE]
instead of (56), where the residual dimensionless parameter is the relative magnitude of inertial to resistive forces.
Conservation of mass is given by
[TABLE]
The appropriate dimensional depth-scale, , is determined from the source flux,
[TABLE]
so that the boundary condition is given by
[TABLE]
Flows over horizontal surfaces decelerate as the basal drag is no longer balanced by a sustained downslope acceleration. Thus, at sufficiently early times the flow speeds and depths are set by source conditions, and after the flow has propagated for sufficient time, the resistive forces become non-negligible and the motion enters a dynamical regime in which the drag force balance the streamwise gradients of the hydrostatic pressure and the inertial forces are negligible. Analogously to Hogg & Woods (2001), simple scaling shows that this regime is fully attained when . In this scenario, we deduce from (68) that
[TABLE]
and consequentially from (69)
[TABLE]
subject to the source condition
[TABLE]
(73) may be integrated numerically to reveal the evolution of the front position as a function of time and the variation of the depth of the current along its length; however, for these currents flowing over a horizontal surface, we may also construct a quasi-analytical similarity solution for the motion.
First, we determine the gearing between spatial and temporal scales that underpins the similarity solution for unsteady flow over a horizontal surface. To do this we scale and balance terms in the governing equation (73) and boundary condition (74). This yields
[TABLE]
Thus we deduce that and . We may then seek a similarity solution of the form
[TABLE]
where is a dimensionless constant to be determined as part of the solution and . On substitution in the governing equation (73), this gives
[TABLE]
where a prime denotes differentiation with respect to . This ordinary differential equation (78) is to be integrated subject to the boundary conditions
[TABLE]
The location is a singular point of the differential equation (78); we therefore start the numerical integration at , noting that
[TABLE]
It is then straightforward to integrate the differential equation (78) numerically and evaluate , , and so . The dimensional expression for distance covered by a horizontal, fluidised current with time is then
[TABLE]
We plot in Figure 22 the similarity solution for the height profile along the current noting that, again, the model predicts a blunt-nosed current that advances along the channel.
The scaled distance against time is shown in figure 23, and the data is collapsed sufficiently for the power-law form of the curve to appear to be reasonable, though the value of the exponent is different from that predicted. However, the shape of the current predicted by the scaled model is very different from the experimental measurements, taking the form of nearly flat current with a snub nose (figure 24).
6.1 Flow within a narrow channel
One of the differences between the experimentally-realised flows over horizontal surfaces and those down slopes is that the former are significantly thicker than the latter, and so it is possible for the side walls to have a strong influence on their development. Here, we analyse the motion of a fluidised current as it flows within a narrow channel of width , between sidewalls for which the streamwise extent of the flow far exceeds the depth of the current , which in turn far exceeds the width of the flow . For this regime, it is possible to simplify the governing equations for the unsteady motion down an incline on the basis that gradients across the flow are much greater than those in any other direction. In this scenario, the dynamical balance is somewhat different from that analysed in §4 and the resulting governing equations for the unsteady evolution of the thickness of the flow are also different (§5).
Our derivation of the governing equation in a narrow channel is developed from the dimensional expressions presented in §3 and then depth-integrated to establish a shallow layer model; however, it will be shown that it leads to a similarity solution with a different gearing between the spatial and temporal variables. Here we only present the governing equations for flows along horizontal channels, but the inclusion of a channel gradient is a straightforward generalisation and could lead to travelling wave solutions analogous to §5.
It is assumed that the solid particles are fully fluidised by an imposed gas flow and attain a state in which the volume fraction is uniform, . Since the flow is relatively thin, vertical accelerations are negligible and the motion is governed by hydrostatic balance given by
[TABLE]
In this expression we have neglected the contribution due to the weight of the gas phase since . In the downslope direction, after neglecting terms proportional to the density and viscosity of the gas, the combined momentum equation of both phases to leading order is given by
[TABLE]
where is the distance across the channel. Here it has been assumed that the normal stresses of the solid phase are isotropic and that gradients across the flow dominate all others. To complete this model, we introduce the granular temperature, which provided the channel is much wider than the grain size, is in local equilibrium between its production and dissipation. Then we may write
[TABLE]
and the constitutive law for the shear stress is given by
[TABLE]
Finally, by eliminating the fluid pressure from the normal force balances of each phase (see (13)), we find that
[TABLE]
The volume fraction, temperature, and therefore the normal stress component, , are independent of to leading order, and thus we deduce that
[TABLE]
This expression determines the average volume fraction as a function of the fluidising gas flux and marks a important departure from the shallow layer model of §5, because to leading order the solid stresses do not contribute to the support of the granular layer.
We progress by assuming that the velocity field of the solid phase exhibits cross-stream dependence, which is identical to that found in fully developed flows with the shear in the vertical plane. Thus we write
[TABLE]
and consequentially at . The dimensional governing equations then express conservation of mass and after sufficient time has passed so that inertia is negligible, a balance between the streamwise pressure gradient and the side wall stresses is given by
[TABLE]
For these flows over a horizontal surface we adopt the dimensionless variables given by
[TABLE]
where the height scale scale, , is determined by
[TABLE]
The dimensionless governing equation then becomes
[TABLE]
subject to
[TABLE]
We may construct a similarity solutions to this governing equation (93) and boundary condition (94) by first deducing the gearing between the spatial and temporal scales. The governing equation demands , while the boundary condition leads to . Thus we deduce that , and that the similarity solutions of the following form may be sought
[TABLE]
where and is a constant to be determined. In dimensional variables,
[TABLE]
On substitution of into the governing equation (93), we deduce that
[TABLE]
subject to and . The similarity differential equation (97) is singular at and the numerical solution may be initiated from a series solution valid close to that value, given by
[TABLE]
when . It is then straightforward to integrate (97) to compute the height profile, shown in figure 25, and the dimensionless constant .
The scaled distance against time for the experiments is shown in figure 26. There is reasonably good collapse of the data (better than in figure 23). Furthermore, it appears to follow a power law, though the exponent of the power law is slightly small than that predicted. The measured, scaled shape of the currents is plotted in figure 27. The predicted shape is now much closer to that of the experiments, although the scaling does not collapse completely all the measured data.
7 Discussion and conclusions
This investigation of fluidised granular currents reveals important distinctions in their dynamical properties from both dry granular flows and static fluidised beds. Most significantly, there is substantial and sustained shear in the velocity profiles. Consequentially particles are driven into each other and this provides a mechanism for the generation of stresses. In the regime we investigated, the inertia of individual grains remains relatively high and thus the particles interact with each other through dissipative collisions, and it these interactions that lead to the shear stress that balances the downslope gravitational acceleration. Granular flows in the absence of fluidisation must generate sufficient normal stresses to support the weight of the flowing layer and thus typically their shear stresses are also relatively large; however, fluidisation changes the balance of forces acting on the current. The fluidising gas flow provides most of the normal support to the flowing layer and thus both normal and shear stresses from the particulate phase are reduced relative to their non-fluidised counterparts, leading to flows that are much more mobile.
In this study we have formed a framework of modelling fluidised currents based on the solid-phase stresses that is generated from collisions between the particles. The degree of agitation in the system is measured through the granular temperature and constitutive laws are employed to determine the stress tensor in terms of the gradients of the velocity field, the granular temperature, and the volume fraction of solids, as well as several material parameters. For the regime studied here in which the flows are many particles thick, a local balance emerges between the generation and dissipation of granular temperature. This leads to an accurate asymptotic model for the complete dynamics, in which the flowing material is essentially modelled by a non-linear, local rheology. Furthermore, this reduction leads to a Bagnold-like expression between the flow depth and the flux of particles carried by the current, with the volume fraction, determined by the fluidising gas flow, contributing to this relationship. Although this approximation is a simplified description of the more complete dynamics, it embodies the key processes of these flows: the fluidised currents are granular flows in which the fluidisation affects the normal support of the layer.
The experimental measurements provide encouraging support for the model. For example, without tuning through empirical factors, the predictions are quite close to the measured flow depths and flow speed in both the uniform steady state and the transient state as it becomes established. Additionally, when applied to flows along horizontal channels, the model is able to predict the unsteady motion to reveal both the progressive deceleration and the growth in flow depth. There are, however, some systematic features in the measurements that are not reproduced in the model. Perhaps the most significant of these is the decrease in particle velocity towards the top of the layer. This feature is absent from the model and presumably corresponds to particles in the ‘free-board’ of the fluidised layer (i.e. the region above the dense current within which the volume fraction of the particles is reduced). Such a dilute layer is subject to slightly different dynamical interactions: the role of particle collisions becomes much reduced and the particles may saltate, and are potentially intermittently suspended above the denser layer below. The model predictions are also dependent upon the material properties that characterise the collisions between the particles and the boundaries. These can be difficult to measure directly, but the specularity coefficient, , and the boundary coefficient of restitution, , only play a significant role for a relatively thin boundary layer in the dynamical regime considered in this study. It is arguable that the boundary conditions require further research to refine and sharpen their formulation.
One important feature that emerges from the modelling framework is the determination of the volume fraction of particles in the fluidised current. Here we have assumed that the flows are relatively dense and that the Ergun equation provides an appropriate representation of the volume fraction dependence of the drag due to the fluidising gas flow. Other expressions could easily be used in its place (see Nott & Jackson, 1992; Agrawal et al., 2001; Oger & Savage, 2013). However, perhaps of greater significance is whether there are ‘bubbles’, or inhomogeneities in the volume fraction within the fluidised current. Patches of increased voidage locally provide paths through which the fluidising gas can more readily flow and thus it is possible for the layer to exhibit fluctuations or instabilities on relatively rapid timescales. Since the local volume fraction affects the mobility of the flowing layer, one might expect fluctuations in volume fraction and velocity to be correlated and consequentially to influence the bulk dynamics. Bubbles may also affect the particle volume fraction in the bed. The classical model of fluidised beds (Toomey & Johnstone, 1952) proposes that all the gas in excess of that necessary to fluidise the particles forms bubbles so that the particle volume fraction in the bulk of the flow is insensitive to : this is contrast with (45). The lack of dependence on of the experimental scaled velocity profiles in figure 13 also suggests that may change less with conditions than might be expected. In contrast to static fluidised beds, the stability of fully developed fluidised flow down inclines has not been assessed (Jackson, 2000) and this appears to be an interesting topic for future research. Indeed it is intriguing that linear shear flows of unfluidised granular materials appear to exhibit transient linearised growth, but asymptotic stability (Savage, 1992; Scmid & Kytomaa, 1994). It would be interesting to investigate whether these properties carry over to sheared fluidised motions.
Our modelling framework and experimental methods could be extended to a number of related flow problems. First, one could investigate fluidised currents that are generated by instantaneous or non-sustained releases, that are not fully fluidised and for which the fluidising gas flow is localised to the region close to the source. These flows would be largely unsteady and, in situations where the fluidisation is not maintained, would introduce additional mechanisms for generating resistive shear stresses as the contact friction begins to become important. Non-monodisperse granular materials would also be interesting to investigate because the onset of full fluidisation is dependent upon grain-size and the proportion of each particle component (Formisani, 1991). It is possible that mixtures of particles segregate according to size and generate an inhomogeneous flowing current in terms of composition and therefore, the average volume fraction (). Finally, we comment that liquid-fluidised systems may pose additional challenges since it is likely that viscous forces at the particle scale are non-negligible and that collisions are strongly affected by lubrication pressure in the the fluid between particles.
Although direct applications have not been the focus of our study, our model formulation could be readily applied to larger-scale flows, either in industrial contexts or in nature. There are a number of practical implications of our results. For example, the transport of granular materials when they are fluidised is likely to be more efficient on even shallowly inclined surfaces than on horizontal surfaces. In addition, the transport is unlikely to be greatly improved by an increase in the gas flow rate, , once the granular materials are fully fluidised. This is because it only directly affects the average volume fraction, , and for practical materials can strongly depend on their characteristics (e.g. the bubble-free expansion seen in small, light Geldart (1973) group A particles), as can the value of the coefficient of restitution .
The assumption at the heart of the modelling framework is that inelastic particulate collisions generate stresses that provide the resistance to motion and this is likely to be the case for larger-scale flows. Of particular note is that for steady flow, no assumption is made for the relative importance of inertial and resisting forces and hence the resulting model is valid for flows of arbitrary scale. The results show that, at least for some granular flows, full understanding of their nature can only be reached if full account is taken of both the interactions between particles and those between particles and the interstitial fluid.
Acknowledgements
The authors thank three anonymous reviewers who helped to improve our manuscript, as well as O. Pouliquen for his handling of it. This work was funded from a grant under the UK NERC Environmental Mathematics and Statistics programme (NER/S/E/2004/12600). This research was supported also in part by the National Science Foundation under Grant No. NSF PHY11-25915 and AJH also acknowledges support from Max Planck Institute for the Physics of Complex Systems (Two-Phase Continuum Models for Geophysical Particle-Fluid Flows). MAG carried out part of this work while holding a University Research Fellowship provided by the Institute of Advanced Study at the University of Bristol. This paper is LabEx Clervolc contribution no. 259.
Appendix A Extended kinetic theory
In this appendix we analyse the consequences for the predicted flow field of employing the extended kinetic theory proposed by Jenkins (2007) and recently used to compute unfluidised flow down inclined planes by Jenkins & Berzi (2010, 2012) and Berzi (2014). In essence, the extension to kinetic theory is based upon the realisation that at higher concentrations, particles begin to form structures in the flow that have a correlation length in excess of their own diameter. Thus the rate of dissipation is reduced - and in terms of the expression of the evolution of granular temperature (17), the dissipation term is now given by . Jenkins (2007) suggested a phenomenological model for the length, , in which its magnitude is proportional to the rate of compression that occurs along at least one axis in shear flows and inversely proportional to the agitation (the granular temperature) that can destroy these structures. Thus in dimensional form for simple shear flows , Jenkins & Berzi (2010) propose
[TABLE]
where is a dimensionless constant of order unity (often ). Jenkins & Berzi (2010) validate this formulation empirically for unfluidised granular flows. We are not aware of any studies that have tested formulae for fluidised flows, but we can nevertheless employ this formulation (98) to compute profiles of the volume fraction of particles, the velocity field and the granular temperature for typical parameter values used in this study (figure 28). For a dimensionless fluidising gas flow rate, equal to and a slope of , we find negligible differences in the profiles apart from very close to the base of the flow. Moreover the dimensionless volume flux per unit width for the ‘standard’ kinetic theory , while for the extended kinetic theory .
For more weakly fluidised flows, there can be a significant difference between the predictions of the two theories, because in these situations the concentration of particles is higher and thus exceeds unity in many parts of the flow. For example when and , we find that extended kinetic theory predicts more energetic and faster moving flows (see figure 29). For these parameter values, the dimensionless volume flux, for the ‘standard’ kinetic theory, whereas for the extended kinetic theory. The flows that we consider in this study are more strongly fluidised than this example and thus we find it unnecessary to include this phenomenon in our analysis in the main body of this paper because it introduces negligible difference to the computed flow.
Appendix B The effects of side-wall stresses
In this appendix we analyse the effects of side-wall resistance on the motion of shallow fluidised flows down inclined channels (see §5) and derive the first-order correction to the prediction of the front speed for flows that are unaffected by side walls. We show that the reduction in front speed is proportional to , where is the scale depth of the current given by (61) and is the channel breadth.
The downslope flows studied experimentally are realised within a channel, the width of which is usually greater than the flow depth, but not far in excess of the depth. Thus, it is feasible that side wall stresses may play a non-negligible role in the overall dynamics and may further retard the motion. Indeed for flows along horizontal surface for which the flow depth is much greater than the width of the channel, we postulate that the side wall stresses may even play a dominant role in the resisting the driving forces (see §6).
We analyse the motion in a channel of width in a regime for which the volume fraction is spatially uniform and given by . On depth- and width-averaging the streamwise balance of momentum (54), we find that the dimensional governing equation is given by
[TABLE]
where denotes the basal shear stress and denotes the side wall stresses.
In general, to include these side wall effects, even if the flow had adjusted to a local balance that was independent of the streamwise coordinate, we would have to resolve the variations of the dependent fields in the plane (see, for example, Oger & Savage, 2013). In this subsection we take a different strategy and develop a model which is appropriate to the regime , where is the scale depth of the flow and given by (61). In this regime, we treat the granular temperature and flow field as predominantly varying with the distance from the basal boundary and assume that these flow fields adopt the form established in §4.1.5 (see Jenkins & Berzi, 2010). This approach was used above to derive the depth-averaged model above, but is now generalised to include lateral gradients in order to model the side wall stresses.
Adopting the dimensionless variables using the scales of (55), we estimate the velocity gradient at the side wall , where is a dimensionless constant of order unity. We may then compute the depth and width averages to deduce a governing equation for a travelling wave solution that features the additional stresses due to the side walls (cf. (63)); it is given by
[TABLE]
The final term of (100) represents the extra stress due to the side walls. Then using the uniform conditions far from the front we deduce that
[TABLE]
Thus the speed, , is reduced by the action of the side-wall stresses and in the regime , we find that
[TABLE]
The value of the coefficient of the second term of the expansion is 0.2.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Adrian (1991) Adrian, R. J. 1991 Particle-imaging techniques for experimental fluid-mechanics. Annual Rev. Fluid Mech. 23 , 261–304.
- 2Adrian & Westerweel (2010) Adrian, R. J. & Westerweel, J. 2010 Particle Image Velocimetry . Cambridge University Press.
- 3Agrawal et al. (2001) Agrawal, K., Loezos, P. N., Syamlal, M. & Sundaresan, S. 2001 The role of meso-scale structures in rapid gas-solid flows. J. Fluid Mech. 445 , 151–185.
- 4Altantzis et al. (2015) Altantzis, C., Bates, R. B. & Ghoniem, A. F. 2015 3d Eulerian modeling of thin rectangular gas-solid fluidized beds: estimation of the specularity coefficient and its effects on bubbling dynamics and circulation times. Powder Technology 270A , 256–270.
- 5Berzi (2014) Berzi, D 2014 Extended kinetic theory applied to dense, granular, simple shear flows. Acta Mech. 225 , 2191–2198, doi:10.1007/s 00707-014-1125-1.
- 6Bokkers et al. (2004) Bokkers, G. A., van Sint Annaland, M. & Kuipers, J. A. M. 2004 Mixing and segregation in a bidisperse gas-solid fluidised bed: a numerical and experimental study. Powder Technol. 140 (3), 176–186.
- 7Botterill & Abdul-Halim (1979) Botterill, J. S. M. & Abdul-Halim, B. H. 1979 The open-channel flow of fluidized solids. Powder Technol. 23 (1), 67–78.
- 8Botterill & Bessant (1973) Botterill, J. S. M. & Bessant, D. J. 1973 The flow properties of fluidized solids. Powder Technol. 8 (5–6), 213–222.
