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 (Reference Bagnold1941) 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 (Reference Smith and Palmer2019), Palmer & Smith (Reference Palmer and Smith2020,  Reference Palmer and Smith2021) and Jolley & Smith (Reference Jolley and Smith2022), 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 Reference Needham and Merkin1983). 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 Reference Gilbarg and Trudinger1998, 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 Reference Gilbarg and Trudinger1998, 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. (Reference Wachspress, Whitehouse, Keller, McClure, Gilmore and Dorsett2008), Govindarajan et al. (Reference Govindarajan, Leishman and Gumerov2013), Porcù et al. (Reference Porcù, Miglio, Parolini, Penati and Vergopolan2020), Tan et al. (Reference Tan, Gao, Barakos, Lin, Zhang and Huang2021,  Reference Tan, Ge, Zhang, Cui and Wang2022a,Reference Tan, Yon, He, Yu and Wangb), Li et al. (Reference Li, Shi, Xu and Wu2023) and Lin et al. (Reference Lin, Xu, Jiang, Hu and Lee2024), 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]$.