4.1.1 Inter-phase drag and constitutive equations

The drag on the solid phase due to the fluidising gas flow is given by $\unicode[STIX]{x1D6FD}(\boldsymbol{u}-\boldsymbol{v})$ , where the drag coefficient, $\unicode[STIX]{x1D6FD}$ , may be written

where $f_{0}$ and $f_{0}^{\ast }$ are given in table 1 (Ergun Reference Ergun1952). 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 $(Re\equiv \unicode[STIX]{x1D70C}_{g}w_{g}d/\unicode[STIX]{x1D707}_{g}<10)$ ; 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. Reference Agrawal, Loezos, Syamlal and Sundaresan2001; Oger & Savage Reference Oger and Savage2013) and these could replace (4.8) within this modelling framework.

Table 1.

The constitutive laws for granular kinetic theory applied to fluidised systems. Here $e$ denotes the coefficient of restitution characterising collisions between particles; $e_{w}$ is the coefficient of restitution between the walls and the particles; $\unicode[STIX]{x1D719}_{m}$ is the volume fraction at maximum packing; $\unicode[STIX]{x1D713}$ is the specularity coefficient (after Johnson & Jackson Reference Johnson and Jackson1987); and $g_{0}$ is the radial basis function which accounts for particle packing (4.15). The coefficient $c^{\ast }=32(1-e)(1-2e^{2})/(81-17e+30e^{2}(1-e))$ (Garzo & Dufty Reference Garzo and Dufty1999).  $f_{0}$ contributes to the inter-phase drag in the viscous regime and $f_{0}^{\ast }$ in the inertial regime (see, for example, van der Hoef, Beetstra & Kuipers Reference van der Hoef, Beetstra and Kuipers2005).  $f_{1}$ , $f_{2}$ model the volume fraction dependence in the collisional contributions to stresses (Garzo & Dufty Reference Garzo and Dufty1999);  $f_{3}$ , $f_{4}$ and $f_{4}^{\ast }$ model contributions to the granular energy balance (Garzo & Dufty Reference Garzo and Dufty1999).  $f_{5}$ to $f_{7}$ determine boundary conditions at the base of the flow: $f_{5}$ contributes to the boundary condition for momentum balances and $f_{6}$ and $f_{7}$ to that for fluctuation energy (see Johnson & Jackson Reference Johnson and Jackson1987; Johnson et al. Reference Johnson, Nott and Jackson1990).

It was noted above that particle interactions are dynamically important because of the momentum transfer arising from particle collisions (Lun et al. Reference Lun, Savage, Jeffrey and Chepurniy1984; Garzo & Dufty Reference Garzo and Dufty1999). Here we examine the collisional stresses and follow Nott & Jackson (Reference Nott and Jackson1992) and Agrawal et al. (Reference Agrawal, Loezos, Syamlal and Sundaresan2001) 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, $T$ , which measures the fluctuations of velocity about the mean, the volume fraction of solids, $\unicode[STIX]{x1D719}$ , and the coefficient of restitution $e$ , 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

where $f_{1}$ and $f_{2}$ are dimensionless functions given in table 1. In this study, we employ the constitutive formulae derived by Garzo & Dufty (Reference Garzo and Dufty1999) and recently used for modelling dense avalanches by Jenkins & Berzi (Reference Jenkins and Berzi2010).

Granular temperature may be generated and ‘conducted’ via the flow processes and dissipated in the collisions. Following Lun et al. (Reference Lun, Savage, Jeffrey and Chepurniy1984) and Garzo & Dufty (Reference Garzo and Dufty1999) amongst others, these effects are encompassed in the following expression of energy balance within the flow

where the flux of granular temperature is given by

and $f_{3}$ , $f_{4}$ and $f_{4}^{\ast }$ 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 Reference Koch and Sangani1999).

Following the formulation of Agrawal et al. (Reference Agrawal, Loezos, Syamlal and Sundaresan2001), dissipation by viscous processes is much smaller than dissipation through inelastic collisions when

Furthermore, the generation of granular temperature by viscous processes is much smaller than that by granular interactions when

The constitutive laws, $f_{1}-f_{4}^{\ast }$ , as well as those involved the boundary conditions ( $f_{5}-f_{7}$ , see § 4.1.3) feature the radial basis function, $g_{0}(\unicode[STIX]{x1D719})$ . Various authors have suggested forms for $g_{0}$ and we employ an expression that is close to the suggestion of Vescovi et al. (Reference Vescovi, Berzi, Richard and Brodu2014), 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 (Reference Torquato1995)) and following Vescovi et al. (Reference Vescovi, Berzi, Richard and Brodu2014) we write

where the weighting function is given by

Thus, when $\unicode[STIX]{x1D719}<\unicode[STIX]{x1D719}_{\ast }$ the radial basis function is given by the formula proposed by Carnahan & Starling (Reference Carnahan and Starling1969), but it exceeds this value when $\unicode[STIX]{x1D719}_{\ast }<\unicode[STIX]{x1D719}$ and diverges as maximum packing is approached. Vescovi et al. (Reference Vescovi, Berzi, Richard and Brodu2014) suggest that $n=2$ and that $\unicode[STIX]{x1D719}_{\ast }=0.4$ . While this choice ensures that $g_{0}$ and its derivative are continuous at $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{\ast }$ , the second derivative is discontinuous. This is problematic for the system of differential equations that we will integrate numerically; therefore, we employ the values $n=3$ and $\unicode[STIX]{x1D719}_{\ast }=0.4$ , which ensure that $g_{0}$ is sufficiently smooth. Moreover, our expression (4.16) with these values is close to those proposed by Torquato (Reference Torquato1995) and Vescovi et al. (Reference Vescovi, Berzi, Richard and Brodu2014) and appears to match the simulation data adequately.

The energetic balance encompassed in (4.11) 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 $\unicode[STIX]{x1D70C}_{s}f_{3}T^{3/2}/d$ ). Recently, however, Jenkins (Reference Jenkins2007) has suggested that at relatively high concentrations, clusters of particles begin to form and thus the correlation length increases to $L_{c}$ ( ${>}d$ ) and then the rate of dissipation is given by $\unicode[STIX]{x1D70C}_{s}f_{3}T^{3/2}/L_{c}$ . This extended kinetic theory has been applied to unfluidised flows of grains down inclined planes by Jenkins & Berzi (Reference Jenkins and Berzi2010, Reference Jenkins and Berzi2012), where an empirical formula for $L_{c}$ , 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, $e$ , and between the particles and the boundary, $e_{w}$ , and the specularity coefficient $\unicode[STIX]{x1D713}$ , which governs the dynamic interaction between the particles and the bottom surface (4.32). These parameters are relatively difficult to measure directly.

The coefficient of restitution, $e$ , plays an important role in continuum models and in discrete particle (or element) models (DPMs), which endeavour to calculate the motion of large ensembles of particles and to resolve individual particle collisions. In continuum models, the difference of $e$ from unity is proportional to the rate at which the collisions dissipate energy (see the definition of $f_{3}$ 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, $e=0.90$ and $0.97$ respectively in the studies of Goldschmidt et al. (Reference Goldschmidt, Beetstra and Kuipers2004) and van der Hoef et al. (Reference van der Hoef, van Sint Annaland, Deen and Kuipers2008)). These values are close to measured values of discrete collisions (see, for example, Kharaz, Gorham & Salman Reference Kharaz, Gorham and Salman2001). When used in kinetic theory models, commonly adopted values of $e$ are rather lower and Jenkins & Zhang (Reference Jenkins and Zhang2002) suggest a means by which 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 3 mm diameter, the measured data of Foerster et al. (Reference Foerster, Louge, Chang and Allia1994) correspond to an effective coefficient of $e=0.85$ if the method of Jenkins & Zhang (Reference Jenkins and Zhang2002) 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 $e_{w}=0.75$ , 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 (Reference Johnson and Jackson1987) and Johnson et al. (Reference Johnson, Nott and Jackson1990). 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

What this means physically is that solid-phase slip is allowed and stress is transmitted in the downslope 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 $\unicode[STIX]{x1D713}$ $(0<\unicode[STIX]{x1D713}<1)$ . Furthermore, fluctuation energy (granular temperature) is generated at the bottom surface and potentially dissipated by inelastic collisions, encompassed through a coefficient of restitution, $e_{w}$ .

At the top of the current $z=h$ , there are the free-surface boundary conditions that fluid shear and normal stresses vanish, given by

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. Reference Johnson, Nott and Jackson1990). Thus we enforce

The boundary conditions (4.17)–(4.19) are of the same character as those employed by researchers in other flow regimes (see, for example, Jenkins & Berzi Reference Jenkins and Berzi2010) 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 downslope momentum equations of each phase (4.2) and (4.4) to eliminate the inter-phase drag so that

In this expression, the key driving force is the downslope 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,

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 (4.11),

whence, the scaling for the velocity field is given by

It is now convenient to introduce dimensionless variables, given by

This set of scalings for the dependent variables of fluidised flows differs from those for flows of unfluidised, collisional granular media (Woodhouse et al. Reference Woodhouse, Hogg and Sellar2010). 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 $gh\cos \unicode[STIX]{x1D703}$ , which is considerably larger than the estimate deduced here (4.24) unless the motion is along relatively steep inclines (i.e. when $\tan \unicode[STIX]{x1D703}\sim 1$ ). 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:

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,

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: (4.13) is satisfied when $\unicode[STIX]{x1D6FF}^{2}/St\ll 1$ and (4.14) when $W_{g}\ll 1$ .

Table 2.

Range of values of the physical parameters in the experiments and the dimensionless groups derived from them as defined by (4.25).

We have five governing equations: mass conservation (4.1), downslope fluid momentum conservation (4.2), the combined normal momentum equation from which the fluid pressure has been eliminated (4.7), the downslope solids momentum equation (4.4) and the equation for the conservation of granular temperature (4.11). Non-dimensionalised, these become

where ${\mathcal{U}}$ measures the magnitude of the dimensionless relative velocity between the phases and is given by ${\mathcal{U}}^{2}=S^{2}(\hat{u} -\hat{v})^{2}/(W_{g}St)^{2}+{\hat{w}}^{2}$ . From (4.17), the dimensionless boundary conditions at the base $(\hat{z}=0)$ are given by

while from (4.18), we enforce at the top surface $(\hat{z}=1)$

Figure 2.

The volume fraction, $\unicode[STIX]{x1D719}(\hat{z})$ , velocity of the solid phase, $\hat{v}(\hat{z})$ , and the granular temperature, $\hat{T}(\hat{z})$ , as functions of the dimensionless depth within the current for parameter values $R=10^{-3}$ , $\unicode[STIX]{x1D713}=0.5$ , $\unicode[STIX]{x1D719}_{m}=0.63$ , $e=0.85$ , $e_{w}=0.75$ , $S=0.1$ , $St=10^{3}\unicode[STIX]{x1D6FF}^{2}$ , $W_{g}=10^{-3}$ and (i) $\unicode[STIX]{x1D6FF}=0.1$ , (ii) $\unicode[STIX]{x1D6FF}=0.01$ and (iii) $\unicode[STIX]{x1D6FF}=0.001$ . Also plotted are the asymptotic solutions (dotted lines), although these are often overlain by the full solution.

Figure 3.

The volume fraction, $\unicode[STIX]{x1D719}(\hat{z})$ , velocity of the solid phase, $\hat{v}(\hat{z})$ , and the granular temperature, $\hat{T}(\hat{z})$ , as functions of the dimensionless depth within the current for parameter values $R=10^{-3}$ , $\unicode[STIX]{x1D713}=0.5$ , $\unicode[STIX]{x1D719}_{m}=0.63$ , $e=0.85$ , $e_{w}=0.75$ , $S=0.1$ , $St=0.1$ , $\unicode[STIX]{x1D6FF}=0.01$ , (i) $W_{g}=0.5\times 10^{-3}$ , (ii) $W_{g}=1.0\times 10^{-3}$ and (iii) $W_{g}=2\times 10^{-3}$ . Also plotted are the asymptotic solutions (dotted lines).

Figure 4.

The volume fraction, $\unicode[STIX]{x1D719}(\hat{z})$ , velocity of the solid phase, $\hat{v}(\hat{z})$ , and the granular temperature, $\hat{T}(\hat{z})$ , as functions of the dimensionless depth within the current for parameter values $R=10^{-3}$ , $\unicode[STIX]{x1D713}=0.3$ , $\unicode[STIX]{x1D719}_{m}=0.63$ , $e=0.85$ , $e_{w}=0.75$ , $St=0.1$ , $\unicode[STIX]{x1D6FF}=0.01$ , $W_{g}=10^{-3}$ , (i) $S=0.1$ ; (ii) $S=0.2$ and (iii) $S=0.3$ . Also plotted are the asymptotic solutions (dotted lines).

The system of governing differential equations (4.27)–(4.31) and boundary conditions (4.32)–(4.33) form a seventh-order differential boundary value problem. We use (4.27) to eliminate ${\hat{w}}$ in favour of $1/(1-\unicode[STIX]{x1D719})$ and we also evaluate $\text{d}^{2}\unicode[STIX]{x1D719}/\text{d}z^{2}$ in (4.31) by explicitly differentiating (4.29). 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, $\unicode[STIX]{x1D719}(\hat{z})$ , the granular temperature, $\hat{T}(\hat{z})$ , and the velocity of the solid phase, $\hat{v}(\hat{z})$ , are plotted in figures 2–4 for various values of the governing dimensionless parameters. We do not plot the gas velocity, $\hat{u} (\hat{z})$ , because outside of thin basal boundary layers, it is indistinguishable from the solids velocity, $\hat{v}(\hat{z})$ . 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 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 $\unicode[STIX]{x1D6FF}$ 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 (4.24).

It is of particular interest to evaluate the dimensionless flux of solids per unit width and the average concentration of particles, respectively given by

and these are plotted as functions of the dimensionless parameters in figures 5–7.

From figure 5, we observe that the dimensionless volume flux per unit width, $\hat{q}$ , and the average volume fraction, $\overline{\unicode[STIX]{x1D719}}$ , do not vary strongly with the relative particle size, $\unicode[STIX]{x1D6FF}$ . 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, $W_{g}$ , are rather more subtle (see figure 6). Increasing $W_{g}$ 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, $\hat{q}$ , does not vary monotonically with $W_{g}$ . Indeed, for the parameters in figure 6, it attains a maximum at a dimensionless gas flow rate $W_{g}=2.3\times 10^{-3}$ and at that value of $W_{g}$ , $\unicode[STIX]{x1D719}=0.43$ . 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 $W_{g}$ and $St$ , but maintain $Re$ constant (see (4.25) and (4.26)). 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 $\hat{q}$ to be maximised at some finite value of $S$ . (For the parameters analysed in figure 7, the local maximum in the flux occurs at $S=0.22$ .)

Figure 5.

The dimensionless volume flux per unit width transported by the flowing layer, $\hat{q}$ , and the depth-averaged volume fraction, $\overline{\unicode[STIX]{x1D719}}$ , as functions of the relative particle size for parameter values $R=10^{-3}$ , $\unicode[STIX]{x1D713}=0.50$ , $\unicode[STIX]{x1D719}_{m}=0.63$ , $e=0.85$ , $e_{w}=0.75$ , $S=0.10$ , $St=10^{3}\unicode[STIX]{x1D6FF}^{2}$ and $W_{g}=10^{-3}$ . Also plotted are the asymptotic solutions (dashed) and the simple approximate solutions (dotted).

Figure 6.

The dimensionless volume flux per unit width transported by the flowing layer, $\hat{q}$ , and the depth-averaged volume fraction, $\overline{\unicode[STIX]{x1D719}}$ , as functions of the dimensionless strength of the fluidising gas flow, $W_{g}$ , for parameter values $R=10^{-3}$ , $\unicode[STIX]{x1D713}=0.50$ , $\unicode[STIX]{x1D719}_{m}=0.63$ , $e=0.85$ , $e_{w}=0.75$ , $S=0.10$ , $St=0.10$ and $\unicode[STIX]{x1D6FF}=10^{-2}$ . Also plotted are the asymptotic solutions (dashed) and the simple approximate solutions (dotted).

Figure 7.

The dimensionless volume flux per unit width transported by the flowing layer, $\hat{q}$ , and the depth-averaged volume fraction, $\overline{\unicode[STIX]{x1D719}}$ , as functions of the slope of the underlying boundary, $S=\tan \unicode[STIX]{x1D703}$ , for parameter values $R=10^{-3}$ , $\unicode[STIX]{x1D713}=0.50$ , $\unicode[STIX]{x1D719}_{m}=0.63$ , $e=0.85$ , $e_{w}=0.75$ , $\unicode[STIX]{x1D6FF}=10^{-2}$ , $W_{g}=9.95\times 10^{-4}(1+S^{2})^{1/2}$ , $St=0.317S^{1/2}(1+S^{2})^{-1/4}$ . Also plotted are the asymptotic solutions (dashed) and the simple approximate solutions (dotted). Note the local maximum at $S=0.21$ .

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 $\unicode[STIX]{x1D6FF}\ll 1$ and $R\ll 1$ . This regime will have a widespread validity as $d\ll h$ in order to use a continuum approach, and for gas–solid flows $\unicode[STIX]{x1D70C}_{g}\ll \unicode[STIX]{x1D70C}_{s}$ . From (4.28), we note that to leading order and away from boundaries, the downslope velocities of the two phases must be equal ( $\hat{u} =\hat{v}+O(1)$ ). Furthermore, from (4.31) there is a local balance between granular temperature production and dissipation such that

provided the volume fraction of particles is not too small (i.e. $\unicode[STIX]{x1D719}\gg \unicode[STIX]{x1D6FF}$ , so that the ‘conductive’ effects of the granular temperature remain negligible).

The governing equations for the normal and perpendicular momentum balances, (4.29) and (4.30), are then given by

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 (4.35). When $\unicode[STIX]{x1D6FF}\ll 1$ , 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 $\hat{T}(1)=0$ and the volume fraction at the base is given by $\unicode[STIX]{x1D719}(0)=\unicode[STIX]{x1D719}_{0}$ , which is determined by substituting for $\unicode[STIX]{x2202}\hat{v}/\unicode[STIX]{x2202}\hat{z}$ from (4.35) into (4.32),

where the constitutive functions $(f_{i})$ are evaluated at $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{0}$ . The basal volume fraction is thus a function of $e$ , $e_{w}$ , $\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D719}_{m}$ .

Rearranging (4.36) and (4.37) and denoting $f=(f_{1}f_{3})^{1/2}$ , we find that

where $\dot{~}$ denotes differentiation with respect to $\unicode[STIX]{x1D719}$ . It is straightforward to integrate numerically these coupled first-order equations subject to the boundary conditions $\hat{T}(1)=0$ and $\unicode[STIX]{x1D719}(0)=\unicode[STIX]{x1D719}_{0}$ . The solutions are plotted in figures 2–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 $\unicode[STIX]{x1D719}(z)=\overline{\unicode[STIX]{x1D719}}$ when

where the constitutive functions are evaluated at $\unicode[STIX]{x1D719}=\overline{\unicode[STIX]{x1D719}}$ . In this case, the temperature gradient is constant, $\unicode[STIX]{x2202}\hat{T}/\unicode[STIX]{x2202}\hat{z}=-\unicode[STIX]{x1D706}$ , with $\unicode[STIX]{x1D706}=\overline{\unicode[STIX]{x1D719}}/f$ . This solution is ‘attracting’ in the sense that trajectories in the phase plane $(\unicode[STIX]{x1D719}(\hat{z}),\hat{T}(\hat{z}))$ approach it when

which in turn demands that $\overline{\unicode[STIX]{x1D719}}>\overline{\unicode[STIX]{x1D719}}_{c}(e)$ and that $S<S_{c}(e)$ (see figure 8 c). If these inequalities are not held then the reduced system evolves towards a state different from a uniform volume fraction $(\unicode[STIX]{x1D719}=\overline{\unicode[STIX]{x1D719}})$ and may not admit solutions at all. Physically, when $S>S_{c}$ , the dissipation of granular temperature, here encapsulated through a coefficient of restitution $e$ , 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 (4.36)). When $e\lessapprox 0.85$ , 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).

Figure 8.

(a) The limiting volume fraction, $\overline{\unicode[STIX]{x1D719}}_{c}$ and slope, $S_{c}$ , for which the uniform volume fraction is the ‘attracting’ solution as functions of the coefficient of restitution; (b) the volume flux per unit width, $\hat{q}$ , (c) the average volume fraction, $\overline{\unicode[STIX]{x1D719}}$ , as functions of the slope, $S$ , for $W_{g}=10^{-3}$ and varying values of the coefficient of restitution and (d) the product of $\overline{\unicode[STIX]{x1D719}}$ and the mobility factor, $F$ , as a function of $\overline{\unicode[STIX]{x1D719}}$ for various values of the coefficient of restitution.

When the dimensionless fluidisation velocity $W_{g}$ and the slope $S$ are set, the average volume fraction within the current, $\overline{\unicode[STIX]{x1D719}}$ , can be calculated using (4.41). Figure 8(a) shows the effect of $e$ on the curves of $\overline{\unicode[STIX]{x1D719}}$ as a function of slope $S$ when the fluidisation flow is constant ( $W_{g}=10^{-3}$ ). The curves in this plot are continued up to the maximum value of the slope, $S_{c}(e)$ , 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 8(b), for a given slope, $S$ , and fluidisation gas flow rate, $W_{g}$ , 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 $\hat{z}=1$ to leading order, we find for the simple approximate solution with uniform volume fraction that the granular temperature is given by

and the velocity field of the solid phase is given by

where $F=(\overline{\unicode[STIX]{x1D719}}^{2}f_{3}/f_{1}^{3})^{1/4}$ . The scaled slip velocity at the wall $Ff_{5}\unicode[STIX]{x1D6FF}$ can be added to (4.44), but when $\unicode[STIX]{x1D6FF}\ll 1$ , it is negligible. The velocity profile (4.44) is similar to the dimensionless ‘Bagnold’ velocity profile up to the factor $F$ , 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 $(\unicode[STIX]{x1D6FF}\ll 1)$ , the approximate solution yields

and this is plotted in figures 5–7, once again illustrating the utility of this asymptotic solution. The quantity $\overline{\unicode[STIX]{x1D719}}F$ thus plays a crucial role in determining the dimensionless flux, $\hat{q}$ and in figure 8(d), we plot its dependence on volume fraction for a range of values of $e$ . We note that $\overline{\unicode[STIX]{x1D719}}F$ is maximised for $\unicode[STIX]{x1D719}\approx 0.41$ (with the precise value weakly dependent on $e$ ) and vanishes both when $\unicode[STIX]{x1D719}$ 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 $(\unicode[STIX]{x1D719}\approx 0.41)$ . Finally, we note that a dimensional estimate of the depth of a fluidised current may be obtained

where $q_{0}$ is the dimensional flux per unit width at the source and the effects of slip at the wall has been neglected.

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, $h_{vis}(\equiv h)$ : $h_{vis}$ is greater than the height at which the maximum velocity is attained. Above $h_{vis}$ 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 $h_{max}$ . The depth of the current can fluctuate a little with time (see figures 18 and 19 a). 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 (4.46).

Figure 9.

Image of a granular current with the measured velocity superimposed onto it (solid line). The image is an average-intensity composite of the images used in the PIV measurements. Note the non-zero (slip) velocity at the base ( $v_{0}$ ). The dashed lines indicate the 95 % confidence interval using (2.2). Also shown are the height where the velocity profile drops to zero, $h_{max}$ , the maximum visible height of the current, $h_{vis}$ , and the height of the peak of the velocity profile, $z(v=v_{max})$ . For this figure $\unicode[STIX]{x1D703}=15^{\circ }$ , $Q=42.5~\text{cm}^{3}~\text{s}^{-1}$ , $w_{g}/u_{mf}=1.5$ .

Table 3.

Various estimates of height in steady-state currents. $h_{vis}$ is measured from photographs of the currents such as that in figure 9. It may be compared with the prediction, $h$ , calculated from (4.46) with $F$ based on $\bar{\unicode[STIX]{x1D719}}_{est}$ (see table 4), $e=0.85$ and $\unicode[STIX]{x1D713}=0.50$ . $h_{max}$ is found directly as the height at which particle velocity drops to zero. In all cases $w_{g}=1.5u_{mf}$ .

The asymptotic solution yielded an approximate formula linking height and source volume flux (4.46), and this is shown in figure 10. In this figure, fixed values of $\overline{\unicode[STIX]{x1D719}}$ were used because from (4.41), $\overline{\unicode[STIX]{x1D719}}$ depends on $\unicode[STIX]{x1D703}$ and it was not possible to find a solution for the full experimental range of $\unicode[STIX]{x1D703}$ for a fixed value of $e=0.85$ . The decreasing effect of $\overline{\unicode[STIX]{x1D719}}$ on mobility as its value approaches $\overline{\unicode[STIX]{x1D719}}_{c}=0.40$ reflects the effect it has on mobility $\overline{\unicode[STIX]{x1D719}}F$ shown in figure 8(d). 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.

Figure 10.

The height of the flowing layer as a function of the scaled volume flux, $q_{0}/\sqrt{g\sin \unicode[STIX]{x1D703}}$ , for a range of channel inclinations. The model curves are for (4.46) with $\unicode[STIX]{x1D719}_{m}=0.63$ and $e=0.85$ .

4.2.2 Velocity profiles

The ensemble-averaged velocity profiles for $10^{\circ }$  and $15^{\circ }$  slopes measured half-way 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 $\unicode[STIX]{x1D703}=10^{\circ }$ , with source fluxes $Q=33~\text{cm}^{3}~\text{s}^{-1}$ and $Q=34~\text{cm}^{3}~\text{s}^{-1}$ 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. Reference Woodhouse, Hogg and Sellar2010); 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 $Q=33~\text{cm}^{3}~\text{s}^{-1}$ leads to small estimates of particle volume fraction with an excessive portion of the flow where $h>h_{vis}$ (see below, § 4.2.3. Based on $h_{vis}$ , $\bar{\unicode[STIX]{x1D719}}=0.33$ ; based on $h_{max}$ , $\bar{\unicode[STIX]{x1D719}}=0.22$ ), which seem physically unlikely, and so this experimental run is not reported further.

Figure 12 shows some instantaneous velocity profiles half-way along the tank for shallower slope angles whose motion may not be steady. The velocity profiles for the $3^{\circ }$  and $5^{\circ }$   slopes are similar in character to the averaged profiles for steeper slopes.

Figure 11.

Ensemble-averaged velocity profiles for steady flows for $\unicode[STIX]{x1D703}=10^{\circ }$ and $\unicode[STIX]{x1D703}=15^{\circ }$ . Error bars are 95 % confidence limits calculated from (2.2) with $m=5$ and $n=16$ .

Figure 12.

Measured velocity profiles from experiment for individual fluidised granular flows that are potentially unsteady when they are at shallow angles. Error bars are calculated from (2.1) with values of $c_{1}$ determined from the correlation functions calculated during the PIV analysis.

4.2.3 Particle volume fraction

The expectation from the model is that the particle volume fraction $\unicode[STIX]{x1D719}$ would be approximately uniform within the fluidised currents. It is not possible to analyse the degree to which $\unicode[STIX]{x1D719}$ 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 $h_{vis}$ , $\unicode[STIX]{x1D719}$ 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, $\bar{\unicode[STIX]{x1D719}}$ , by integrating the particle velocity profiles and dividing the measured particle flow rate $\unicode[STIX]{x1D719}_{0}Q$ by the result. The results $\unicode[STIX]{x1D719}_{meas}$ are shown in table 4, using $h_{vis}$ as the overall depth of the current, to give measured volume fraction of the currents ( $\bar{\unicode[STIX]{x1D719}}_{meas}$ ). The measured values of particle volume fraction can be compared with estimated values, $\unicode[STIX]{x1D719}_{est}$ , which have been calculated using (4.41).

There can be quite good agreement between $\bar{\unicode[STIX]{x1D719}}_{est}$ and $\bar{\unicode[STIX]{x1D719}}_{meas}$ despite several inherent uncertainties in their computation. Overall, the values of $\bar{\unicode[STIX]{x1D719}}_{meas}$ were comparable to the values of $\bar{\unicode[STIX]{x1D719}}_{est}$ and with the typical values of $\unicode[STIX]{x1D719}=0.50{-}0.60$ for static fluidised beds (Epstein & Young Reference Epstein and Young1962). In between $h_{vis}$ and $h_{max}$ particles are present, but in practice the fall off of velocity above $h_{vis}$ is sufficiently rapid that this makes very little difference to calculations of $\unicode[STIX]{x1D719}_{meas}$ : if $\unicode[STIX]{x1D719}_{meas}$ is calculated on the basis of the top of the current being $h_{max}$ rather than $h_{vis}$ , then its value decreases by less than 0.02, except when $\unicode[STIX]{x1D703}=15^{\circ }$ and $Q=11.2~\text{cm}^{3}~\text{s}^{-1}$ , when it reduces by 0.09.

Table 4.

Measured estimates for $Q$ and $\bar{\unicode[STIX]{x1D719}}$ from integration of the velocity profiles. Values of $\bar{\unicode[STIX]{x1D719}}_{est}$ are found from solving (4.41), and those for $\bar{\unicode[STIX]{x1D719}}_{meas}$ from integration of the measured velocity profiles up to $h_{vis}$ . The errors are calculated using 95 % confidence limits calculated from (2.2).

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. Reference Goldschmidt, Beetstra and Kuipers2004). Their value is sufficiently badly defined that in their investigations of bubbling fluidised beds of glass particles Altantzis, Bates & Ghoniem (Reference Altantzis, Bates and Ghoniem2015) used values between $10^{-4}$ and $0.5$ and Li, Grace & Bi (Reference Li, Grace and Bi2010) from $0$ to $0.5$ . 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 $\unicode[STIX]{x1D713}$ from the directly measured slip velocities and gradients using the boundary condition (4.17) and the definition of $f_{5}$ in table 1, so that in terms of dimensional variables

The results are shown in table 5, and it can be seen that a measured average value of $\unicode[STIX]{x1D713}$ is 0.28. The values for $\unicode[STIX]{x1D713}$ 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 $\unicode[STIX]{x1D713}$ ; however, despite these uncertainties, the value of $\unicode[STIX]{x1D713}$ is reasonably consistent. The effect of $\unicode[STIX]{x1D713}$ on the velocity profiles is to introduce a slip velocity proportional to $\unicode[STIX]{x1D6FF}f_{5}$ . Even for the relatively large values of $\unicode[STIX]{x1D713}$ estimated here, the magnitude of this dimensionless term is relatively small.

Table 5.

Measured slip velocities at the wall, $v(0)$ , and velocity gradients $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}z|_{z=0}$ with the resulting estimate for slip length and for $\unicode[STIX]{x1D713}$ using (4.47). The slip length is expressed in terms of particle diameters and is defined as $v(0)/\text{d}(\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}z)|_{z=0}$ . $e=0.85$ , $\unicode[STIX]{x1D719}_{m}=0.63$ .

4.2.5 Scaling of the velocity profiles

The measured velocity profiles scaled as $\hat{v}$ and $\hat{z}$ are plotted in figure 13. With the exception of the lowest flow rate when $\unicode[STIX]{x1D703}=15^{\circ }$ , the data collapse 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 $\bar{\unicode[STIX]{x1D719}}$ and hence $F$ , but this is not reflected in the experimental curves for which the scaling appears to eliminate the effect of $S$ . $F$ is also affected by the value of $e$ . Increasing $e$ causes a decrease in the predicted dimensionless velocity: the theoretical curves shift towards the left and the difference between the curves for $\unicode[STIX]{x1D703}=10^{\circ }$ and $15^{\circ }$ becomes less (see inset, figure 13).

Figure 13.

Velocity profiles with measurements scaled as (4.24) compared with the model scaled velocity profile under steady state (4.44) for each angle where $e=0.85$ . The curves correspond to $\bar{\unicode[STIX]{x1D719}}=0.51$ when $\unicode[STIX]{x1D703}=10^{\circ }$ and $\bar{\unicode[STIX]{x1D719}}=0.41$ when $\unicode[STIX]{x1D703}=15^{\circ }$ . $W_{g}=8.09\times 10^{-4}$ when $\unicode[STIX]{x1D703}=10^{\circ }$ and $W_{g}=8.24\times 10^{-4}$ when $\unicode[STIX]{x1D703}=15^{\circ }$ . The inset graph shows the effect of the value of $e$ on the model solutions with the chain-dot curves corresponding to $e=0.95$ and the dotted curves to $e=0.75$ .