A mathematical model for wind-generated particle–fluid flow fields with an application to the helicopter cloud problem

Abstract

We develop a model for the interaction of a fluid flowing above an otherwise static particle bed, with generally the particles being entrained or detrained into the fluid from the upper surface of the particle bed, and thereby forming a fully two phase fluidized cloud above the particle bed. The flow in this large-scale fluidized region is treated as a two-phase flow, whilst the key processes of entrainment and detrainment from the particle bed are treated by examining the local dynamical force balances on the particles in a thin transition layer at the interface between the fully fluidized region and the static particle bed. This detailed consideration leads to the formation of an additional macroscopic boundary condition at this interface, which closes the two-phase flow problem in the bulk fluidized region above. We then introduce an elementary model of the well-known helicopter brownout problem, and use the theory developed in the first part of the paper to fully analyse this model, both analytically and numerically.

1. Introduction

In this paper, we develop a rational (that is, based directly on the fundamental Newtonian laws of mechanics) physically based mathematical model for the generation of wind-driven particle flow fields from a static particle bed. The physical model is composed of three structured regions: the fluidized region in which the particles lifted from the surface of the static particle bed are in suspension in the fluid, forming a fully developed two-phase flow; the key interfacial layer in which the interaction mechanisms between the local near-surface fluid flow and the particle bed enable surface particles to be entrained by the local flow from the static bed of particles into the fluidized region and/or detrained from the local flow in the fluidized region into the static bed of particles; and the static bed region where the particles are stored.

A particular example of such wind-driven particle flow fields is the scenario in which a helicopter may descend towards a bed of sand, the helicopter cloud problem (often referred to as the ‘brownout problem’ in the aero-engineering literature). In this case, the local interaction between the downdraught and swirling flow, generated by the helicopter rotor, and the upper surface of the otherwise static sand bed, entrains sand particles from a thin interfacial layer into the fluidized region, which then flow in the form of a high-velocity particle-laden cloud around the helicopter body, which we refer to as the helicopter cloud, a consequence of which is generally a significant deterioration in visibility for the helicopter pilot – when this becomes too severe, it is referred to as brownout. This problem has received much attention in the engineering literature, but most work to date has been driven primarily by observation and experimentation. Attempts to model the physics of brownout include those by Wachspress et al. (2008), Phillips & Brown (2009), Phillips, Kim & Brown (2011), Govindarajan, Leishman & Gumerov (2013), Porcù et al. (2020), Tan et al. (2021,  2022a,b), Li et al. (2023) and Lin et al. (2024). In addition, a very recent thesis by Rovere (2022) reviews particle tracking simulations together with experiments and high-speed photographic evidence related to the helicopter cloud problem, and provides an excellent survey of the engineering literature in this area. These modelling approaches broadly split into two classes, which we now delineate.

The first approach is purely computational, and treats each particle within the fluid flow as an individual entity, with no effort being made to address the entrainment or detrainment of particles into the particle-laden flow via local interaction of the flow with the upper surface of the otherwise static particle bed. Indeed, the particle bed is removed from consideration, and instead, the particles are simply placed into the flow in a random spatial distribution at an initial reference time. Subsequently, each particle is tracked as it moves within the fluid flow according to its own dynamical equation of motion under the action of the locally induced fluid interaction forces and gravity (the following, amongst many others, are papers adopting this principal strategy as a key modelling component: Govindarajan et al. 2013; Porcù et al. 2020; Tan et al. 2021,  2022a,b; Li et al. 2023; Lin et al. 2024). Although this approach does give rational information about particulate behaviour once the particles have made their way into the fluidized flow, the lack of consideration of the entrainment and detrainment process at the interface between the fluidized region and the static particle bed is a very severe drawback for this modelling approach – it is this interaction process that is the fundamental and key process in initiating and driving the whole phenomenon of the particle cloud. In addition, given that particle clouds in the fluidized region are generally concentrated rather than dilute, with a low fluid voidage (high particle concentration), the approach of treating each particle in the fluidized region as an individual entity is exceptionally inefficient in computational time, and a two-phase flow approach is much more natural and efficient.

The second modelling approach is to introduce a continuum particle density field in the fluidized region (measuring particle volume per unit spatial volume of the two-phase flow), and then to postulate that this field satisfies a suitable advection–diffusion partial differential equation (PDE) throughout the flow field. This approach does go some way to address the entrainment and detrainment of particles from the static particle bed. Specifically, this is done in a phenomenological way by the introduction of a localized particle mass source term into the advection–diffusion PDE. This source term is localized in space, so that it acts only in a thin neighbourhood of the interface between the fluidized region and the static particle bed, and is designed to represent the localized input/ouput of particles from the static bed into the fluidized region (see, in particular, Phillips et al. 2011). The nature of this source term is purely phenomenological, empirically based on the very particular situation under immediate consideration, and must be recalibrated in every specific example; as such, the ability of this approach to capture, in general, a decent representation of the key rational mechanism of entrainment/detrainment at the interface of the fluidized region and the static bed is very limited, and depends entirely, and critically, upon the specific choice of this term. This serious limitation is addressed well in the paper of Phillips et al. (2011). In their conclusions, they state that further developments in the theory for the brownout problem depend crucially on an understanding ‘particularly of how sediment leaves the ground and enters the flow around the helicopter in the first place’. Moreover, these authors state (Phillips et al. 2011, p. 126) that ‘the physics of sediment entrainment from the ground … is no doubt the area which could most benefit from further research and, indeed, where a breakthrough in understanding could contribute most to improving our confidence in the ability of predictive methods to capture the detailed mechanisms at the origin of the brownout cloud’. The key message is that although this approach does address local particle entrainment and detrainment from the static particle bed, a thorough understanding and implementation of the fundamental interactive mechanics governing these two principal processes, and their coupling with the fluidized region, is essential for improving predictability in the problem. It is this serious defect that we address from fundamental first principles in this paper, which leads to a natural macroscopic boundary condition on void fraction to be applied at the interface of the fluidized region and the static bed, which closes the continuum scale problem two-phase flow in the fluidized region.

In more general terms, the development of fluidized clouds in gas or liquid flows above otherwise static particle beds has many applications – for example, clouds generated by the motion of desert vehicles, motion of sub-sea vehicles close to the ocean bed, and large-scale particle clouds raised by localized desert or sub-surface oceanic storms and disturbances. Indeed, the motion of desert dunes via bed-load and suspended-load sand transport, and the related physical processes, were originally considered in the classic text of Bagnold (1941), in which it was acknowledged that for a generally applicable predictive theory for both of these transport phenomena, a detailed understanding of the local interactive mechanisms in the layer adjacent to the interface between the fluid and the otherwise static particle bed would be essential. In the absence of such an understanding Bagnold (1941), whilst acknowledging their shortcomings, resorted to the development of empirically fitted sediment transport relations. The physical aspects of wind-blown sand have been addressed more recently in the extensive review article of Kok et al. (2012).

In this paper, we address these issues directly and generically. We consider the general situation of an incompressible fluid overlying an otherwise static particle bed, and allow for the general motion of the fluid to interact with the particle bed in a thin transition layer, in which the local force balance between the fluid flow and the incipiently fluidized particles is fully accounted for on the thickness scale of this thin transition region, which rationally accounts for both the entrainment and detrainment of particles from the static bed and into the fully fluidized region, and vice versa. The resolution of the dynamics on the scale of this thin transition layer leads to an additional macroscopic boundary condition on local void fraction that must be applied at the interface of the fully fluidized region and the particle-bed surface. This closes the macroscopic two-phase flow problem in the fully fluidized region, with now the entrainment and detrainment mechanism rationally accounted for, and the need for phenomenological parameters removed. This closed model can now be used to formulate the problem in the fluidized region for any particular circumstances under consideration. To illuminate this theoretical model, we apply it in detail to a simple version of the helicopter cloud problem in the latter part of the paper. Although direct experimental comparisons for the results from our model are not possible at present without a specifically designed experimental programme, it is worthwhile referring to the high-speed photography in Rovere (2022), through the images labelled as figures 1.1, 1.3 and in particular 2.6, which gives a good qualitative comparison with our theoretical examples presented in § 7 of this paper.

An outline of the paper is as follows. In § 2, we give a detailed description of the physical model, and in particular our approach to rationally modelling the dynamics of particle entrainment and detrainment in an interfacial layer. This is followed in § 3 by the formulation of the physical model as a mathematical model, with a detailed treatment of the local dynamics in the interfacial layer, and the subsequent development of the corresponding macroscopic boundary conditions that close the full two-phase flow problem in the fully fluidized region. This is followed in § 4 by bringing together these components in the formulation of the full closed mathematical problem in the fluidized region. This mathematical problem is made dimensionless in § 5, where the formulation is reduced to the determination of the voidage field and a fluid velocity potential. The application to the helicopter cloud problem is developed in § 6, where the effect of the helicopter close to the surface is modelled by a suitably chosen fluid dipole and fluid line-vortex. The boundary value problems that are associated with this application are formulated in a form suitable for numerical solution via finite-difference approximation; in § 7 we present numerical examples demonstrating how varying certain parameters may change qualitative features of the voidage field, with precise details of our numerical method and its accuracy given in Appendix B.

2. The physical model

The physical model that we develop here is aimed at describing the particle and fluid flow fields generated above a ground surface composed of a packed bed of solid particles, when surface particles are raised into suspension via the flow of an incompressible fluid above the surface. The physical model is developed in three structured regions, as follows.

1. (1) The fluidized region in which the particles lifted from the ground surface are in suspension in the fluid. The flow in this region is a fully developed two-phase flow.

2. (2) The interfacial layer in which ground surface particles are transferred from the static bed of particles into the fluidized region and/or are returned from the fluidized region into the static bed of particles.

3. (3) The static particle-bed region where the particles are stored in the static particle bed.

A schematic diagram of the above structure is shown in figure 1. In the fluidized region, the two-phase flow varies on a length scale $l$ that is very much greater than the particle-spacing length scale $l_s$ and the particle radius $a$, so that

(2.1a,b)\begin{equation} a \ll l,\quad l_s \ll l. \end{equation}

Under (2.1a,b), the particle and fluid fields in the fluidized region can be modelled as two inter-penetrating continua, and the flow is governed by conservation of mass and momentum in the two interacting phases. The interfacial layer is a thin layer, of approximately uniform thickness, separating the static particle bed from the fully developed fluidized region above, and accommodates particle mass transfer to and from the static bed and the fluidized region. The thickness of the interfacial layer is typically of the order of ten particle diameters. Thus with $\delta$ being the interfacial layer thickness, we have

(2.2)\begin{equation} \delta \approx 20a, \end{equation}

so that formally, $a\ll \delta$, and the interfacial layer is thin compared to the two-phase flow continuum length scale $l$, but sufficiently wide to contain enough particles for local average particle velocities to be used. Thus we formally require

(2.3)\begin{equation} a\ll \delta \ll l. \end{equation}

The fluid flow in the thin interfacial layer is taken as being tangential to the surface of the static particle bed when viewed relative to the normal motion of the surface of the static particle bed. The particle motion within the interfacial layer, based on (2.3), is modelled as a force balance between the drag and lift on the particles induced by the surface tangential fluid flow in the interfacial layer, and the gravitational and particle collisional forces on the particles, together with the fluid shearing force at the upper surface of the interfacial layer and the particle–particle frictional force at the lower surface of the interfacial layer. Mass transfer between the static particle bed and the fluidized region, via the interfacial layer, then determines the surface deformation of the static particle bed. With $h$ representing the transverse static particle-bed surface deformation length scale (see figure 1), we consider the situation when

(2.4)\begin{equation} \delta \ll h \ll l. \end{equation}

Connection between the interfacial layer and the bulk fluidized region is made by requiring continuity of the normal particle and fluid mass flow fields and the voidage field at a well-defined interface. Of course, this interface is a model device to replace the rapid transition of the normal length scale from the interfacial layer into the bulk fluidized region above. Provided that such a transition is passive, its modelling as an interface is justified, and the theory bears this out. The static particle-bed region is taken as a particle store, with the packed particles in static equilibrium (without strong cohesion). For the specific problem addressed in § 6, namely the helicopter cloud problem, the continuum length scale $l$ will be taken as the geometric dimension of the helicopter blade length. The hydrodynamic effect of the helicopter is modelled as a fluid dipole, aligned with the axis of blade rotation and located at the helicopter blade rotation point, to represent the fluid downdraught, together with a line-vortex through the axis of rotation of the blades, to represent the fluid swirl.

Figure 1. Schematic diagram of the physical model.

3. The mathematical model

Based on the physical model described in § 2, we derive the dynamical equations of motion in the fluidized region and the interfacial layer. We begin in the fluidized region. Throughout, reference will be made to a fixed Cartesian coordinate system with coordinates $x$, $y$ and $z$, and with $z$ pointing vertically upwards. The usual associated unit vectors are $\boldsymbol {i}$, $\boldsymbol {j}$ and $\boldsymbol {k}$, respectively. The position vector is $\boldsymbol {r} = (x,y,z)$, and $\boldsymbol {r}_h = (x,y)$ corresponds to the position vector in the horizontal $xy$-plane. The location of the lower surface of the thin interfacial layer (which is also the upper surface of the static particle bed) is taken as being at

(3.1)\begin{equation} z = \xi (\boldsymbol{r}_h,t), \end{equation}

for $\boldsymbol {r}_h\in \mathbb {R}^2$ and $t\geq 0$, with $t$ being time. It follows from (3.1) that the unit normal vector field at the thin interfacial layer, pointing into the fluidized region, is given by

(3.2)\begin{equation} \hat{\boldsymbol{n}}(\boldsymbol{r}_h,t) = \frac{(-\boldsymbol{\nabla}_h\xi,1)}{(1+|\nabla_h\xi|^2)^{1/2}}, \end{equation}

for $\boldsymbol {r}_h\in \mathbb {R}^2$ and $t\geq 0$. Here,

(3.3)\begin{equation} \boldsymbol{\nabla}_h(\cdot) = \frac{\partial}{\partial x} (\cdot)\,\boldsymbol{i} + \frac{\partial}{\partial y}(\cdot)\,\boldsymbol{j}\end{equation}

is the usual horizontal gradient operator. We now consider the equations of motion in the fluidized region.

3.1. The fluidized region

In the fluidized region, the fluid and particle phases are modelled as inter-penetrating continua, with the fluid being incompressible and inviscid on the two-phase continuum length scale, but viscous on the particle length scale. Conservation of mass in the fluid and particle phases then requires, respectively,

(3.4)$$\begin{gather} E_t + \boldsymbol{\nabla}\boldsymbol{\cdot} (E\boldsymbol{u}) = 0, \end{gather}$$

(3.5)$$\begin{gather}-E_t + \boldsymbol{\nabla}\boldsymbol{\cdot} ((1-E)\boldsymbol{v}) = 0, \end{gather}$$

with $\boldsymbol {r}$ in $D(t)$ and $t > 0$. Here, $D(t)$ is the interior of the fluidized region occupied by the inter-penetrating continua at time $t\geq 0$, $E(\boldsymbol {r},t)$ is the voidage field (representing the volume of fluid per unit volume in $D(t)$ and taking values between the packing voidage of the particles $E_s$ (${<}1$) and unity), $\boldsymbol {u}(\boldsymbol {r},t)$ is the fluid velocity field, and $\boldsymbol {v}(\boldsymbol {r},t)$ is the particle velocity field. Subscript $t$ represents partial differentiation with respect to time, and $\boldsymbol {\nabla }$ is the usual gradient operator:

(3.6)\begin{equation} \boldsymbol{\nabla}(\cdot) = \frac{\partial}{\partial x}(\cdot)\,\boldsymbol{i} + \frac{\partial}{\partial y}(\cdot)\,\boldsymbol{j} + \frac{\partial}{\partial z}(\cdot)\,\boldsymbol{k}. \end{equation}

Conservation of momentum in the fluid and particle phases also requires, respectively,

(3.7)$$\begin{gather} \rho_fE(\boldsymbol{u}_t + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}) ={-} \boldsymbol{\nabla} p - \rho_f g E \boldsymbol{k} - \boldsymbol{D}, \end{gather}$$

(3.8)$$\begin{gather}\rho_s (1-E)(\boldsymbol{v}_t + (\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{v}) ={-}\boldsymbol{\nabla} p_s - \rho_s g (1-E) \boldsymbol{k} + \boldsymbol{D}, \end{gather}$$

with $\boldsymbol {r}$ in $D(t)$ and $t>0$. Here, $p=p(\boldsymbol {r},t)$ is the fluid pressure field, $p_s=p_s(E(\boldsymbol {r},t))$ is the particle–particle collisional pressure field, and $\boldsymbol {D}$ is the total drag force per unit volume exerted by the fluid on the particles (which we anticipate, at the single-particle scale in the fluidized flow field, will significantly dominate the single particle lift force). In addition, $\rho _f$ and $\rho _s$ are the constant fluid and particle material densities, respectively, and $g$ is the acceleration due to gravity. It remains to relate the particle collision pressure to the voidage field, and determine the form of the drag force. In general, the particle–particle collisional pressure is a decreasing function of voidage, which approaches zero as the voidage approaches unity. In the absence of detailed experimental evidence, we adopt, as a first approximation, the simple linear form

(3.9)\begin{equation} p_s=p_s(E) = p_0(1-E), \end{equation}

with $p_0 > 0$ being a material constant. For the drag force per unit volume, we write

(3.10)\begin{equation} \boldsymbol{D} = \beta(E)\,n_p (6{\rm \pi} a \mu (\boldsymbol{u}-\boldsymbol{v})), \end{equation}

where $n_p$ is the number of particles per unit volume, so that $n_p = (1-E)/((4/3){\rm \pi} a^3)$, whilst the bracketed term is the Stokes drag on a single particle, and $\beta =\beta (E)$ represents the effect of neighbouring particles on the Stokes drag. In general, $\beta (E)\geq 1$ is a decreasing function of $E\in [E_s,1]$, and $\beta (1)=1$. Thus we have

(3.11)\begin{equation} \boldsymbol{D} = \frac{9}{2}\,\frac{\mu}{a^2}\,(1-E)\,\beta(E)\,(\boldsymbol{u}-\boldsymbol{v}), \end{equation}

with $\mu$ being the kinematic viscosity of the fluid. On substituting from (3.9) and (3.11) into (3.7) and (3.8), we have

(3.12)$$\begin{gather} \rho_f E(\boldsymbol{u}_t + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}) ={-}\boldsymbol{\nabla} p - \rho_f g E \boldsymbol{k} -(9\mu/2a^2) (1-E)\,\beta(E)\,(\boldsymbol{u}-\boldsymbol{v}), \end{gather}$$

(3.13)$$\begin{gather}\rho_s (1-E)(\boldsymbol{v}_t + (\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{v}) = p_0\,\boldsymbol{\nabla} E - \rho_s g (1-E) \boldsymbol{k} + (9\mu/2a^2) (1-E)\,\beta(E)\,(\boldsymbol{u}-\boldsymbol{v}), \end{gather}$$

with $\boldsymbol {r}$ in $D(t)$ and $t > 0$. In obtaining (3.7) and (3.8) (and consequently (3.12) and (3.13)), for the micro-scale particle–fluid interaction, we have only included the drag term, which we expect to dominate over lift and buoyancy in the bulk fluidized region, where the fluid flow is primarily inviscid and irrotational on the two-phase continuum length scale (except in the interfacial layer, where ground effect vorticity will enhance lift; see § 3.2). Thus the dynamics in the fluidized region is represented by the four PDEs (3.4) and (3.5) with (3.12) and (3.13). We now consider the dynamics in the interfacial layer.

3.2. The interfacial layer

Within the thin interfacial layer, we consider the dynamics on the thickness scale $\delta$ to be in local equilibrium relative to the slow (compared to fluid velocity scale) normal motion of the interfacial layer, with the fluid velocity $\boldsymbol {u}_I$ and particle velocity $\boldsymbol {v}_I$ being tangential to the interfacial layer relative to the normal motion of the interfacial layer, and without variation across the interfacial layer. Normal to the interfacial layer, the local lift force balances the gravity force on the particles, whilst tangential to the interface, the drag and fluid shearing forces on the particles are balanced by the gravity force and the static particle-bed friction force on the particles. We consider a cylindrical cross-sectional element, of radius O $(\delta )$, passing through the interfacial layer, and centred at $\boldsymbol {r} = (\boldsymbol {r}_h,\xi (\boldsymbol {r}_h,t))$, as shown in figure 2.

Figure 2. Cylindrical element across the interfacial layer. The element is centred at $\boldsymbol {r}=(\boldsymbol {r}_h,\xi (\boldsymbol {r}_h,t))$ which is represented by $\bullet$. Here, $\varPi$ is the interfacial layer tangent plane at $\boldsymbol {r}=(\boldsymbol {r}_h,\xi (\boldsymbol {r}_h,t))$. The cross-sectional area of the cylindrical element is $A$.

A force balance on the particles in the cylindrical element in figure 2 gives, in equilibrium,

(3.14)\begin{equation} \boldsymbol{F}_g + \boldsymbol{F}_L + \boldsymbol{F}_D + \boldsymbol{F}_s + \boldsymbol{F}_f = 0, \end{equation}

where $\boldsymbol {F}_g$ is the gravitational force on the particles, $\boldsymbol {F}_L$ is the lift force on the particles, and $\boldsymbol {F}_D$ is the drag force on the particles, within the element. In addition, $\boldsymbol {F}_s$ is the fluid shearing force exerted on the particles over the upper surface of the element, whilst $\boldsymbol {F}_f$ is the frictional force exerted on the particles over the lower surface of the element. As we are considering the local dynamics in the interfacial layer, on the length scale $\delta$ ($\ll h \ll l$), to be in equilibrium, $\boldsymbol {u}_I$ and $\boldsymbol {v}_I$, together with the voidage in the interfacial layer $E_I\in [E_s,1]$, may be considered as fixed throughout the cylindrical element. We then have

(3.15)\begin{equation} \boldsymbol{F}_g ={-}(1-E_I) \rho_s g\delta A((\boldsymbol{k}\boldsymbol{\cdot}\hat{\boldsymbol{n}})\hat{\boldsymbol{n}} + (\boldsymbol{k} - (\boldsymbol{k}\boldsymbol{\cdot}\hat{\boldsymbol{n}})\hat{\boldsymbol{n}})), \end{equation}

noting that the vector $\boldsymbol {k} - (\boldsymbol {k}\boldsymbol {\cdot }\hat {\boldsymbol {n}})\hat {\boldsymbol {n}}$ lies in the tangent plane $\varPi$ to the interfacial layer. Next, we have

(3.16)\begin{equation} \boldsymbol{F}_L = n_p\delta A\varPhi(E_I)\,L_p\hat{\boldsymbol{n}}, \end{equation}

where $L_p$ is the lift force on a single particle in the cylindrical element, and $\varPhi =\varPhi (E_I)$ accounts for the effect of neighbouring particles. For a single particle in the interfacial layer, moving through the fluid as a rolling grain, we have

(3.17)\begin{equation} L_p = 4{\rm \pi} a^2 \rho_f\,| \boldsymbol{u}_I - \boldsymbol{v}_I|^2, \end{equation}

which is the Magnus force induced on an individual rolling particle grain as it moves along the interfacial layer. We anticipate that in this incipiently fluidized interfacial layer, this Magnus force dominates the Saffman lift force due to the local shear in the fluid flow. The function $\varPhi (E_I)$ is anticipated to be, in general, monotone increasing in $E_I$, with $\varPhi (E_s) <1$ and $\varPhi (1)=1$. Without detailed experimental evidence, we adopt the simple linear-form approximation

(3.18)\begin{equation} \varPhi(E_I) = E_I. \end{equation}

On substituting from (3.17) and (3.18) into (3.16), with $n_p = 3(1-E_I)/(4{\rm \pi} a^3)$, we obtain

(3.19)\begin{equation} \boldsymbol{F}_L = \left(\frac{3\rho_f}{a}\right)\delta AE_I(1-E_I)\,| \boldsymbol{u}_I - \boldsymbol{v}_I|^2\,\hat{\boldsymbol{n}}. \end{equation}

Similarly, adopting (3.11), with $E=E_I$, $\boldsymbol {u}=\boldsymbol {u}_I$ and $\boldsymbol {v}=\boldsymbol {v}_I$, we obtain

(3.20)\begin{equation} \boldsymbol{F}_D = \frac{9}{2}\,\frac{\mu}{a^2}\,\delta A(1-E_I)\,\beta(E_I)\,(\boldsymbol{u}_I - \boldsymbol{v}_I), \end{equation}

with $\beta (\cdot )$ as introduced in § 3.1. As discussed earlier, we recall that the vector $(\boldsymbol {u}_I - \boldsymbol {v}_I)$ is taken to lie in the tangent plane $\varPi$. For the upper surface shearing force, we may thus write

(3.21)\begin{equation} \boldsymbol{F}_s = A (1-E_I) \tau\,| \boldsymbol{u}_I - \boldsymbol{v}_I|\,(\boldsymbol{u}_I - \boldsymbol{v}_I),\end{equation}

where $\tau >0$ is the turbulent fluid shear coefficient. Finally, the particle–particle friction force at the lower surface is taken as

(3.22)\begin{equation} \boldsymbol{F}_f = A (1-E_I) F \boldsymbol{v}_I, \end{equation}

where $F > 0$ is the dynamic friction coefficient. We now substitute from (3.15), (3.19), (3.20), (3.21) and (3.22) into (3.14), and separate components normal and tangential to the tangent plane $\varPi$. We obtain the scalar equation

(3.23)\begin{equation} \varLambda E_I\,|\boldsymbol{u}_I - \boldsymbol{v}_I|^2 - (\boldsymbol{k}\boldsymbol{\cdot}\hat{\boldsymbol{n}}) = 0, \end{equation}

and the tangent plane vector equation

(3.24)\begin{equation} -\rho_s g \delta (\boldsymbol{k} - (\boldsymbol{k}\boldsymbol{\cdot}\hat{\boldsymbol{n}})\hat{\boldsymbol{n}}) + \frac{9}{2}\,\frac{\delta \mu}{a^2}\,\beta(E_I)\,(\boldsymbol{u}_I - \boldsymbol{v}_I) + \tau\,|\boldsymbol{u}_I - \boldsymbol{v}_I|\,(\boldsymbol{u}_I - \boldsymbol{v}_I) - F \boldsymbol{v}_I = \boldsymbol{0}_T, \end{equation}

where $\boldsymbol {0}_T$ is the zero vector in the tangent plane $\varPi$, and

(3.25)\begin{equation} \varLambda = \frac{3\rho_f}{a \rho_s g}. \end{equation}

Since the interfacial layer is thin ($\delta \ll h \ll l$), following (2.2), we conclude that the dominant terms on the left-hand side of (3.24) are those representing upper surface shear and lower surface friction (surface forces dominating body forces). With $u_I^s$ and $v_I^s$ being typical scales for $\boldsymbol {u}_I$ and $\boldsymbol {v}_I$, a balance between upper surface shear $\boldsymbol {F}_s$ and lower surface friction $\boldsymbol {F}_f$ gives

(3.26)\begin{equation} v_I^s = \textit{O}\left(\frac{\tau}{F}\,(u_I^s)^2\right) \end{equation}

and

(3.27)\begin{equation} \delta \ll \min \left(\frac{\tau(u_I^s)^2}{\rho_s g}, \frac{\tau a^2 u_I^s}{\mu} \right). \end{equation}

With this approximation, the dominant form of (3.24) is

(3.28)\begin{equation} \tau\,|\boldsymbol{u}_I - \boldsymbol{v}_I|\,(\boldsymbol{u}_I - \boldsymbol{v}_I) - F \boldsymbol{v}_I = \boldsymbol{0}_T. \end{equation}

Also, we anticipate that

(3.29)\begin{equation} v_I^s \ll u_I^s, \end{equation}

which follows from (3.26), when $\tau u_I^s \ll F$. Thus

(3.30)\begin{equation} |\boldsymbol{v}_I| \ll |\boldsymbol{u}_I|, \end{equation}

which allows us to approximate both of (3.23) and (3.28), which become

(3.31)$$\begin{gather} \varLambda E_I\,|\boldsymbol{u}_I|^2 - (\boldsymbol{k}\boldsymbol{\cdot}\hat{\boldsymbol{n}}) = 0, \end{gather}$$

(3.32)$$\begin{gather}\tau\,|\boldsymbol{u}_I|\,\boldsymbol{u}_I - F \boldsymbol{v}_I = \boldsymbol{0}_T. \end{gather}$$

A further simplification of (3.31) can be made in recalling from (2.4) that $h\ll l$, so the interfacial layer slope, which is O $(h/l)$, is small. In particular,

(3.33)\begin{equation} \hat{\boldsymbol{n}} = \boldsymbol{k} + \textit{O}(h/l), \end{equation}

so we may take

(3.34)\begin{equation} \hat{\boldsymbol{n}} \sim \boldsymbol{k}, \end{equation}

after which (3.31) becomes

(3.35)\begin{equation} \varLambda E_I\,|\boldsymbol{u}_I|^2 - 1 = 0. \end{equation}

Since $E_s \leq E_I \leq 1$, with $E_s$ being the packing voidage for the particles, it follows from (3.35) that the interfacial layer collapses for $|\boldsymbol {u}_I|^2 < 1/\varLambda$ (with $E_I = 1$), whilst it becomes saturated (with $E_I=E_s$) for $|\boldsymbol {u}_I|^2 > 1/(\varLambda E_s)$. Between these limiting values of $|\boldsymbol {u}_I|^2$, we have, from (3.32) and (3.35), that in the interfacial layer,

(3.36)$$\begin{gather} \boldsymbol{v}_I = (\tau/F)\,|\boldsymbol{u}_I|\,\boldsymbol{u}_I, \end{gather}$$

(3.37)$$\begin{gather}E_I = {\frac{1}{\varLambda\,|\boldsymbol{u}_I|^2}}. \end{gather}$$

In (3.37), the local voidage in the interfacial layer $E_I$ is related to the local tangential fluid velocity in the interfacial layer, $\boldsymbol {u}_I$, whilst (3.36) relates the local tangential particle velocity in the interfacial layer $\boldsymbol {v}_I$ to the local tangential fluid velocity in the interfacial layer $\boldsymbol {u}_I$, and the form is similar to a mobile bed-load sediment transport function (see, for example, Bagnold (1941) for the classical discussion of bed-load sediment transport formulae in a wind-blown sand environment). It should also be noted that throughout the above determination of the overall force balance on particles in the small continuum (containing many particles) element across the interfacial layer, we have taken the primary single-particle representative forces for drag, friction, shear and lift to be as simple as possible, whilst retaining their principal features, and this is sufficient for our present purpose. However, it is recognized that if required, these basic forms could be refined, and thus improved, by introducing higher-order contributions, using, for example, the recent results of Smith & Palmer (2019), Palmer & Smith (2020,  2021) and Jolley & Smith (2022), amongst others, arising from detailed studies of the motion of small free bodies located in the thin viscous boundary layer next to a wall, in high and moderate Reynolds number flows. We next consider the motion of the interfacial layer.

3.3. Interfacial layer motion

We consider a cylindrical volume from $z=-d$ in the static bed region ($z=-d$ forming a reference level in the static particle-bed region) and extending vertically upwards to the upper surface of the interfacial layer at $z=\xi + \delta$. The situation and nomenclature are illustrated in figure 3.

Figure 3. A qualitative schematic diagram representing a small element as a vertical slice through the static bed and the interfacial layer, with an illustration of the notation used in analysing this element.

Now the rate of change with respect to time $t$ of the total particle mass in the cylindrical element must be equal to the total mass flux of particles into or out of the cylindrical element. The total particle mass in the cylindrical element is $M$, where

(3.38)\begin{align} M &= \int\int_{A} \int_{z={-}d}^{z=\xi} \rho_s (1-E_s)\,\mathrm{d} z\,\mathrm{d} A + \int\int_{A} \int_{z=\xi}^{z=\xi+\delta} \rho_s (1-E_I)\,\mathrm{d} z\,\mathrm{d} A \nonumber\\ &= \rho_s(1-E_s) \int\int_{A} (\xi+d)\,\mathrm{d} A + \rho_s \delta \int\int_{A} (1-E_I)\,\mathrm{d} A, \end{align}

recalling that $E_s$ is the packing voidage of the particles in the static particle bed, and $E_I$ is independent of $z$ in the interfacial layer. The total mass flux of particles into the cylindrical element is $Q$, where

(3.39)\begin{align} Q &={-}\int\int_{S} \rho_s (1-E_I) \boldsymbol{v}_I \boldsymbol{\cdot} \hat{\boldsymbol{N}} \,\mathrm{d} s - \int\int_{A'} \rho_s (1-E_I)\,\boldsymbol{v}(\boldsymbol{r}_h,\xi+\delta,t)\boldsymbol{\cdot}\hat{\boldsymbol{n}} \, \mathrm{d} A' \nonumber\\ &={-}\int_{\mathcal{C}} \rho_s \delta (1-E_I) \boldsymbol{v}_I \boldsymbol{\cdot} \hat{\boldsymbol{N}} \, \mathrm{d} \ell - \int\int_{A'} \rho_s (1-E_I)\,\boldsymbol{v}(\boldsymbol{r}_h,\xi+\delta,t)\boldsymbol{\cdot}\hat{\boldsymbol{n}} \, \mathrm{d} A' \nonumber\\ &={-}\rho_s\delta\int\int_{A} \boldsymbol{\nabla}_h\boldsymbol{\cdot}((1-E_I)\boldsymbol{v}_I) \, \mathrm{d} A - \rho_s \int\int_{A'} (1-E_I)\,\boldsymbol{v}(\boldsymbol{r}_h,\xi+\delta,t)\boldsymbol{\cdot}\hat{\boldsymbol{n}} \, \mathrm{d} A', \end{align}

up to O $(h/l)$, via Green's theorem in the plane, with $\boldsymbol {\nabla }_h(\cdot )$ defined as in (3.3), and making use of (3.33). Finally, we may rewrite the surface integral over $A'$ as a surface integral over $A$ to obtain, on using (3.33),

(3.40)\begin{equation} Q = \int\int_{A} (- \rho_s\delta\,\boldsymbol{\nabla}_h\boldsymbol{\cdot}((1-E_I)\boldsymbol{v}_I) - \rho_s (1-E_I)\,\boldsymbol{v}(\boldsymbol{r}_h,\xi+\delta,t)\boldsymbol{\cdot}\hat{\boldsymbol{n}})\, \mathrm{d} A. \end{equation}

It now follows from the earlier conservation statement, via (3.38) and (3.40), that

(3.41)\begin{align} & (1-E_s)\xi_t - \delta {E_I}_t + \delta\,\boldsymbol{\nabla}_h \boldsymbol{\cdot}((1-E_I)\boldsymbol{v}_I) \nonumber\\ &\quad ={-} (1-E_I)\,\boldsymbol{v}(\boldsymbol{r}_h,\xi+\delta,t)\boldsymbol{\cdot}\hat{\boldsymbol{n}},\quad \boldsymbol{r}_h\in\mathbb{R}^2,\enspace t>0. \end{align}

A further simplification can be made here, based upon (2.4). The ratio of the second two terms on the left-hand side of (3.41) to the first term on the left-hand side of (3.41) is of O $(\delta /h)$, so (3.41) may be approximated by

(3.42)\begin{equation} (1-E_s)\xi_t ={-} (1-E_I)\,\boldsymbol{v}(\boldsymbol{r}_h,\xi,t)\boldsymbol{\cdot}\hat{\boldsymbol{n}},\quad \boldsymbol{r}_h\in\mathbb{R}^2,\enspace t>0, \end{equation}

which finally describes the dynamics of the interfacial layer motion.

3.4. Boundary conditions between the fluidized region and the interfacial layer

First, the normal fluid velocity in the fluidized region at $z=\xi +\delta$ must agree with the normal velocity of the interfacial layer. This requires

(3.43)\begin{equation} \xi_t = \boldsymbol{u}(\boldsymbol{r}_h,\xi,t)\boldsymbol{\cdot}\hat{\boldsymbol{n}},\quad \boldsymbol{r}_h\in\mathbb{R}^2,\enspace t>0, \end{equation}

where $\xi +\delta$ is approximated as $\xi$ via (2.4). Second, the tangential fluid velocity in the fluidized region at $z=\xi +\delta$ must agree with the tangential fluid velocity in the interfacial layer. This requires, after using (2.4), that

(3.44)\begin{equation} \boldsymbol{u}_I(\boldsymbol{r}_h,t) = \boldsymbol{u}(\boldsymbol{r}_h,\xi,t) - (\boldsymbol{u}(\boldsymbol{r}_h,\xi,t)\boldsymbol{\cdot}\hat{\boldsymbol{n}})\hat{\boldsymbol{n}} \equiv \boldsymbol{u}_h(\boldsymbol{r}_h,t),\quad \boldsymbol{r}_h\in\mathbb{R}^2,\enspace t>0. \end{equation}

Finally, the voidage field in the fluidized region at $z=\xi +\delta$ must agree with the voidage field in the interfacial layer, which, on using (2.4), requires

(3.45)\begin{equation} E_I(\boldsymbol{r}_h,t) = E(\boldsymbol{r}_h,\xi,t),\quad \boldsymbol{r}_h\in\mathbb{R}^2,\enspace t>0. \end{equation}

The three equations (3.43)–(3.45) provide boundary conditions relating conditions in the fluidized region to conditions in the interfacial layer.

4. Full problem in the fluidized region

Our main objective is to describe the dynamics in the fluidized region, together with the motion of the interfacial layer. To this end, we can, via elimination, obtain a decoupled problem in the fluidized region. Via (3.4) and (3.5) with (3.12) and (3.13), we have

(4.1)$$\begin{gather} E_t + \boldsymbol{\nabla} \boldsymbol{\cdot} (E\boldsymbol{u}) = 0, \end{gather}$$

(4.2)$$\begin{gather}-E_t + \boldsymbol{\nabla}\boldsymbol{\cdot} ((1-E)\boldsymbol{v}) = 0, \end{gather}$$

(4.3)$$\begin{gather}\rho_f E(\boldsymbol{u}_t + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}) ={-}\boldsymbol{\nabla} p - \rho_f g E \boldsymbol{k} - (9\mu/2a^2) (1-E)\,\beta(E)\,(\boldsymbol{u}-\boldsymbol{v}), \end{gather}$$

(4.4)$$\begin{gather}\rho_s (1-E)(\boldsymbol{v}_t + (\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{v}) = p_0\,\boldsymbol{\nabla} E - \rho_s g (1-E) \boldsymbol{k} + (9\mu/2a^2) (1-E)\,\beta(E)\,(\boldsymbol{u}-\boldsymbol{v}), \end{gather}$$

for $\boldsymbol {r}$ in $D(t)$, $t > 0$. Also, we must have

(4.5)\begin{equation} ((1-E_s)\boldsymbol{u} + (1-E) \boldsymbol{v})\boldsymbol{\cdot}\hat{\boldsymbol{n}} = 0\quad \mbox{on }z=\xi, \mbox{ with }\boldsymbol{r}_h\in\mathbb{R}^2, t>0, \end{equation}

via (3.42), (3.43) and (3.45), whilst (3.45) together with (3.37) and (3.44) requires

(4.6)\begin{equation} E = \left\{\begin{array}{ll} 1, & \mbox{for }|\boldsymbol{u}_h|^2<1/\varLambda, \\ \displaystyle{\dfrac{1}{\varLambda|\boldsymbol{u}_h|^2}}, & \mbox{for }1/\varLambda\leq|\boldsymbol{u}_h|^2\leq 1/(\varLambda E_s), \\ E_s, & \mbox{for }|\boldsymbol{u}_h|^2 > 1/(\varLambda E_s), \end{array}\right.\quad \mbox{on }z=\xi, \mbox{ with }\boldsymbol{r}_h\in\mathbb{R}^2, t>0, \end{equation}

where

(4.7)\begin{equation} \boldsymbol{u}_h(\boldsymbol{r}_h,t) = \boldsymbol{u}(\boldsymbol{r}_h,\xi,t) - (\boldsymbol{u}(\boldsymbol{r}_h,\xi,t)\boldsymbol{\cdot}\hat{\boldsymbol{n}})\hat{\boldsymbol{n}},\quad \boldsymbol{r}_h\in\mathbb{R}^2,\enspace t>0, \end{equation}

is the tangential component of the fluid velocity in the fluidized region at $z=\xi$. Finally, we have, via (3.43), that

(4.8)\begin{equation} \xi_t = \boldsymbol{u}\boldsymbol{\cdot}\hat{\boldsymbol{n}}\quad \mbox{on }z=\xi, \mbox{ with }\boldsymbol{r}_h\in\mathbb{R}^2, t>0. \end{equation}

Thus the decoupled problem in the fluidized region is composed of the PDEs (4.1)–(4.4) together with the boundary conditions (4.5)–(4.8) on $z=\xi$. In particular, we observe that the detailed consideration of the dynamics in the interfacial layer has led to the introduction of the key boundary condition (4.6) on the voidage at the interface. It is now convenient to write (4.1)–(4.8) in dimensionless form.

5. The dimensionless problem

Let $u_s$ be a typical velocity scale in the fluidized region (for both the fluid and particle phases), with the length scales $l$ and $h$ as introduced in § 2. The time scale in the fluidized region is then $l/u_s$, whilst a balance between pressure and drag in the momentum equation (4.3) leads to a pressure scale

(5.1)\begin{equation} p_s = \frac{\mu u_s l}{a^2}. \end{equation}

Based upon these scales, we introduce the dimensionless variables

(5.2a–h)\begin{equation} x=lx',\quad y=ly',\quad z=lz',\quad \xi=l\xi',\quad t=\frac{l}{u_s}t',\quad \boldsymbol{u}=u_s\boldsymbol{u}',\quad \boldsymbol{v}=u_s\boldsymbol{v}',\quad p=p_s p'. \end{equation}

On substitution from (5.2a–h) and (5.1) into (4.1)–(4.8), we obtain the dimensionless problem in the fluidized region as, after dropping primes for convenience,

(5.3)$$\begin{gather} E_t + \boldsymbol{\nabla} \boldsymbol{\cdot} (E\boldsymbol{u}) = 0, \end{gather}$$

(5.4)$$\begin{gather}-E_t + \boldsymbol{\nabla}\boldsymbol{\cdot} ((1-E)\boldsymbol{v}) = 0, \end{gather}$$

(5.5)$$\begin{gather}\bar{R}_fE(\boldsymbol{u}_t + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}) ={-}\boldsymbol{\nabla} p - G E \boldsymbol{k} - (1-E)\,\bar{\beta}(E)\,(\boldsymbol{u}-\boldsymbol{v}), \end{gather}$$

(5.6)$$\begin{gather}\rho \bar{R}_f (1-E)(\boldsymbol{v}_t + (\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{v}) = \alpha\,\boldsymbol{\nabla} E - \rho G (1-E) \boldsymbol{k} + (1-E)\,\bar{\beta}(E)\,(\boldsymbol{u}-\boldsymbol{v}), \end{gather}$$

for $\boldsymbol {r}$ in $D(t)$, $t > 0$, together with

(5.7)$$\begin{gather} ((1-E_s)\boldsymbol{u} + (1-E) \boldsymbol{v})\boldsymbol{\cdot}\hat{\boldsymbol{n}} = 0, \end{gather}$$

(5.8)$$\begin{gather}E = \displaystyle{\left\{ \begin{array}{ll} 1, & \mbox{for }|\boldsymbol{u}_h|^2<1/\gamma, \\ \displaystyle{\dfrac{1}{\gamma|\boldsymbol{u}_h|^2}}, & \mbox{for }1/\gamma\leq|\boldsymbol{u}_h|^2\leq 1/(\gamma E_s), \\ E_s, & \mbox{for }|\boldsymbol{u}_h|^2 > 1/(\gamma E_s), \end{array}\right.} \end{gather}$$

(5.9)$$\begin{gather}\xi_t = \boldsymbol{u} \boldsymbol{\cdot} \hat{\boldsymbol{n}}, \end{gather}$$

on $z=\xi$, with $\boldsymbol {r}_h\in \mathbb {R}^2$, $t>0$. Here, we have introduced the dimensionless parameters $\rho =\rho _s/\rho _f$, giving the density ratio of the particulate material to the fluid material, and

(5.10)\begin{equation} \bar{R}_f=\epsilon^2 R_f \end{equation}

where $\epsilon = a/l$ is the ratio of the microscopic scale particle radius and the macroscopic continuum length scale, and in general can be considered very small, whilst the fluid phase Reynolds number $R_f=\rho _f u_s l/\mu$, based on the macroscopic length scale $l$, is expected to be large, justifying our neglect of fluid phase macroscopic viscous terms in (5.5). Correspondingly, $\bar {R}_f=\epsilon ^2 R_f$ measures the ratio of the inertia terms to the drag terms in the momentum equations (5.5) and (5.6). In what follows, we will restrict attention to the situation when the macroscopic Reynolds number is large, but the ratio of the microscopic to macroscopic length scales is sufficiently small so that

(5.11)\begin{equation} 1 \ll R_f \ll \epsilon^{{-}2}, \end{equation}

hence $\epsilon ^2 \ll \bar {R}_f\ll 1$. As a consequence of this, in the fluid phase macroscopic momentum balance, inertia will play a secondary role. This is not uncommon in two-phase gas/solid particle flow, as seen in models of gas fluidized beds (see, for example, Needham & Merkin 1983). In addition, $G = \rho _f g a^2 / \mu u_s$ measures the ratio of the gravity terms to the drag terms in the momentum equations (5.5) and (5.6), whilst $\alpha =p_0 a^2 / l\mu u_s$ measures the ratio of the particle–particle pressure term to the drag term in the momentum equation (5.6). Finally, $\gamma =3u_s^2 \rho _f / a g \rho _s$ measures the ratio of the lift force to the gravity force acting on the particles in the interfacial layer. For convenience, we have also written

(5.12)\begin{equation} \bar{\beta}(E) = \tfrac{9}{2}\,\beta(E). \end{equation}

The PDEs (5.3)–(5.6), together with the boundary conditions (5.7)–(5.9), govern the dynamics in the fluidized region, and the dynamics of the location of the interfacial layer. For the problems that we wish to examine (in particular for the helicopter cloud problem, where the length scale $l$ is based upon the helicopter rotary circle radius, and the velocity scale $u_s$ is based upon the downdraught generated by the rotor, with the fluid being air and the particles being sand grains), we estimate that

(5.13a–c)\begin{equation} \rho \gg 1,\quad \bar{R}_f\ll \frac{\alpha}{\rho} \ll 1,\quad G \ll \frac{\alpha}{\rho} \ll 1. \end{equation}

Based upon these estimates of the dimensionless parameters for the cloud problem, we may approximate (5.5) and (5.6) by

(5.14)$$\begin{gather} \boldsymbol{0} ={-} \boldsymbol{\nabla} p - (1-E)\,\bar{\beta}(E)\,(\boldsymbol{u}-\boldsymbol{v}), \end{gather}$$

(5.15)$$\begin{gather}\boldsymbol{0} = \alpha\,\boldsymbol{\nabla} E + (1-E)\,\bar{\beta}(E)\,(\boldsymbol{u}-\boldsymbol{v}), \end{gather}$$

for $\boldsymbol {r}$ in $D(t)$, $t > 0$, which henceforth replace (5.5) and (5.6). In both phases, the fundamental balance in (5.14) and (5.15) is between drag and pressure gradient. We can now make some direct progress with (5.14) and (5.15). An addition of (5.14) and (5.15) gives

(5.16)\begin{equation} \boldsymbol{\nabla}(p-\alpha E) = \boldsymbol{0},\quad \boldsymbol{r}\in D(t),\enspace t > 0, \end{equation}

which gives, on integration,

(5.17)\begin{equation} p(\boldsymbol{r},t) = \alpha\,E(\boldsymbol{r},t) + F(t),\quad \boldsymbol{r}\in \bar{D}(t),\enspace t \geq 0, \end{equation}

where $\bar {D}(t)$ is the closure of $D(t)$, which replaces (5.14). Here, $F(t)$ is an arbitrary function of $t\geq 0$. We note from (5.17) that fluid isobars and isovoids coincide in $\bar {D}(t)$. A rearrangement of (5.15) then leads to

(5.18)\begin{equation} \boldsymbol{v}(\boldsymbol{r},t) = \boldsymbol{u}(\boldsymbol{r},t) + \frac{\alpha\,\boldsymbol{\nabla}(H(E(\boldsymbol{r},t)))}{1-E(\boldsymbol{r},t)},\quad \boldsymbol{r}\in \bar{D}(t),\enspace t \geq 0,\end{equation}

where for later convenience we have introduced

(5.19)\begin{equation} H(E) = \int_{E_s}^E \frac{\mathrm{d}\theta}{\bar{\beta}(\theta)},\quad E_s \leq E \leq 1, \end{equation}

as sketched in figure 4.

Figure 4. A typical graph of the function $H(E)$ against $E$.

We next add (5.3) and (5.4) to obtain

(5.20)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} (E\boldsymbol{u} + (1-E)\boldsymbol{v}) = 0,\quad \boldsymbol{r}\in D(t),\enspace t >0, \end{equation}

which replaces (5.4). We also observe from (5.18) that

(5.21)\begin{equation} \boldsymbol{\nabla} \times \boldsymbol{v} = \boldsymbol{\nabla} \times \boldsymbol{u},\quad \boldsymbol{r}\in D(t),\enspace t >0. \end{equation}

Attention is now restricted to irrotational flow in the fluid phase, so that $\boldsymbol {\nabla } \times \boldsymbol {u} = \boldsymbol {0}$ in $\bar {D}(t)$ for all $t\geq 0$. It then follows from (5.21) that $\boldsymbol {\nabla } \times \boldsymbol {v} = \boldsymbol {0}$ in $\bar {D}(t)$ for all $t\geq 0$, so the flow in the particle phase is also irrotational. Under these conditions, there exist scalar potentials $\phi$ and $\psi$ such that

(5.22a,b)\begin{equation} \boldsymbol{u} = \boldsymbol{\nabla} \phi,\quad \boldsymbol{v} = \boldsymbol{\nabla} \psi, \end{equation}

for all $\boldsymbol {r}$ in $\bar {D}(t)$, $t\geq 0$. It follows from (5.18) that

(5.23)\begin{equation} \boldsymbol{\nabla} \psi = \boldsymbol{\nabla}\phi + \frac{\alpha}{1-E}\,\boldsymbol{\nabla}(H(E)),\quad \boldsymbol{r}\in \bar{D}(t),\enspace t\geq 0, \end{equation}

after which substituting into the two remaining equations (5.20) and (5.3) gives

(5.24)$$\begin{gather} \nabla^2(\phi + \alpha\,H(E)) = 0, \end{gather}$$

(5.25)$$\begin{gather}E_t + \boldsymbol{\nabla} \phi \boldsymbol{\cdot}\boldsymbol{\nabla} E - \alpha E\,\nabla^2(H(E)) = 0, \end{gather}$$

for $\boldsymbol {r}$ in $D(t)$, $t> 0$, which are two coupled nonlinear PDEs determining $\phi$ and $E$. Then $\psi$ and $p$ are obtained from (5.23) and (5.17). Finally, we substitute from (5.22a,b) and (5.23) into the boundary conditions (5.7)–(5.9), which become

(5.26)$$\begin{gather} E = \displaystyle{\left\{\begin{array}{ll} 1, & \mbox{for }|\nabla_h \phi|^2<1/\gamma, \\ \displaystyle{\dfrac{1}{\gamma\,|\nabla_h \phi|^2}}, & \mbox{for }1/\gamma\leq|\nabla_h \phi|^2\leq 1/(\gamma E_s), \\ E_s, & \mbox{for }|\nabla_h \phi|^2 > 1/(\gamma E_s), \end{array}\right.} \end{gather}$$

(5.27)$$\begin{gather}(2-E_s-E)\,\boldsymbol{\nabla}\phi\boldsymbol{\cdot}\hat{\boldsymbol{n}}+ \alpha\,\boldsymbol{\nabla} H(E)\boldsymbol{\cdot}\hat{\boldsymbol{n}} = 0, \end{gather}$$

(5.28)$$\begin{gather}\xi_t = \boldsymbol{\nabla}\phi \boldsymbol{\cdot} \hat{\boldsymbol{n}}, \end{gather}$$

on $z=\xi$, with $\boldsymbol {r}_h$ in $\mathbb {R}^2$, $t>0$. We observe from (5.24) that $\phi +\alpha H(E)$ is harmonic in $D(t)$ for all $t > 0$, whilst $E$ satisfies a convection–diffusion equation (5.25) in $D(t)$ for all $t>0$, which is nonlinear, with diffusion coefficient

(5.29)\begin{equation} \alpha = \frac{p_0 a^2}{l \mu u_s}. \end{equation}

The boundary conditions (5.26)–(5.28) are coupled and nonlinear.

With a view to the helicopter cloud problem, in general we expect that the wind velocity scale will be sufficiently large so that

(5.30)\begin{equation} u_s \gg \frac{p_0 a^2}{l\mu}. \end{equation}

Therefore the parameter $\alpha$ may be regarded as small. It then follows that we may write

(5.31)\begin{equation} \phi = \bar{\phi} + \textit{O}(\alpha), \end{equation}

with $\bar {\phi }={\textit {O}}(1)$ as $\alpha \rightarrow 0$, whilst from (5.27) and (5.28) we may write

(5.32)\begin{equation} \xi = \alpha\bar{\xi} + \textit{O}(\alpha^2), \end{equation}

with $\bar {\xi }={\textit {O}}(1)$ as $\alpha \rightarrow 0$, from which it follows that $\hat {\boldsymbol {n}} = \boldsymbol {k} + {\textit {O}}(\alpha )$ as $\alpha \rightarrow 0$. On substituting from (5.31) and (5.32) into (5.24) and (5.27), the leading-order problem for $\bar {\phi }$ decouples, and is given by

(5.33)\begin{equation} \nabla^2 \bar{\phi} = 0,\quad \boldsymbol{r}\in D_0,\enspace t>0, \end{equation}

subject to

(5.34)\begin{equation} \boldsymbol{\nabla} \bar{\phi} \boldsymbol{\cdot} \boldsymbol{k} = 0,\quad z=0,\quad \boldsymbol{r}_h\in \mathbb{R}^2, \enspace t>0. \end{equation}

Here, $D_0$ is the fixed domain

(5.35)\begin{equation} D_0 = \mathbb{R}^2 \times (0,\infty). \end{equation}

With $E= {\textit {O}}(1)$ as $\alpha \rightarrow 0$, the remaining problem for the voidage $E$ is given, from (5.25) and (5.26), as

(5.36)\begin{equation} E_t + \boldsymbol{\nabla} \bar{\phi} \boldsymbol{\cdot} \boldsymbol{\nabla} E - \alpha E\,\nabla^2(H(E)) = 0,\quad \boldsymbol{r}\in D_0,\enspace t>0, \end{equation}

subject to

(5.37)\begin{equation} E = \left\{\begin{array}{ll} 1, & \mbox{for }|\nabla_h \bar{\phi}|^2<1/\gamma, \\ \displaystyle{\dfrac{1}{\gamma\,|\nabla_h \bar{\phi}|^2}}, & \mbox{for }1/\gamma\leq|\nabla_h \bar{\phi}|^2\leq 1/(\gamma E_s), \\ E_s, & \mbox{for }|\nabla_h \bar{\phi}|^2 > 1/(\gamma E_s), \end{array}\right.\quad z=0, \ \boldsymbol{r}_h\in \mathbb{R}^2, \ t>0. \end{equation}

It should be noted here that to maintain the spatial uniformity of the approximation in $E$ when $0 < \alpha \ll 1$, the dominant terms in both convection and diffusion have been retained in (5.36). Finally, on substitution from (5.32), (3.33) and (5.27) into (5.28), we obtain, at leading order,

(5.38)\begin{equation} \bar{\xi}_t ={-}\frac{1}{2-E_s-E}\,\boldsymbol{\nabla} H(E)\boldsymbol{\cdot}\boldsymbol{k},\quad z=0, \boldsymbol{r}_h\in \mathbb{R}^2, \ t>0. \end{equation}

We are now in a position to formulate the helicopter cloud problem in detail.

6. The helicopter cloud problem

In this section, we introduce an elementary model for the helicopter cloud problem, and examine this model in detail using the generic framework established in the preceding sections. Specifically, we model a helicopter hovering steadily, with rotor blades rotating in a horizontal plane, in otherwise still air above a sand bed, which has undisturbed level at $z=0$. The effect of the helicopter, with rotor blade length $l$, hovering steadily with the rotor at a vertical distance $lz_d$ above the undisturbed particle bed level at $z=0$, as shown in figure 5, is modelled as follows.

Figure 5. The parameter $z_d$ measures the ratio of the hovering height of the helicopter rotor to the rotor blade radius – though precise dimensions will vary, we might expect $z_d\approx 0.5$ to correspond to the helicopter resting on the sand bed.

1. (i) A half-space fluid dipole is placed at location

   (6.1)\begin{equation} \boldsymbol{r} = (0,0,lz_d), \end{equation}
   with $z_d>0$ (dimensionless). The axis of the dipole aligns with the unit vector $-\boldsymbol {k}$, and the moment of the dipole is $M_s>0$. This represents the downdraught effect of the helicopter rotor motion. The associated fluid velocity scale is then
   (6.2)\begin{equation} u_s = \frac{M_s}{l^3}. \end{equation}
   The dipole is placed at the centre point of the axis of a cylindrical shell aligned with the $z$-axis, which has axial length and radius $\varDelta l$, with $\varDelta \ll 1$. The cylindrical shell is permeable to fluid and particles, with the normal fluid and particle velocities taken as equal on the surface of the cylindrical shell (which represents drag domination close to the dipole, where the fluid velocity becomes unbounded approaching the dipole) and equal to the normal velocity field associated with the dipole at the centre of the shell. In dimensionless variables, the dipole is located at $\boldsymbol {r}=(0,0,z_d)$, and the cylindrical shell length and radius is $\varDelta$.

2. (ii) A fluid line-vortex is placed with its core along the positive $z$-axis. The strength of the line vortex is taken as $\varGamma _s>0$. This represents the swirl generated by the helicopter rotor motion.

We might expect this simplistic model of the flow field to at least represent the overall qualitative features of that flow field generated by a helicopter in hovering mode, at least away from the immediate neighbourhood of the helicopter body. The key objective is to provide a relevant case to illustrate an application of the generic theory developed in the preceding sections. Representing the interior of the cylindrical shell containing the dipole as $C$, with boundary $\partial C$, the fluidized region now occupies the domain $D_c = D_0 \setminus \bar {C}$. The modelling structures (i) and (ii) require us to introduce the boundary conditions, in dimensionless form,

(6.3)$$\begin{gather} \boldsymbol{\nabla}\bar{\phi} \boldsymbol{\cdot} \hat{\boldsymbol{N}} = \boldsymbol{\nabla} \phi_d \boldsymbol{\cdot} \hat{\boldsymbol{N}},\quad \boldsymbol{r} \in \partial C, \end{gather}$$

(6.4)$$\begin{gather}\boldsymbol{\nabla}\bar{\phi} = \frac{\varOmega}{2{\rm \pi} R}\,\hat{\boldsymbol{\theta}} + \textit{O}(1) \quad \mbox{as } R \rightarrow 0 \mbox{ uniformly in }\bar{D}_c. \end{gather}$$

Here, $(R,\theta,z)$ are the usual cylindrical polar coordinates in relation to the Cartesian coordinates $(x,y,z)$ that were introduced earlier, with unit vectors $\hat {\boldsymbol {R}}$, $\hat {\boldsymbol {\theta }}$ and $\boldsymbol {k}$. In (6.3), $\phi _d$ represents the harmonic half-space dipole, given by

(6.5)\begin{equation} \phi_d = \sum_{i=1}^2 \frac{({-}1)^{i-1} (z-z_i)}{4{\rm \pi}(R^2 + (z-z_i)^2)^{3/2}} \quad \mbox{in }\bar{D}_c, \end{equation}

with

(6.6a,b)\begin{equation} z_1 = z_d,\quad z_2 ={-}z_d, \end{equation}

whilst the dimensionless parameter $\varOmega$ in (6.4) is given by

(6.7)\begin{equation} \varOmega = \frac{\varGamma_s l^2}{M_s}. \end{equation}

Here, $\varOmega$ measures the ratio of the fluid velocity scale induced by the line-vortex to that induced by the dipole. In addition to (6.3) and (6.4), it follows from (i) and (5.22a,b) with (5.23) that we must have

(6.8)\begin{equation} \boldsymbol{\nabla} E \boldsymbol{\cdot} \hat{\boldsymbol{N}} = 0,\quad \boldsymbol{r}\in\partial C. \end{equation}

Throughout the above, $\hat {\boldsymbol {N}}$ is the unit outward normal on the cylindrical shell. It remains to consider the far field boundary conditions at large distances from the cylindrical shell surface $\partial C$. We require

(6.9)\begin{align} E \rightarrow 1\quad \mbox{as } (R^2 + z^2) \rightarrow \infty \mbox{ uniformly in }\bar{D}_c, \end{align}

(6.10)\begin{align} \boldsymbol{\nabla}\bar{\phi} \rightarrow 0\quad \mbox{as } (R^2 + z^2) \rightarrow \infty \mbox{ uniformly in }\bar{D}_c, \\ \bar{\xi} \rightarrow 0\quad \mbox{as } R\rightarrow \infty \mbox{ uniformly in } 0 \leq \theta < 2{\rm \pi}. \nonumber \end{align}

To conclude, we note that the dynamics of the helicopter rotor will generally relate the dipole strength $M_s$ and the line-vortex strength $\varGamma _s$ in the model. The simplest form, which we adopt here, is a linear relation

(6.11)\begin{equation} \varGamma_s = k M_s, \end{equation}

with $k>0$ a constant of proportionality, depending upon the particular helicopter rotor under consideration. The dimensionless parameter $\varOmega$ is then given by

(6.12)\begin{equation} \varOmega = k l^2, \end{equation}

independent of the dipole and line-vortex strengths. In its final form, the steady helicopter cloud model has reduced to the following problem for the fluid velocity potential $\bar {\phi }$, using (5.33), (5.34), (6.3), (6.4) and (6.10):

(6.13)\begin{equation} \left.\begin{gathered} \nabla^2 \bar{\phi} = 0,\quad \boldsymbol{r}\in D_c, \\ \boldsymbol{\nabla} \bar{\phi} \boldsymbol{\cdot} \boldsymbol{k} = 0, \quad z=0, \ \boldsymbol{r}_h\in \mathbb{R}^2, \\ \boldsymbol{\nabla}\bar{\phi} \boldsymbol{\cdot} \hat{\boldsymbol{N}} = \boldsymbol{\nabla} \phi_d \boldsymbol{\cdot} \hat{\boldsymbol{N}}, \quad \boldsymbol{r} \in \partial C,\\ \boldsymbol{\nabla}\bar{\phi} = \frac{\varOmega}{2{\rm \pi} R}\,\hat{\boldsymbol{\theta}} + \textit{O}(1) \quad \mbox{as } R \rightarrow 0 \mbox{ uniformly in }\bar{D}_c, \\ \boldsymbol{\nabla}\bar{\phi} \rightarrow 0\quad \mbox{as } |\boldsymbol{r}| \rightarrow \infty \mbox{ uniformly in }\bar{D}_c. \end{gathered}\right\} \end{equation}

The unique solution to this problem is readily established to be $\bar {\phi }: \bar {D}_c \rightarrow \mathbb {R}$, given by

(6.14)\begin{equation} \bar{\phi}(R,\theta,z) = \phi_d(R,z) + \frac{\varOmega}{2{\rm \pi}}\,\theta, \end{equation}

for $\boldsymbol {r}(R,\theta,z)\in \bar {D}_c$. The steady problem for $E: \bar {D}_c\rightarrow \mathbb {R}$ is then, via (5.36), (5.37), (6.8) and (6.9), given by

(6.15)\begin{equation} \boldsymbol{\nabla} \bar{\phi} \boldsymbol{\cdot}\boldsymbol{\nabla} E - \alpha E\,\nabla^2(H(E)) = 0,\quad \boldsymbol{r}\in D_c, \end{equation}

with

(6.16)\begin{equation} E = \left\{\begin{array}{ll} 1, & \mbox{for }|\nabla_h \bar{\phi}|^2<1/\gamma, \\ \displaystyle{\dfrac{1}{\gamma\,|\nabla_h \bar{\phi}|^2}}, & \mbox{for }1/\gamma\leq|\nabla_h \bar{\phi}|^2\leq 1/(\gamma E_s), \\ E_s, & \mbox{for }|\nabla_h \bar{\phi}|^2 > 1/(\gamma E_s), \end{array}\right.\quad z=0, \ \boldsymbol{r}_h\in \mathbb{R}^2, \ t>0, \end{equation}

and

(6.17)$$\begin{gather} \boldsymbol{\nabla} E \boldsymbol{\cdot} \hat{\boldsymbol{N}} = 0, \quad \boldsymbol{r}\in\partial C, \end{gather}$$

(6.18)$$\begin{gather}E \rightarrow 1\quad \mbox{as } |\boldsymbol{r}| \rightarrow \infty \mbox{ uniformly in }\bar{D}_c. \end{gather}$$

We now observe that $\bar {\phi }:\bar {D}_c\rightarrow \mathbb {R}$, as given in (6.14) and (6.5), is such that

(6.19)\begin{equation} \boldsymbol{\nabla} \bar{\phi}(R,z) = a(R,z)\,\hat{\boldsymbol{R}} + \frac{\varOmega}{2{\rm \pi} R}\,\hat{\boldsymbol{\theta}} + b(R,z)\,\boldsymbol{k} \end{equation}

for $\boldsymbol {r}\in \bar {D}_c$, with

(6.20)\begin{align} a(R,z) &= \frac{3R}{4{\rm \pi}} \left[\frac{z+z_d}{(R^2 + (z + z_d)^2)^{5/2}} - \frac{z-z_d}{(R^2 + (z - z_d)^2)^{5/2}} \right], \end{align}

(6.21)\begin{align} b(R,z) &= \frac{3}{4{\rm \pi}} \left[\frac{(z+z_d)^2}{(R^2 + (z + z_d)^2)^{5/2}} - \frac{(z-z_d)^2}{(R^2 + (z - z_d)^2)^{5/2}} \right] \nonumber\\ &\quad + \frac{1}{4{\rm \pi}} \left[\frac{1}{(R^2 + (z - z_d)^2)^{3/2}} - \frac{1}{(R^2 + (z + z_d)^2)^{3/2}} \right]. \end{align}

Thus via (6.19)–(6.21), we have

(6.22)\begin{equation} |\nabla_h \bar{\phi}(R,0)|^2 = \frac{9z_d^2 R^2}{4{\rm \pi}^2(R^2+z_d^2)^5} + \frac{\varOmega^2}{4{\rm \pi}^2 R^2}, \end{equation}

for $R> 0$. Using (6.22), the problem (6.15)–(6.18) may be rewritten as

(6.23)\begin{equation} \left.\begin{gathered} \boldsymbol{\nabla} \bar{\phi} \boldsymbol{\cdot} \boldsymbol{\nabla} E - \alpha E\,\nabla^2(H(E)) = 0, \quad \boldsymbol{r}\in D_c, \\ E = g (|\boldsymbol{r}_h|)\quad \mbox{on } z=0,\ \boldsymbol{r}_h\in\mathbb{R}^2, \\ \boldsymbol{\nabla} E \boldsymbol{\cdot} \hat{\boldsymbol{N}} = 0\quad \mbox{on } \boldsymbol{r}\in\partial C, \\ E \rightarrow 1\quad \mbox{as } |\boldsymbol{r}| \rightarrow \infty \mbox{ uniformly in }\bar{D}_c. \end{gathered}\right\} \end{equation}

which we henceforth refer to as [BVP]. The function $g:[0,\infty )\rightarrow \mathbb {R}$ is continuous and piecewise smooth, and is given by, via (6.16) and (6.22),

(6.24)\begin{equation} g(X) = \left\{\begin{array}{ll} E_s, & \mbox{for }|\nabla_h \bar{\phi}(X,0)|^2 > (\gamma E_s)^{{-}1}, \\ \displaystyle{\gamma^{{-}1}\left( \dfrac{9z_d^2 X^2}{4{\rm \pi}^2(X^2+z_d^2)^5} + \dfrac{\varOmega^2}{4{\rm \pi}^2 X^2} \right)^{{-}1}}, & \mbox{for }\gamma^{{-}1} \leq|\nabla_h \bar{\phi}(X,0)|^2\leq (\gamma E_s)^{{-}1}, \\ 1, & \mbox{for }|\nabla_h \bar{\phi}(X,0)|^2<\gamma^{{-}1}. \\ \end{array}\right. \end{equation}

We note that it is straightforward to refine the far field boundary condition (6.23) in [BVP] to

(6.25)\begin{equation} E(\boldsymbol{r}) = 1 + \textit{O}(|\boldsymbol{r}|^{{-}2})\quad \mbox{as } |\boldsymbol{r}| \rightarrow \infty \mbox{ uniformly in }\bar{D}_c. \end{equation}

On observing that

(6.26)\begin{equation} E_s \leq g(X) \leq 1 \quad \mbox{for all } X\in[0,\infty), \end{equation}

an application of the strong elliptic maximum principle (see, for example, Gilbarg & Trudinger 1998, chapter 9) establishes that any solution $E:\bar {D}_c\rightarrow \mathbb {R}$ to [BVP] must satisfy the inequality

(6.27)\begin{equation} E_s < E(\boldsymbol{r}) < 1\quad \mbox{for all } \boldsymbol{r}\in D_c. \end{equation}

It then follows from (6.27) and (6.25) that [BVP] has a classical unique solution (see, for example, Gilbarg & Trudinger 1998, chapter 9). An immediate consequence of uniqueness for [BVP], together with the rotational symmetry of [BVP] about the $z$-axis, is that the solution to [BVP] is axisymmetric in the cylindrical polar coordinates $(R,\theta,z)$; that is, if $E:\bar {D}_c\rightarrow \mathbb {R}$ is the solution to [BVP], then $E=E(R,z)$.

A final simplification can be made to [BVP]. In general $\bar {\beta }:[E_s,1]\rightarrow \mathbb {R}$ is monotone decreasing, with $\bar {\beta }(1)=9/2$. In the absence of further information, we will take, as a first approximation, the simple linear form

(6.28)\begin{equation} \bar{\beta}(E) = \bar{\beta}_0 - (\bar{\beta}_0 - 9/2) E \end{equation}

for $E\in [E_s,1]$, with $\bar {\beta }_0 > 9/2$ a material constant. It then follows, after an integration, that

(6.29)\begin{equation} H(E) = (\bar{\beta}_0 - 9/2)^{{-}1} \log(\bar{\beta}_0(\bar{\beta}_0 - (\bar{\beta}_0 - 9/2) (E-E_s))^{{-}1}) \end{equation}

for $E\in [E_s,1]$. We anticipate that, in general, the variation in $\bar {\beta }(E)$ over $E\in [E_s,1]$ will be small, so we write

(6.30)\begin{equation} \bar{\beta}_0 = \tfrac{9}{2} + \delta \beta_0, \end{equation}

with $0<\delta \beta _0\ll 1$, after which (6.29) may be approximated as

(6.31)\begin{equation} H(E) = \tfrac{2}{9} (E-E_s) + \textit{O}(\delta\beta_0) \end{equation}

as $\delta \beta _0\rightarrow 0$ uniformly for $E\in [E_s,1]$. The complete form of [BVP] when written in terms of the cylindrical polar coordinates $(R,\theta,z)$ is presented in Appendix A, and its domain $\omega$ is illustrated in figure 6.

Figure 6. The domain $\omega$ in the $(R,z)$ plane.

Before proceeding further with [BVP], we must consider the detailed form of $g:[0,\infty )\to \mathbb {R}$, as defined in (6.24). First, we examine $|\nabla _h \bar {\phi }(X,0)|^2$ for $X\in (0,\infty )$, the form of which depends upon the two parameters $\varOmega$ and $z_d$. We begin by observing that

(6.32)\begin{equation} |\nabla_h \bar{\phi}(X,0)|^2 \sim \frac{\varOmega^2}{4{\rm \pi}^2 X^2} \end{equation}

as $X\to 0^+$ and as $X\to \infty$. Recalling (6.22), we have

(6.33)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} X}|\nabla_h \bar{\phi}(X,0)|^2 = \frac{9 X^4 z_d^2(z_d^2- 4X^2)-\varOmega^2(X^2+z_d^2)^6}{2{\rm \pi}^2 X^3 (X^2+z_d^2)^6}, \end{equation}

hence $|\nabla _h \bar {\phi }(X,0)|^2$ is strictly monotone decreasing for $X\geq z_d/2$. For fixed $z_d>0$, we see from (6.32) and (6.33) that there is a critical value $\varOmega =\varOmega _c(z_d)$ such that $|\nabla _h \bar {\phi }(X,0)|^2$ will also be strictly monotone decreasing for $X\in [0,z_d/2]$, unless $\varOmega <\varOmega _c(z_d)$, where

(6.34)\begin{equation} \varOmega_c(z_d) = \max_{X\in[0,z_d/2]} f(X), \end{equation}

with

(6.35)\begin{equation} f(X):=\frac{3 X^2 z_d(z_d^2- 4X^2)^{1/2}}{(X^2+z_d^2)^3}, \quad X\in[0,z_d/2]. \end{equation}

We see immediately that $f(X)>0$ for $X\in (0,z_d/2)$, with $f(0)=f(z_d/2)=0$, and it is straightforward to calculate

(6.36)\begin{equation} f'(X) = \frac{36 X z_d}{(X^2+z_d^2)^4(z_d^2-4X^2)^{1/2}} \left(\left(\frac{2z_d^2}{3}-X^2\right)^2-\frac{5z_d^4}{18}\right), \end{equation}

hence $f'(X)=0$ if and only if $X=0$ (a local minimum) or

(6.37)\begin{equation} X= X_c(z_d) := \sqrt{\frac{4-\sqrt{10}}{6}}\,z_d \approx 0.373658\,z_d, \end{equation}

a local maximum, with $f'(X)>0$ for $X\in (0,X_c(z_d))$, and $f'(X)<0$ for $X\in (X_c(z_d),z_d/2)$. At $X=X_c(z_d)$, substitution into (6.35) and a little algebraic manipulation reveals that

(6.38)\begin{equation} \varOmega_c(z_d) = f(X_c(z_d)) = \frac{2}{5z_d^2}\,\sqrt{\frac{-505+188\sqrt{10}}{405}} \approx \frac{0.188046}{z_d^2}.\end{equation}

We summarize the cases $\varOmega \geq \varOmega _c(z_d)$ and $\varOmega <\varOmega _c(z_d)$ as follows.

1. (a) For ${\varOmega \geq \varOmega _c(z_d)}$, $|\nabla _h \bar {\phi }(X,0)|^2$ is strictly monotone decreasing in $X>0$, and for $\varOmega >\varOmega _c(z_d)$, it has the qualitative form shown in figure 7(a). If $\varOmega =\varOmega _c(z_d)$, then there is a point of inflection at $X=X_c$.

2. (b) For ${0< \varOmega < \varOmega _c(z_d)}$, $|\nabla _h \bar {\phi }(X,0)|^2$ has a local minimum and a local maximum in $X>0$, and has the qualitative form shown in figure 7(b).

Figure 7. A qualitative sketch of $|\nabla _h \bar {\phi }(X,0)|^2$ against $X$ for (a) $\varOmega > \varOmega _c(z_d)$, and (b) $0< \varOmega < \varOmega _c(z_d)$.

We can now construct $g:[0,\infty )\to \mathbb {R}$, via (6.24). We have

(6.39)\begin{equation} g(X) = \left\{\begin{array}{ll} E_s, & X\in I_1, \\ \displaystyle{\gamma^{{-}1}\left(\dfrac{9z_d^2 X^2}{4{\rm \pi}^2(X^2+z_d^2)^5} + \dfrac{\varOmega^2}{4{\rm \pi}^2 X^2} \right)^{{-}1}}, & X\in I_2, \\ 1, & X\in I_3. \\ \end{array}\right. \end{equation}

In case (a),

(6.40a–c)\begin{equation} I_1 = [0,R_1^-],\quad I_2 = (R_1^-,R_1^+), \quad I_3 = [R_1^+,\infty). \end{equation}

However, in case (b), depending upon the choice of $z_d$ and $\varOmega$, there are four possibilities, as follows:

1. (i) $I_1 = [0,R_1^-]$, $I_2 = (R_1^-,R_1^+)$, $I_3 = [R_1^+,\infty )$;

2. (ii) $I_1 = [0,R_1^-]$, $I_2 = (R_1^-,R_1^+)\cup (R_2^+,R_3^+)$, $I_3 = [R_1^+,R_2^+]\cup [R_3^+,\infty )$;

3. (iii) $I_1 = [0,R_1^-]\cup [R_2^-,R_3^-]$, $I_2 = (R_1^-,R_2^-)\cup (R_3^-,R_1^+)$, $I_3 = [R_1^+,\infty )$;

4. (iv) $I_1 = [0,R_1^-]\cup [R_2^-,R_3^-]$, $I_2 = (R_1^-,R_1^+)\cup (R_2^+,R_2^-) \cup (R_3^-,R_3^+)$, $I_3 = [R_1^+,R_2^+]\cup [R_3^+,\infty )$.

Here, $R_i^-$, $i=1,2,3$, are the consecutive positive zeros of the algebraic equation

(6.41)\begin{equation} |\nabla_h \bar{\phi}(X,0)|^2 = (\gamma E_s)^{{-}1}, \end{equation}

whilst $R_i^+$, $i=1,2,3$, are the consecutive positive zeros of the corresponding algebraic equation

(6.42)\begin{equation} |\nabla_h \bar{\phi}(X,0)|^2 = \gamma^{{-}1}. \end{equation}

We note that in case (b), with $z_d$, $\varOmega$ and $E_s$ fixed, small values of $\gamma$ correspond to case (i), whilst increasing $\gamma$ moves through cases (ii)–(iv).

In order to understand more precisely for which choices of $\gamma$, $E_s$, $z_d$ and $\varOmega$ we will see each of cases (i)–(iv), we now investigate a little further how the local maximum and minimum of $|\nabla _h \bar {\phi }(X,0)|^2$ in case (b) depend on $z_d$ and $\varOmega$. Defining $x_{max}$ and $x_{min}$, respectively, as the values of $X$ at which the local maximum and minimum of $|\nabla _h \bar {\phi }(X,0)|^2$ are achieved, we see from (6.33), with the change of variables

(6.43a,b)\begin{equation} \bar{x}=\left(\frac{X}{z_d}\right)^2, \quad \bar{\varOmega}=\varOmega^2 z_d^4, \end{equation}

that $\bar {x}_{max}:=(x_{max}/z_d)^2$ and $\bar {x}_{min}:=(x_{min}/z_d)^2$ are the solutions $\bar {x}\in (0,1/4)$ of

(6.44)\begin{equation} 9\bar{x}^2(1-4\bar{x})=\bar{\varOmega}(1+\bar{x}^6), \end{equation}

where if

(6.45)\begin{equation} \bar{\varOmega}> \bar{\varOmega}_c:=\varOmega_c^2 z_d^4 =\frac{4}{25} \left(\frac{-505+188\sqrt{10}}{405}\right) \approx 0.0353613, \end{equation}

then (6.44) has no solutions, if $\bar {\varOmega }=\bar {\varOmega }_c$, then (6.44) has one solution at $\bar {x}_*=(4-\sqrt {10})/ 6\approx 0.139620$, and if $0<\bar {\varOmega }<\bar {\varOmega }_c$, then (6.44) has precisely two solutions, $\bar {x}_{max}$ and $\bar {x}_{min}$. As $\bar {\varOmega }\to 0$, we see that $\bar {x}_{min}\to 0$ and $\bar {x}_{max}\to 1/4$. A little algebraic manipulation then reveals that

(6.46)\begin{equation} |\nabla_h \bar{\phi}(X,0)|^2 = \frac{9\bar{x}(2-3\bar{x})}{4{\rm \pi}^2(1+\bar{x})^6z_d^6}\quad \mbox{for }X=x_{max},x_{min}, \end{equation}

where $\bar {x}$ solves (6.44). From this, we see that

(6.47a,b)\begin{equation} |\nabla_h \bar{\phi}(x_{min},0)|^2 \to 0,\quad |\nabla_h \bar{\phi}(x_{max},0)|^2 \to \frac{576}{3125{\rm \pi}^2z_d^6}\approx \frac{0.0186755}{z_d^6}, \quad \mbox{as }\varOmega\to 0, \end{equation}

and

(6.48)\begin{equation} |\nabla_h \bar{\phi}(x_{min},0)|^2 = |\nabla_h \bar{\phi}(x_{max},0)|^2 = \frac{955+424\sqrt{10}}{10125{\rm \pi}^2z_d^6}\approx \frac{0.0229742}{z_d^6}\quad \mbox{for }\varOmega=\varOmega_c(z_d), \end{equation}

hence the range of possible values of $|\nabla _h \bar {\phi }(x_{min},0)|^2$ and $|\nabla _h \bar {\phi }(x_{max},0)|^2$ is rather narrow.

We now consider the parameters in [BVP]. There are six dimensionless parameters, namely, $E_s$, $\alpha$, $\gamma$, $\varOmega$, $z_d$ and $\varDelta$. We can give order of magnitude estimates of these parameters for the helicopter cloud problem using typical helicopter parameters as recorded, for example, in the papers by Wachspress et al. (2008), Govindarajan et al. (2013), Porcù et al. (2020), Tan et al. (2021,  2022a,b), Li et al. (2023) and Lin et al. (2024), together with the standard properties of air and sand, to approximate $l,M_s,\gamma _s$. The parameter $\varDelta$ measures the ratio of the helicopter rotor core radius to the rotor blade radius, and typically, $\varDelta \approx 10^{-2}$. The parameter $z_d$ measures the ratio of the hovering height of the helicopter rotor to the rotor blade radius (recall figure 5), and typically, $z_d\approx 10^{-1}$–$10^1$. For the fluid, the ratio of the hovering helicopter induced swirl velocity to the induced downwash velocity is measured by $\varOmega$, and using typical values extracted from the above papers, we have $\varOmega \approx 10^{-1}$–$10^0$. The parameter $E_s$ is the particle packing voidage, and for sand in air this is typically of magnitude $10^{-2}$. The ratio of the lift force per unit volume to the gravity force per unit volume on particles in the interfacial layer is given by $\gamma$, and again using the typical values extracted from the above papers, we have $\gamma \approx 10^{2}$–$10^3$. Finally, $\alpha$ measures the ratio of particle collisional pressure to drag in the particle flow phase, and we tentatively estimate, for sand fluidized in air, that $\alpha \approx 10^{-1}$. With these estimates of all dimensionless parameters appearing in [BVP], in the next section we consider the numerical solution to [BVP].

7. Numerical examples

The accurate and efficient numerical solution of [BVP] requires a degree of non-trivial consideration, and as such, we present the technical details and assessment in Appendix B, where we also establish convergence of our numerical scheme (see example B.1). Here, we now consider how the solution behaviour depends on our parameter choices. In examples 7.1–7.3, we fix

(7.1a–c)\begin{equation} E_s = 0.01,\quad \alpha = 0.1,\quad \varDelta = 0.04, \end{equation}

and consider the nine possible combinations of

(7.2a,b)\begin{equation} \gamma = 1000, 500, 100,\quad \varOmega = 1, 0.5, 0.1 \end{equation}

(with $\gamma =500$, $\varOmega =0.5$, corresponding to example B.1). As described in Appendix B, the accuracy of our numerical scheme depends on a number of discretization parameters. From our assessment there of the scheme, we here choose the parameters $N=8$, $L_R=32N$, $L_z=16N$ (compared to $L_R=12N$ and $L_z=6N$ in example B.1), and $\epsilon =1/(10N^2)$, in order to ensure that the computational domain is large enough that all relevant solution behaviour will be captured across this wide range of values of $\gamma$ and $\varOmega$. These parameter choices lead to the total number of degrees of freedom in our numerical solution being $\mbox {DOF}=2\,317\,056$ for $z_d = 0.5$, $\mbox {DOF}=2\,481\,536$ for $z_d = 0.9$, and $\mbox {DOF}=2\,646\,016$ for $z_d = 1.3$, for each of examples 7.1–7.3. In each of the figures below, we truncate the computational domains where appropriate for presentational purposes, noting that the solution satisfies $E\approx 1$ everywhere outside the plotted range – note that all plots are over the same range within each figure, but that the plotted range varies for examples 7.1–7.3 ($R\in [0,6]$, $z\in [0,5]$ for example 7.1; $R\in [0,5]$, $z\in [0,4]$ for example 7.2; $R\in [0,3]$, $z\in [0,2]$ for example 7.3), though the computational domain for any given value of $z_d$ is identical for each example.

For each combination of $\gamma$ and $\varOmega$, we choose $z_d=1.3$, $z_d=0.9$ and $z_d=0.5$, representing the helicopter landing on the bed of sand (recall figure 5). Recalling (6.38), the critical value $\varOmega _c$ of $\varOmega$ that determines whether the boundary data is in case (a) or case (b), as shown in figure 7, depends only on $z_d$, and for the three choices of $z_d$ considered here, we have

(7.3a–c)\begin{equation} \varOmega_c(0.5)\approx 0.7522, \quad \varOmega_c(0.9)\approx 0.2322, \quad \varOmega_c(1.3)\approx 0.1113. \end{equation}

For all parameter choices in examples 7.1 and 7.2, and all in examples 7.3 except one (detailed below), we are in either case (a) or case (b)(i). This is illustrated for the case $\gamma =500$, $\varOmega =0.5$, and $z_d=0.9$ or $z_d=0.5$ (as in example B.1) in figure 8: for $z_d=0.9$, $\varOmega =0.5>\varOmega _c(0.9)$, hence case (a) holds; for $z_d=0.5$, $\varOmega =0.5<\varOmega _c(0.5)$, hence case (b) holds – however, in this case, both the local maximum and the local minimum take values greater than $(\gamma E_s)^{-1}$ and $\gamma ^{-1}$, hence case (b)(i) holds. The boundary data $g$ (shown in figures 8c,d) is qualitatively comparable in each of these cases.

Figure 8. Boundary data $g(X)$ and $|\nabla _h\bar {\phi }(X,0)|^2$, $\gamma ^{-1}=0.002$ and $(\gamma E_s)^{-1}=0.2$ for $\gamma =500$, $\varOmega =0.5$ and (a,c) $z_d=0.9$ (case (a)), and (b,d) $z_d=0.5$ (case (b)(i)). Note that the dotted line $\gamma ^{-1}=0.002$ almost overlays the $X$-axis.

Finally, we note that our expectation, from the physical interpretation of the parameters, is that there will be more sand in the air for large $\gamma$, since large $\gamma$ means that lift dominates gravity in the lifting layer (recalling that $\gamma$ is the ratio of the lift force per unit volume to the gravity force per unit volume on particles in the interfacial layer), so that larger $\gamma$ corresponds to more sand being entrained into the fluidized region. We recall also that $\varOmega$ measures the ratio of the swirl velocity to the downdraught velocity, and recalling the definition of $g$ in (6.24), we expect that this will also lead to more sand being entrained into the air for large $\varOmega$.

Example 7.1 (Large $\gamma$)

We first consider $\gamma =1000$. We plot the voidage field $E$ for $\varOmega =1$ in figures 9(a,d,g), for $\varOmega =0.5$ in figures 9(b,e,h), and for $\varOmega =0.1$ in figures 9(c,f,i), each for $z_d=1.3$ (figures 9a–c), $z_d=0.9$ (figures 9d–f) and $z_d=0.5$ (figures 9g–i). We see more sand in the air (lower values of $E$) corresponding to larger values of $\varOmega$, and to lower values of $z_d$. For $\varOmega =1$, and for $\varOmega =0.5$ and $z_d=0.9, 1.3$, the boundary data fits case (a) and is qualitatively comparable to that seen in figure 8c. For $\varOmega =0.1$, and for $\varOmega =0.5$ and $z_d=0.5$, the boundary data fits case (b)(i) and is qualitatively comparable to that seen in figure 8d. In either case, this tallies with the single plume that we see rising up around the helicopter as $z_d$ decreases.

Figure 9. Voidage field $E$, example 7.1, plotted for $\gamma =1000$ and (a,d,g) $\varOmega =1$, (b,e,h) $\varOmega =0.5$, and (c,f,i) $\varOmega =0.1$, each for (a–c) $z_d=1.3$, (d–f) $z_d=0.9$, and (g–i) $z_d=0.5$. In each plot, $R\in [0,6]$, $z\in [0,5]$.

Example 7.2 (Medium $\gamma$)

Next, we consider $\gamma =500$. We plot the voidage field $E$ for $\varOmega =1$ in figures 10(a,d,g), for $\varOmega =0.5$ in figures 10(b,e,h), and for $\varOmega =0.1$ in figures 10(c,f,i), each for $z_d=1.3$ (figures 10a–c), $z_d=0.9$ (figures 10d–f) and $z_d=0.5$ (figures 10g–i). Again, we see more sand in the air (lower values of $E$) corresponding to larger values of $\varOmega$, and to lower values of $z_d$. Comparing to figure 9, we see less sand in the air for the lower value of $\gamma$. Note that the plotted range is smaller in figure 10 than in figure 9, though the computational domain is identical. Exactly as for example 7.1, for $\varOmega =1$, and for $\varOmega =0.5$ and $z_d=0.9, 1.3$, the boundary data fits case (a), whilst for $\varOmega =0.1$, and for $\varOmega =0.5$ and $z_d=0.5$, the boundary data fits case (b)(i), with similar qualitative solution behaviour.

Figure 10. Voidage field $E$, example 7.2, plotted for $\gamma =500$ and (a,d,g) $\varOmega =1$, (b,e,h) $\varOmega =0.5$ and (c,f,i) $\varOmega =0.1$, each for (a–c) $z_d=1.3$, (d–f) $z_d=0.9$ and (g–i) $z_d=0.5$. In each plot, $R\in [0,5]$, $z\in [0,4]$.

Example 7.3 (Small $\gamma$)

Next, we consider $\gamma =100$. We plot the voidage field $E$ for $\varOmega =1$ in figures 11(a,d,g), for $\varOmega =0.5$ in figures 11(b,e,h) and for $\varOmega =0.1$ in figures 11(c,f,i), each for $z_d=1.3$ (figures 11a–c), $z_d=0.9$ (figures 11d–f) and $z_d=0.5$ (figures 11g–i). Again, we see more sand in the air (lower values of $E$) corresponding to larger values of $\varOmega$, and to lower values of $z_d$. Comparing to figures 9 and 10, we see less sand in the air for the lower value of $\gamma$. Note that the plotted range is smaller in figure 11 than in both figures 10 and 9, though the computational domain is identical in each case. As for examples 7.1 and 7.2, for $\varOmega =1$, and for $\varOmega =0.5$ and $z_d=0.9, 1.3$, the boundary data fits case (a), whilst for $\varOmega =0.5$ and $z_d=0.5$, and for $\varOmega =0.1$ and $z_d=0.9,1.3$, the boundary data fits case (b)(i), with similar qualitative solution behaviour. However, for $\varOmega =0.1$ and $z_d=0.5$, we are now in case (b)(iii), as illustrated in figure 12. In this case, we see that the local maximum is greater than $(\gamma E_s)^{-1}$, whilst the local minimum is less than $(\gamma E_s)^{-1}$. As a result, the boundary data $g$ has a small hump near $X\approx 0.05$, though it is hard to see any discernible effect from this on the qualitative solution behaviour (figure 11i).

Figure 11. Voidage field $E$, example 7.3, plotted for $\gamma =100$ and (a,d,g) $\varOmega =1$, (b,e,h) $\varOmega =0.5$ and (c,f,i) $\varOmega =0.1$, each for (a–c) $z_d=1.3$, (d–f) $z_d=0.9$ and (g–i) $z_d=0.5$. In each plot, $R\in [0,3]$, $z\in [0,2]$.

Figure 12. Boundary data $g(X)$ and $|\nabla _h\bar {\phi }(X,0)|^2$, $\gamma ^{-1}=0.01$ and $(\gamma E_s)^{-1}=1$ for $\varOmega =0.1$, $z_d=0.5$, example 7.3. Note that the dotted line $\gamma ^{-1}=0.01$ almost overlays the $X$-axis.

Example 7.4 (Two-plume example)

In our final example, we choose

(7.4a–d)\begin{equation} E_s = 0.01,\quad \varDelta = 0.04, \quad \gamma=70, \quad \varOmega=0.1, \end{equation}

with first $\alpha = 0.1$ (as above), and then $\alpha =0.02$. Our choice of $\gamma =70$ is smaller than in examples 7.1–7.3, hence $\gamma ^{-1}$ is larger, enabling us to demonstrate case (b)(ii) (recalling figure 7), which provides new and qualitatively different phenomena in the solution behaviour compared to that seen in examples 7.1–7.3.

Here, we will start with $z_d=1.1$, and gradually reduce it, with the transition between cases (b)(i) and (b)(ii) being rather sensitive to the choice of $z_d$, for $z_d\in [0.9,1.1]$. This is illustrated in figure 13, where we plot $|\nabla _h\bar {\phi }(X,0)|^2$, the value $\gamma ^{-1}$ and the boundary data $g$ for a range of values of $z_d$. We see that for $z_d=1.1$, the local maximum is below $\gamma ^{-1}$ (case (b)(i)). As $z_d$ decreases, the local maximum becomes greater than $\gamma ^{-1}$, causing a dip in the boundary data, where $g<1$ between two regions where $g=1$ (case (b)(ii)). As $z_d$ decreases further, to $z_d=0.98$ and below, the local minimum rises above $\gamma ^{-1}$, and we return to case (b)(i). The presence of two separated regions where $g=1$ drives the formation of a qualitatively new phenomenon, compared to examples 7.1–7.3, namely a second ‘plume’.

Figure 13. Boundary data $g(X)$, and $|\nabla _h\bar {\phi }(X,0)|^2$, $\gamma ^{-1}=1/70$, example 7.4, plotted for $z_d=1.1, 1.05, 1.02, 1.0, 0.98, 0.95, 0.9, 0.8$. For $z_d=1.0,1.02,1.05$, we are in case (b)(ii); for all other values of $z_d$, we are in case (b)(i). The line $(\gamma E_s)^{-1}=0.7$ is not shown on the plot.

Specifically, for $z_d=1.1$, we see a single concentration of sand below the helicopter. As the helicopter descends to $z_d=1.05$, a second concentration of sand forms; as the helicopter descends further, this second concentration grows and starts to form a second ‘plume’ (compared to the single plume that we saw in examples 7.1–7.3), until at approximately $z_d=1.0$, the two concentrations of sand merge, forming a single plume again. This behaviour is illustrated for $\alpha =0.1$ in figure 14. All experiments in this example were carried out with $N=8$, and with $L_R=16N$, $L_z=8N$ (lower values than for examples 7.1–7.3, with the key phenomena more localized near the origin), and again $\epsilon =1/(10N^2)$. As for the examples above, we again truncate the computational domains for presentational purposes, noting that the solution satisfies $E\approx 1$ outside the plotted range.

Figure 14. Voidage field $E$, example 7.4, $\alpha =0.1$, $z_d\in [0.4,1.1]$.

Finally, we repeat the same experiment with $\alpha =0.02$ instead of $\alpha =0.1$. This smaller value of $\alpha$, compared to all other examples above, represents the case where convection is dominating more than diffusion, and we see in figure 15 that this leads to a ‘wispier’ (thinner) plume, compared to the earlier examples. Again, we see the formation of a second ‘plume’, for the same reasons as explained above for the case $\alpha =0.1$.

Figure 15. Voidage field $E$, example 7.4, $\alpha =0.02$, $z_d\in [0.4,1.1]$.

8. Conclusion

In this paper, we have developed a model for the formulation and analysis of problems relating to the flow of an incompressible fluid above an otherwise static particle bed, and the consequent development and evolution of a fluidized particle cloud in the fluid flow. With the flow in the bulk fluidized region modelled as a two-phase flow, the principal contribution of the paper has been to develop a rational theory for the key processes of particle entrainment and detrainment from the otherwise static particle bed into the fully fluidized region above. This leads to a natural additional macroscopic boundary condition at the interface of the fully fluidized region and the particle bed, which renders a closed problem in the fully fluidized region, for the voidage field and a fluid velocity potential.

To illustrate the applicability of this theory, we have formulated it in relation to a simple model of the helicopter cloud problem, leading to a fully determined nonlinear elliptic boundary value problem for the voidage field throughout the fluidized region. This provides a rational and efficient approach to computing the voidage field for this application. Careful numerical solution via a finite-difference approximation scheme has demonstrated through a series of experiments how the qualitative behaviour of the voidage field varies as key parameters change, in a way that is entirely consistent with our understanding of the underlying physical processes. It may be of interest to look at this example in relation to controlled scale model experiments in order to specifically determine the qualitative and quantitative accuracy of the rational model that has been developed in this paper, but such considerations are left to future work. However, we do note that qualitative structural agreement between the high-speed photography presented in Rovere (2022), and particularly relating to figure 2.6 therein, and the theoretical figures presented here in § 7, is very encouraging.