2.1. Preliminary considerations

The cold flow surrounding fire whirls, to be described using cylindrical polar coordinates $(r^*,\theta ,z^*)$ and associated velocity components $(u^*,v^*,w^*)$, is induced by the entrainment of the turbulent plume that extends vertically above the flame along the axis of symmetry. The volumetric entrainment rate per unit length, taken from investigations free from fire-whirl swirl, increases with the two-thirds power of the vertical distance $z^*$ according to $\Phi = 2{\rm \pi} C B^{1/3} z^{*2/3}$, where $B$ is the specific buoyancy flux (Batchelor Reference Batchelor1954) and $C$ is a dimensionless factor, which approximately assumes the value $C=0.041$, as suggested by experimental results (Rouse, Yih & Humphreys Reference Rouse, Yih and Humphreys1952; List Reference List1982). Correspondingly, the flow induced in the axial plane has velocities decaying with the radial distance according to $C(B/r^*)^{1/3}$. For $r^* \gg [\nu /(CB^{1/3})]^{3/2}$ the associated Reynolds number $C B^{1/3} r^{*2/3}/\nu$ is large, resulting in nearly inviscid motion, which, in the absence of swirl, is described by a self-similar potential solution that is due to Taylor (Reference Taylor1958). The corresponding slip velocity at the wall is given by

(2.1)\begin{equation} u^*_w={-}A_T C (B/r^*)^{1/3}, \end{equation}

involving the numerical factor

(2.2)\begin{equation} A_T= \frac{4 (2)^{1/3} {\rm \pi}^2}{3{\rm \Gamma}^3(1/3)} \simeq 0.8624, \end{equation}

where ${\rm \Gamma}$ denotes the gamma function. The potential solution fails in the boundary layer, where the radial velocity, also self-similar, is given by $u^*/u^*_w=f_{T}'$ in terms of the derivative of the reduced streamfunction $f_{T}(\varsigma )$, a function of the rescaled vertical distance $\varsigma =(A_T C B^{1/3}/\nu )^{1/2} r^{* -2/3} z^*$ determined from the boundary-value problem

(2.3a,b)\begin{equation} f_{T}'''-\frac{4}{3}f_{T} f_{T}''+\frac{1}{3} (1-f_{T}^{\prime 2})=0, \quad f_{T}(0)=f_{T}'(0)=f_{T}'(\infty)-1=0, \end{equation}

the subscript $T$ referring to the boundary layer accompanying Taylor's potential flow. In the notation employed throughout the paper the prime denotes differentiation of functions of one variable (e.g. in the above description, it represents differentiation with respect to the self-similar coordinate $\varsigma$). It is worth mentioning that the flow structure surrounding turbulent plumes, relevant to fire whirls, is fundamentally different from that surrounding laminar plumes, characterized by small entrainment rates $\Phi \sim \nu$ and associated flow Reynolds numbers of order unity, the case analysed by Schneider (Reference Schneider1981), who found an exact self-similar solution of the first kind for the swirl-free flow in the axial plane. As shown recently by Coenen et al. (Reference Coenen, Rajamanickam, Weiss, Sánchez and Williams2019b), the accompanying circulation in this viscous case is described by a self-similar solution of the second kind.

As previously mentioned, the three-dimensional inviscid motion surrounding fire whirls is affected by the manner in which swirl is imparted to the flow. Instead of focusing on a specific configuration, for generality in the following analysis the inviscid flow outside the boundary layer will be described using as a canonical model the exact axisymmetric solution of the Euler equations resulting from combining Taylor's potential solution for the meridional flow and a line vortex of circulation $2 {\rm \pi}\Gamma$ for the azimuthal flow. Correspondingly, at the outer edge of the boundary layer the radial velocity approaches the value given in (2.1) while the azimuthal component approaches the value

(2.4)\begin{equation} v^*_w=\Gamma/r^*. \end{equation}

The presence of swirl alters the flow across the boundary layer, so that, even for this model problem, a self-similar description, which is available in the absence of swirl, as described above in (2.3a,b), does not exist when $v^*_w \ne 0$. The fundamental lack of flow similarity can be illustrated by considering the flow at large radial distances, where the radial motion becomes dominant, as can be inferred from the different decay rates present in (2.1) and (2.4). Correspondingly, as $r^* \rightarrow \infty$ the self-similar function $u^*/u^*_w=f_{T}'(\varsigma )$ appears to be the appropriate leading-order representation for the radial velocity across the boundary layer, while the azimuthal velocity $v^*=\Gamma g_{T}(\varsigma )/r^*$ should be determined by the accompanying problem

(2.5a,b)\begin{equation} g_{T}''-\frac{4}{3} f_{T} g_{T}'=0, \quad g_{T}(0)=g_{T}(\infty)-1=0, \end{equation}

obtained at leading order from the axisymmetric boundary-layer form of the azimuthal momentum equation. This presumed self-similar structure fails, however, because the last problem has no solution, which can be seen by investigating the behaviour as $\varsigma \to \infty$ of the first integral $g_{T}'=g_{T}'(0) \exp [\frac {4}{3} \int _0^\varsigma f_{T} \,\textrm {d} \tilde {\varsigma } ]$ to show that $g_{T}'(0)=0$ to avoid divergence, so that the only possible solution is $g_{T}=$ constant, which cannot satisfy simultaneously both boundary conditions $g_{T}(0)=g_{T}(\infty )-1=0$; this lack of similarity is also encountered when the outer flow is driven solely by a potential vortex (Gol'dshtik Reference Gol'dshtik1960) (see also King & Lewellen (Reference King and Lewellen1964) for a discussion of boundary-layer self-similarity when the outer azimuthal velocity varies with a general power of the radial distance).

Progress in understanding can be achieved by investigating the development of the boundary layer from a given radial location, as was done in the previous analysis of the boundary layer on a disk of radius $a$ (Burggraf et al. Reference Burggraf, Stewartson and Belcher1971). The same approach is to be considered below, with the ratio of the radial-to-azimuthal velocity

(2.6)\begin{equation} \sigma=\frac{u^*_w(a)}{v^*_w(a)} \end{equation}

at the disk edge arising as the only controlling parameter in the resulting description. This idealized disk problem may, for example, be considered to provide an approximate description of the main features of the boundary-layer flow in the region between the fire and swirl-producing vanes at radius $a$ in laboratory fire-whirl experiments, with the velocity ratio $\sigma$ being directly related to the angle of inclination of the vanes. In particular, the terminal velocity profile at $r^* \ll a$ can be anticipated to provide a realistic representation for the flow surrounding localized fire whirls, with the parameter $\sigma$ measuring the level of swirl introduced by the collective effect of the flow-deflecting obstacles, located at radial distances much larger than the characteristic size of the fuel source feeding the fire.

2.2. Problem formulation

Following Burggraf et al. (Reference Burggraf, Stewartson and Belcher1971), the problem is scaled using $a$ and

(2.7)\begin{equation} \delta=\frac{a}{\sqrt{{\textit{Re}}}} \end{equation}

for the radial and axial coordinates, with

(2.8)\begin{equation} {\textit{Re}}=\frac{\Gamma}{\nu} \end{equation}

representing the relevant Reynolds number. Correspondingly, the azimuthal and radial velocity components $u^*$ and $v^*$ are scaled with $\Gamma /a$, corresponding to a radial pressure gradient scaling with $\rho \Gamma ^2/a^3$ ($\rho$ representing the density), while the axial component $w^*$ is scaled with $(\Gamma /a)/\sqrt {{\textit {Re}}}$, resulting in the dimensionless variables

(2.9a–e)\begin{equation} r = \frac{r^*}{a}, \quad z= \frac{z^*}{\delta}, \quad u = \frac{u^*}{\Gamma/a}, \quad v = \frac{v^*}{\Gamma/a}, \quad w= \frac{w^*}{\Gamma/(a \sqrt{{\textit{Re}}})}. \end{equation}

Neglecting terms of order ${\textit {Re}}^{-2} \ll 1$ reduces the conservation equations, written in their steady axisymmetric form for a constant-density fluid, to their boundary-layer form

(2.10)\begin{gather} u{\frac{\partial u}{\partial r}} + w{\frac{\partial u}{\partial z}} - \frac{v^2}{r} ={-}\frac{\sigma^2}{3r^{5/3}} -\frac{1}{r^3} + {\frac{\partial ^2 u}{\partial z^2}}, \end{gather}

(2.11)\begin{gather}u{\frac{\partial }{\partial r}} (rv) + w{\frac{\partial }{\partial z}} (rv) = {\frac{\partial ^2}{\partial z^2}}(rv), \end{gather}

(2.12)\begin{gather}{\frac{\partial }{\partial r}}(ru) + {\frac{\partial }{\partial z}}(rw) =0. \end{gather}

The first two terms on the right-hand side of (2.10) arise from the radial pressure gradients imposed by the external Taylor and potential-vortex flows, respectively. These equations are to be integrated for decreasing values of $r$ with the boundary conditions

(2.13a,b)\begin{equation} u \rightarrow -\sigma r^{{-}1/3}, \quad v\rightarrow \frac{1}{r} \quad \text{as} \ z\to \infty \end{equation}

and

(2.14)\begin{equation} u=v=w=0 \quad \text{at} \ z=0 \end{equation}

for $r < 1$ and the initial velocity profiles $u=-\sigma$ and $v=1$ at $r=1$, consistent with (2.13a,b). The only parameter in the description is the initial flow inclination $\sigma$. As expected, for $\sigma =0$ the problem reduces exactly to that addressed by Burggraf et al. (Reference Burggraf, Stewartson and Belcher1971).

2.3. Sample numerical results

The problem (2.10)–(2.14) was integrated numerically by marching from $r=1$ with decreasing values of $r$ for selected values of $\sigma$. Sample results are shown in figure 2 for $\sigma =1$. The radial and azimuthal velocities are uncoupled for $1-r \ll 1$, when the effects of the centripetal acceleration $-v^2/r$ and pressure gradient $-\sigma ^2 r^{-3/5}/3-r^{-3}$ can be neglected in (2.10) at leading order, reducing the solution with $\sigma \ne 0$ to $-u/(\sigma r^{-1/3})=r v=f_{B}'(\zeta )$, where $f_{B}$ is the Blasius streamfunction, obtained by integration of $f_{B}'''+f_{B} f_{B}''/2=0$ with boundary conditions $f_{B}(0)=f_{B}'(0)=f_{B}'(\infty )-1=0$, with the prime denoting here differentiation with respect to the local self-similar coordinate $\zeta =z/(1-r)^{1/2}$. The asymptotic predictions for $1-r \ll 1$ are compared in figure 2 with the profiles obtained numerically at $r=0.9$.

Figure 2. Boundary-layer profiles of radial and azimuthal velocity at various radial locations for $\sigma =1$. The dashed curves represent the asymptotic predictions $-u/(\sigma r^{-1/3})=r v=f_{B}'(\zeta )$ for $1-r \ll 1$ evaluated at $r=0.9$.

The effect of the azimuthal motion on the radial flow is no longer negligible as $1-r$ increases to values of order unity, leading to an overshoot in the radial velocity, as is already evident in the results of figure 2 for $r=0.5$. This overshoot becomes more prominent as $r \rightarrow 0$, with the peak value of $r u$ reaching a near-unity value at small distances $z \sim r$. A detailed view of this near-wall region is shown in figure 3, where the dashed curves represent analytic results, to be developed below.

Figure 3. Comparison of the near-wall boundary-layer profiles obtained from the asymptotic predictions $ru = - \psi _0'(\eta )$ and $rv=C_1 r^\lambda \gamma _1(\eta )$ (dashed curves) with those determined numerically at various radial locations by integration of (2.10)–(2.12) for $\sigma = 1$ (solid curves).

The profiles of $r u$ and $r v$ shown in figure 2 are seen to approach a terminal shape as $r \rightarrow 0$. Although the specific shape of these terminal profiles depends on the value of $\sigma$, all solutions show a common multi-layered asymptotic structure, which, with the exception of the outermost layer, is fundamentally similar to that described in Burggraf et al. (Reference Burggraf, Stewartson and Belcher1971) for the potential vortex ($\sigma =0$). A detailed analysis of the different layers is given below, and the associated asymptotic solutions are combined to generate composite expansions for the radial and azimuthal velocity components, providing an accurate boundary-layer description for $r \ll 1$ (i.e. dimensional distances $r^* \ll a$).

3.1. The lower viscous sublayer

At leading order, the solution in the viscous sublayer, where $rv \ll -ru$ (as is apparent from figure 3), is independent of $\sigma$ and corresponds to that described by Burggraf et al. (Reference Burggraf, Stewartson and Belcher1971). With circulation neglected, the boundary–layer equation (2.10) can be expressed in terms of the self-similar coordinate $\eta = {z}/{r}$ and accompanying streamfunction $\psi = r \psi _0(\eta )$, defined such that $ru = - \psi _0'$ and $rw=\psi _0-\eta \psi '_0$, to give at leading order the problem

(3.1a,b)\begin{equation} \psi_0''' - \psi_0 \psi_0'' - \psi_0'^2 + 1 =0, \quad \psi_0(0)=\psi_0'(0)=\psi_0'(\infty) -1 =0. \end{equation}

The solution, expressible in terms of the parabolic cylinder functions (Mills Reference Mills1935), provides the asymptotic behaviour $\psi _0 \simeq \eta -1.0864$ and

(3.2)\begin{equation} rw \rightarrow -1.0864 \end{equation}

for $\eta \gg 1$. The accompanying weak azimuthal motion is described by a self-similar solution of the second kind of the form $rv \propto r^\lambda \gamma _1(\eta )$, where the eigenfunction $\gamma _1(\eta )$ obeys the linear equation

(3.3)\begin{equation} \gamma_1'' - \psi_0 \gamma_1' + \lambda \psi_0'\gamma_1 = 0 \end{equation}

stemming from (2.11). A non-trivial solution satisfying the non-slip condition $\gamma _1=0$ at $\eta =0$ and exhibiting algebraic growth $\gamma _1 \propto \eta ^\lambda$ (as opposed to exponential growth) as $\eta \rightarrow \infty$, as needed to enable matching with the outer profile, exists only for a discrete set of values of the eigenvalue $\lambda$, with the smallest eigenvalue, corresponding to the dominant eigenfunction for small $r$, found to be $\lambda = 0.6797$.

The solution to the above eigenvalue problem determines the near-wall azimuthal velocity $r v = C_1 r^\lambda \gamma _1(\eta )$, aside from a multiplicative factor $C_1$, a function of $\sigma$ to be determined by matching with the outer inviscid solution. For definiteness, it is convenient to use $\gamma _1 = \eta ^\lambda$ as $\eta \gg 1$ as the normalization condition for the eigenfunction $\gamma _1$, resulting in the asymptotic behaviour

(3.4)\begin{equation} r v \rightarrow C_1 z^\lambda \quad \textrm{as} \ \eta \rightarrow \infty, \end{equation}

to be employed below in the matching procedure. The near-wall asymptotic predictions $ru = - \psi _0'(\eta )$ and $r v = C_1 r^\lambda \gamma _1(\eta )$ for $r \ll 1$ are compared in figure 3 with the results of numerical integration for $\sigma =1$. The predictions for $r v = C_1 r^\lambda \gamma _1(\eta )$ are computed with $C_1=0.9065$, the value obtained below by matching with the outer inviscid results when $\sigma =1$. The asymptotic predictions and the numerical results are seen to be virtually indistinguishable for $r=0.01$ in figure 3.

3.2. The main inviscid layer

In the intermediate layer $z \sim 1$, the expansions for the velocity components at $r \ll 1$ take the form

(3.5)\begin{equation} \left.\begin{array}{c@{}} ru= F_0(z) + rF_1(z) +\cdots,\\ rv= G_0(z) + rG_1(z) + \cdots,\\ rw= H_1(z) + \cdots. \end{array}\right\} \end{equation}

The leading-order functions $F_0(z)$ and $G_0(z)$ must satisfy $F_0 \rightarrow -1$ and $z^{-\lambda } G_0 \rightarrow C_1$ as $z \rightarrow 0$, corresponding to matching with the viscous sublayer, and as $z \rightarrow \infty$ they must approach the outer values $F_0=G_0-1=0$, consistent at this order with the velocity found outside the boundary layer. The two functions $F_0$ and $G_0$ are related by the equation

(3.6)\begin{equation} F_0^2 + G_0^2=1, \end{equation}

which follows from the leading-order $1/r^3$ terms in (2.10), but, other than that, their specific shape depends on the development of the boundary layer for $0< r < 1$, yielding different profiles for different values of $\sigma$. The additional functions $F_1$, $G_1$ and $H_1$ appearing in (3.5) are related to the leading-order functions by

(3.7a–c)\begin{equation} H_1=1.0864 F_0,\quad F_1={-} 1.0864 F'_0,\quad G_1={-} 1.0864 G'_0 \end{equation}

as can be seen by carrying the asymptotic solution to a higher order, with the numerical factor $1.0864$ selected to ensure inner matching with the vertical velocity (3.2).

To determine $G_0$ and $F_0=-\sqrt {1-G_0^2}$, the numerical integration of (2.10)–(2.12) was extended to extremely small radial distances $r \sim 10^{-4}$, and the asymptotic predictions $ru= F_0(z)- 1.0864 r F'_0(z)$ and $rv= G_0(z)- 1.0864 r G'_0(z)$ were used to extrapolate the result to $r=0$. The solution was further corrected to remove the viscous sublayer by replacing the solution at $z \ll 1$ with the near-wall behaviour

(3.8a,b)\begin{equation} F_0 ={-}1 + \frac{1}{2} C_{1}^2 z^{2\lambda}, \quad G_0 = C_{1} z^{\lambda}, \end{equation}

arising from matching with (3.4), with the constant $C_1$ obtained from the numerical integrations by evaluating $z^{-\lambda } r v$ at small distances $z \sim r$ from the wall, yielding for instance $C_1=(0.7125,0.8275,0.9065,0.9596,1.0181,1.6187)$ for $\sigma =(5,2,1,0.5,0.2,0.1)$. The observed evolution for decreasing values of $\sigma$ appears to be in agreement with the limiting value $C_1=1.6518$ reported by Burggraf et al. (Reference Burggraf, Stewartson and Belcher1971) for $\sigma =0$. Sample profiles of F ₀ and G ₀ are shown in figure 4 for selected values of $\sigma$.

Figure 4. Profiles of $G_0$ and $F_0$ for various $\sigma$. Shown as dashed line are the profiles of Burggraf et al. (Reference Burggraf, Stewartson and Belcher1971) for $\sigma =0$.

While the asymptotic behaviour $G_0=\sqrt {1-F_0^2} \propto z^{\lambda }$ as $z \rightarrow 0$ shown in (3.8a,b) applies for both $\sigma =0$ and $\sigma \ne 0$, the solution as $z \rightarrow \infty$ is qualitatively different in these two cases, with the exponential decay found by Burggraf et al. (Reference Burggraf, Stewartson and Belcher1971) for $\sigma =0$ being replaced for $\sigma \ne 0$ by an algebraic decay of the form

(3.9a,b)\begin{equation} F_0 ={-}C_{2} z^{-\mu} + \cdots, \quad G_0 = 1 - \frac{1}{2} C_{2}^2 z^{{-}2\mu} + \cdots, \end{equation}

where the factor $C_2$ and the exponent $\mu$ depend on $\sigma$. Their values were obtained from the numerical results by examining the decay with vertical distance of the near-axis terminal profiles of $ru$ and $rv$, yielding for instance $C_2=(2.85,1.50,0.88,0.47)$ and $\mu =(1.21,1.19,1.16,1.11)$ for $\sigma =(5,2,1,0.5)$. At a given $r \ll 1$, the range of $z$ over which (3.9a,b) applies decreases for decreasing $\sigma$, thereby hindering the precise evaluation of $C_2$ and $\mu$ for $\sigma < 0.5$. The observed evolution of the approximate values, computed with use of the velocity profiles at the smallest radial distance reached in the boundary-layer computations (i.e. $r \simeq 10^{-4}$), indicates that the exponent $\mu$ decreases with decreasing $\sigma$ to approach unity as $\sigma \rightarrow 0$, with the accompanying value of $C_2$ vanishing in this limit.

3.3. The upper transition region

The asymptotic expansion (3.5) fails in a transition region corresponding to $z\sim r^{-2/(3\mu )} \gg 1$ where, according to (3.9a,b), the value of $F_0$ becomes of order $F_0 \sim r^{2/3}$, comparable to the limiting value $ru=-\sigma r^{2/3}$ found at $z=\infty$. This transition region can be described in terms of the order-unity similarity coordinate $\xi =(\sigma /C_2)^{1/\mu }r^{{2}/{3\mu }}z$ and associated rescaled velocity variables

(3.10)\begin{equation} \left.\begin{array}{c@{}} ru= \sigma r^{2/3} f(\xi),\\ rv= 1 + \sigma^2 r^{4/3}g(\xi),\\ rw=C_2^{1/\mu}\sigma^{{(\mu-1)}/{\mu}}r^{-{(2+\mu)}/{3\mu} } h(\xi). \end{array}\right\} \end{equation}

At leading order, the boundary-layer equations (2.10)–(2.12) in this inviscid outer transition region simplify to

(3.11)\begin{gather} \left(\frac{2\xi}{3\mu}f + h\right)f' = \frac{f^2-1}{3} + 2g, \end{gather}

(3.12)\begin{gather}\left(\frac{2\xi}{3\mu}f + h\right)g' ={-}\frac{4fg}{3}, \end{gather}

(3.13)\begin{gather}\left(f + \frac{\xi}{\mu}f'\right) + \frac{3h'}{2} = 0. \end{gather}

The solution can be reduced to a quadrature as follows. Dividing (3.11) by (3.12) and integrating the resulting equation using the boundary conditions $f(\infty )\rightarrow -1$ and $g(\infty )\rightarrow 0$, which follow from matching with the outer potential solution, provides

(3.14)\begin{equation} g = \frac{1}{2} (1-f^2). \end{equation}

Substitution of this result into (3.11) and elimination of $h$ with use of (3.13) leads to the autonomous equation

(3.15)\begin{equation} (\,f^2 - 1)f'' - \frac{\mu+1}{\mu} ff'^2 =0, \end{equation}

which can be integrated once using the boundary condition $f \rightarrow -\xi ^{-\mu }$ as $\xi \rightarrow 0$, obtained by matching with (3.9a,b), to give $f'=\mu (\,f^2-1)^{(\mu +1)/(2\mu )}$, finally yielding

(3.16)\begin{equation} \mu \xi = \int_{-\infty}^f \frac{\textrm{d} \tilde{f}}{(\,\tilde{f}^2-1)^{{(\mu+1)}/{(2\mu)}}}, \end{equation}

with $\tilde {f}$ being a dummy integration variable. The above integral, which is expressible in terms of incomplete beta function, provides, together with the previous equation (3.14), the radial and azimuthal velocity distributions $f(\xi )$ and $g(\xi )$.

Inspection of (3.16) reveals that, for the values $\mu >1$ that apply in our description, the function $f$ is a front solution that reaches the boundary value $f=-1$ at a finite location $\xi =\xi _o$ given by

(3.17)\begin{equation} \xi_o = \frac{1}{2\mu}\textrm{B}\left(\frac{1}{2\mu}, \frac{\mu-1}{2\mu}\right), \end{equation}

as follows directly from (3.16), with $\textrm {B}$ representing here the beta function. Note that, since the value of $\mu -1$ remains relatively small for $\sigma \sim 1$, the front is always located at large distances $\xi _o \simeq 1/(\mu -1)$.

3.4. The composite expansion

The separate solutions found in the different regions can be used to generate the composite expansions

(3.18)\begin{equation} \left.\begin{array}{c@{}} ru ={-} \psi^{\prime}_0(z/r) + 1 + F_0(z) + \sigma r^{2/3} f[(\sigma/C_2)^{1/\mu}r^{{2}/{(3\mu)}}z] + C_2 z^{-\mu},\\ rv = C_1 r^\lambda \gamma_1(z/r) - C_1 z^\lambda + G_0(z) + \sigma^2 r^{4/3} g[(\sigma/C_2)^{1/\mu}r^{{2}/{(3\mu)}}z] + \tfrac{1}{2} C_2^2 z^{{-}2\mu}, \end{array}\right\} \end{equation}

which describe the radial and azimuthal velocity profiles as $r\to 0$ with small errors of order $r$. The accuracy of these expansions is tested in figure 5 by comparing the asymptotic predictions with the results of numerical integrations of the boundary-layer equations (2.10)–(2.12). The degree of agreement displayed in the figure is clearly satisfactory, with the composite expansion being virtually indistinguishable from the numerical results at $r=0.01$.

Figure 5. Profiles of radial and azimuthal velocity for $\sigma =0.5$ (a–c) and $\sigma = 2.0$ (d–f) obtained at different radial locations $r\ll 1$ by numerical integration of (2.10)–(2.12) (solid curves) and by evaluation of the composite expansion (3.18) (dashed curves).

4.1. The rescaled problem

In the non-slender collision region, of characteristic size $\delta =a/\sqrt {{\textit {Re}}}$, all three velocity components have comparable magnitudes $u^* \sim v^* \sim w^* \sim \Gamma /\delta$. Correspondingly, the analysis of this region necessitates introduction of rescaled velocity components $\tilde {u}=u^*/(\Gamma /\delta )$, $\tilde {v}=v^*/(\Gamma /\delta )$ and $\tilde {w}=w^*/(\Gamma /\delta )$ along with a rescaled radial coordinate $\tilde {r}=r^*/\delta$, while the accompanying vertical coordinate $\tilde {z}=z^*/\delta =z$ is that used in the boundary-layer analysis. With these scales, the steady, axisymmetric continuity and momentum equations take the dimensionless form

(4.1)\begin{gather} \frac{1}{\tilde{r}} \frac{\partial }{\partial \tilde{r}}(\tilde{r} \tilde{u})+\frac{\partial \tilde{w}}{\partial \tilde z}=0, \end{gather}

(4.2)\begin{gather}\tilde{u} \frac{\partial \tilde{u}}{\partial \tilde{r}}-\frac{\tilde{v}^2}{\tilde{r}}+\tilde{w} \frac{\partial \tilde{u}}{\partial \tilde z}={-}\frac{\partial \tilde{p}}{\partial \tilde{r}}+\frac{1}{{\textit{Re}}} \left[ \frac{\partial}{\partial \tilde{r}} \left( \frac{1}{\tilde{r}} \frac{\partial}{\partial \tilde{r}}(\tilde{r} \tilde{u}) \right)+\frac{\partial^2 \tilde{u}}{\partial \tilde z^2}\right], \end{gather}

(4.3)\begin{gather}\tilde{u} \frac{\partial}{\partial \tilde{r}} (\tilde{r} \tilde{v}) +\tilde{w} \frac{\partial }{\partial \tilde z}(\tilde{r}\tilde{v})=\frac{1}{{\textit{Re}}} \left[\tilde{r} \frac{\partial}{\partial \tilde{r}} \left( \frac{1}{\tilde{r}} \frac{\partial}{\partial \tilde{r}}(\tilde{r} \tilde{v}) \right)+\frac{\partial^2 }{\partial \tilde z^2}(\tilde{r} \tilde{v})\right], \end{gather}

(4.4)\begin{gather}\tilde{u} \frac{\partial \tilde{w}}{\partial \tilde{r}}+\tilde{w} \frac{\partial \tilde{w}}{\partial \tilde z}={-}\frac{\partial \tilde{p}}{\partial \tilde z}+\frac{1}{{\textit{Re}}} \left[\frac{1}{\tilde{r}}\frac{\partial}{\partial \tilde{r}} \left( \tilde{r}\frac{\partial \tilde{w}}{\partial \tilde{r}} \right)+\frac{\partial^2 \tilde{w}}{\partial \tilde z^2}\right], \end{gather}

where $\tilde {p}$ denotes the spatial pressure variations scaled with the characteristic dynamic pressure $\rho (\Gamma /\delta )^2$. The distribution of radial and azimuthal velocity at large radial distances is given by the terminal profiles (3.18) written in terms of the rescaled variables. The solution depends on the Reynolds number ${\textit {Re}} =\Gamma /\nu$, which appears explicitly in the equations, and on the ambient swirl level, through the parameter $\sigma$ present in the boundary velocity profiles (3.18).

For the large values ${\textit {Re}} \gg 1$ considered here, the flow can be expected to be nearly inviscid, although rotational, with viscous effects largely confined to a near-wall layer and to a near-axis core region, both of characteristic size $\delta _v=\delta /\sqrt {{\textit {Re}}} \ll \delta$. The inviscid solution is to be investigated in detail below for different values of $\sigma$. To check for consistency of the large-Reynolds-number structure, additional attention is given to the accompanying near-wall boundary layer. The computations reveal that boundary-layer separation occurs at a finite distance $r^* \sim \delta$ regardless of the value of $\sigma$, indicating that the description of the corner region should, in principle, account for viscous effects. This finding motivates additional Navier–Stokes computations for moderately large values of ${\textit {Re}}$, which allow us to investigate the extent of the separation region and its dependence on the Reynolds number.

The results to be presented below, extending the previous studies by using boundary velocity profiles that are directly relevant to fire-whirl applications, pertain to cold flow only. Before proceeding with the analysis, it is worth discussing the relevance of the results in connection with localized fire whirls. If one considers, for definiteness, the case of fire whirls developing above liquid-fuel pools, then the fuel-pool diameter $D$ emerges as relevant characteristic length, to be compared with the size of the collision region $\delta$. In the relevant distinguished limit $D \sim \delta$ the fire whirl develops in the collision region, driven by the approaching boundary-layer profile described in (3.18). Since ${\textit {Re}} \gg 1$, the flame, developing from the liquid-pool rim, would be confined initially to the viscous layer $z^* \sim \delta _v=\delta /\sqrt {{\textit {Re}}}$ found in the immediate vicinity of the pool surface and, upon flow deflection near the origin, to the near-axis viscous core found at $r^* \sim \delta /\sqrt {{\textit {Re}}}$. The inviscid results given below provide in this case the velocity profile found outside the thin reactive regions as well as the associated imposed pressure gradient, both along the liquid-pool surface and along the vertical axis. The fire whirl, driven by the fast upward flow resulting from the deflection, would continue to develop vertically over distances larger than $\delta$, eventually transitioning to a turbulent plume. Depending on the flow conditions, buoyancy effects, which eventually drive the turbulent plume, can become important already in the reactive boundary layer developing near the liquid-pool surface, possibly helping to prevent boundary-layer separation. A recent attempt to describe this layer (Li, Yao & Law Reference Li, Yao and Law2019) has employed a constant-density model along with the self-similar velocity profile computed from (3.1a,b). Clearly, more accurate numerical computations, accounting for variable-density and buoyancy effects and using as boundary condition the wall-velocity distribution obtained in the inviscid analysis of the collision region, are worth pursuing in future work.

4.2. The reduced inviscid formulation

As first shown by Hicks (Reference Hicks1899), the inviscid equations that follow from removing the viscous terms involving the factor ${\textit {Re}}^{-1}$ from the momentum equations (4.1)–(4.4) can be combined into a single equation for the streamfunction $\tilde \psi$. As explained by Batchelor (Reference Batchelor1967), the development uses the condition that the circulation per unit azimuthal angle $\tilde {C}=\tilde {r} \tilde {v}$ and the total head $\tilde {H}=\tilde {p}+(\tilde {u}^2+\tilde {v}^2+\tilde {w}^2)/2$ remain constant along any given streamline, allowing the azimuthal component of the vorticity to be written in the form

(4.5)\begin{equation} \frac{\partial \tilde u}{\partial \tilde z}-\frac{\partial \tilde w}{\partial \tilde{r}}=\frac{\tilde{C}}{\tilde{r}} \frac{\textrm{d} \tilde{C}}{\textrm{d} \tilde \psi}-\tilde{r} \frac{\textrm{d} \tilde{H}}{\textrm{d} \tilde \psi}, \end{equation}

finally yielding

(4.6)\begin{equation} {\frac{\partial ^2\tilde \psi}{\partial \tilde r^2}} - \frac{1}{\tilde r}{\frac{\partial \tilde \psi}{\partial \tilde r}} + {\frac{\partial ^2\tilde \psi}{\partial \tilde z^2}} ={-} \tilde{C}\frac{\textrm{d}\tilde{C}}{\textrm{d}\tilde \psi}+\tilde r^2 \frac{\textrm{d} \tilde{H}}{\textrm{d}\tilde \psi}, \end{equation}

upon substitution of the expressions $\tilde {r} \tilde {w}=\partial \tilde \psi /\partial \tilde {r}$ and $\tilde {r} \tilde {u}=-\partial \tilde \psi /\partial \tilde z$. As shown below, the functions $\tilde {C}(\tilde \psi )$ and $\tilde {H}(\tilde \psi )$ are to be evaluated using the terminal velocity profiles (3.18) written in the simplified form $\tilde r \tilde u=F_0(\tilde z)$ and $\tilde r \tilde v=G_0(\tilde z)$, corresponding to ${\textit {Re}} \rightarrow \infty$, with the functions $F_0$ and $G_0$ carrying the dependence on $\sigma$, as shown in figure 4.

Since the streamlines lie parallel to the wall as $\tilde r \to \infty$, the equation $\tilde {r} \tilde {u}=-\partial \tilde \psi /\partial \tilde z$ provides $\tilde \psi =-\int _0^{\tilde z} F_0\,\textrm {d} \tilde z$, which can be used to determine the boundary distribution

(4.7)\begin{equation} \tilde \psi=\tilde \psi_\infty(\tilde z)={-}\int_0^{\tilde z} F_0 \,\textrm{d} \tilde z \quad \textrm{as} \ \tilde r \to \infty, \end{equation}

and implicitly through

(4.8)\begin{equation} \tilde \psi={-}\int_0^{\tilde z_\infty} F_0 \,\textrm{d} \tilde z, \end{equation}

the height $\tilde z_\infty (\tilde \psi )$ at which a given streamline originates. While the head tends to a uniform value as the velocity decays far from the axis, so that $\textrm {d} \tilde {H}/\textrm {d} \tilde \psi =0$ in (4.6), the circulation $\tilde {C}$ varies between streamlines, yielding a contribution to (4.6) that can be evaluated by using $\tilde {C}\,\textrm {d}\tilde {C}/\textrm {d}\tilde \psi =F'_0[\tilde z_\infty (\tilde \psi )]$, derived with use of (3.6). The problem then reduces to that of integrating the nonlinear equation

(4.9)\begin{equation} {\frac{\partial ^2\tilde \psi}{\partial \tilde r^2}} - \frac{1}{\tilde r}{\frac{\partial \tilde \psi}{\partial \tilde r}} + {\frac{\partial ^2\tilde \psi}{\partial \tilde z^2}} = m(\tilde \psi), \end{equation}

where $m(\tilde \psi )=-F'_0[\tilde z_\infty (\tilde \psi )]$, with boundary conditions

(4.10a,b)\begin{equation} \tilde \psi(0,\tilde z) = \tilde \psi(\tilde r,0)=\tilde \psi(\infty,\tilde z) -\tilde \psi_\infty(\tilde z) =0 \quad \text{and} \quad {\frac{\partial \tilde \psi}{\partial \tilde z}} = 0 \quad \text{as} \ \tilde z\rightarrow \infty. \end{equation}

The solution depends on $\sigma$ through the derivative and antiderivative of the function $F_0$, which appear on the right-hand side of (4.9) and in the boundary distribution $\tilde \psi _\infty$ given in (4.7), respectively. Note that, despite the slow algebraic decay $F_0 \simeq -C_{2} \tilde z^{-\mu }$ indicated in (3.9a,b), the condition $\mu > 1$ guarantees that the antiderivative $\int _0^{\tilde z} F_0 \,\textrm {d} \tilde z$ approaches a finite value as $\tilde z \rightarrow \infty$ for all values of $\sigma$. Once $\tilde \psi (\tilde r,\tilde z)$ has been determined, the distribution of azimuthal velocity can be evaluated with use of $\tilde r \tilde v=\tilde {C}(\tilde \psi )=G_0[\tilde z_\infty (\tilde \psi )]$ supplemented by (4.8), as follows in the inviscid limit from conservation of circulation along streamlines. This reduced description is the basis of many of the early vortex-core studies (Rotunno Reference Rotunno1980; Fiedler & Rotunno Reference Fiedler and Rotunno1986; Wilson & Rotunno Reference Wilson and Rotunno1986).

4.3. Sample results

Numerical solutions to the problem defined in (4.9) and (4.10a,b) were obtained using a finite-element method (Hecht Reference Hecht2012). No convergence problems were encountered for any value of $\sigma$. The streamlines and circulation distribution corresponding to $\sigma = 1$ are shown in figure 7(a). The integration provides in particular the slip velocity along the wall $\tilde u_w = - (1/\tilde r){{\partial \tilde \psi }/{\partial \tilde z}}|_{\tilde z=0}$ and the associated radial pressure gradient $-\partial \tilde {p}/\partial \tilde {r}=\tilde u_w \,\textrm {d}\tilde u_w/ \textrm {d} \tilde {r}$, with the latter shown in figure 7(b) for selected values of $\sigma$. In all cases, the pressure gradient, whose magnitude increases with decreasing $\sigma$, is favourable far from the axis and adverse near the axis. The value of $\tilde u_w \,\textrm {d}\tilde u_w/ \textrm {d} \tilde {r}$ is seen to decrease linearly with the radial distance on approaching the origin, a behaviour that is consistent with the local stagnation-point solution $\tilde \psi \sim \tilde r^2\tilde z$ that prevails at $\tilde r^2+\tilde z^2\ll 1$, as follows from a local analysis of (4.9).

Figure 7. The inviscid structure of the collision region calculated from (4.9), including the streamlines $\tilde \psi =(0.25,0.5,0.75,1.0,1.25,1.5,1.75,2.0)$ (solid curves) and circulation $\tilde {r} \tilde {v}$ (colour map) for $\sigma = 1$ (a), the negative radial pressure gradient on the wall obtained from $\tilde u_w(\tilde r) =- (1/\tilde r){{\partial \tilde \psi }/{\partial \tilde z}}|_{\tilde z=0}$ for $\sigma =(0.1,1.0,2.0)$ (b) and the corresponding profiles of axial velocity $\tilde {w}$ (c) and circulation $\tilde {r} \tilde {v}$ (d) approached as $\tilde {z}\rightarrow \infty$.

The deflected streamlines become aligned with the axis for $\tilde z \gg 1$, when the streamfunction approaches the limiting distribution $\Psi (\tilde r)=\tilde \psi (\tilde r,\infty )$, to be determined from integration of

(4.11a,b)\begin{equation} \Psi''-\Psi'/\tilde r=m(\Psi), \quad \Psi(0)=\Psi(\infty)+\int_0^\infty F_0 \,\textrm{d}\tilde z=0, \end{equation}

the one-dimensional counterpart of (4.9). The corresponding distributions of axial velocity $\tilde {w}=\Psi '/\tilde r$ and circulation $\tilde r \tilde {v}=\tilde C[\Psi (\tilde r)]$, which provide the initial conditions for studying the development of the flow above the collision region, are shown in figures 7(c) and 7(d) for the three representative flow inclinations $\sigma =(0.1,1,2)$, for which the boundary values of the streamfunction are $\Psi (\infty )=(1.62,4.16,5.57)$, respectively. As can be seen, the rising jet is found to be wider for increasing $\sigma$, a consequence of the shape of the boundary velocity distributions $F_0$ and $G_0$. The integration provides, in particular, the peak axial velocity $\tilde {w}_0=\tilde {w}(0)$, given by $\tilde w_0=(1.22,1.01,0.86)$ for $\sigma =(0.1,1,2)$. Near the axis, where $m=-C_1^2 \lambda \Psi ^{2\lambda -1}$ and $\tilde {C}=C_1 \Psi ^\lambda$ with $\lambda =0.6797$, as follows from (3.8a,b), the solution takes the form

(4.12a,b)\begin{equation} \tilde{w}=\tilde{w}_0-\frac{C_1^2 \lambda \tilde{w}_0^{2 \lambda-1} \tilde r^{4\lambda-2}}{2^{2\lambda}(2 \lambda-1)} \quad \textrm{and} \quad \tilde r \tilde{v}=C_1 \left(\tilde{w}_0 \tilde r^2/2 \right)^\lambda. \end{equation}

Since $\lambda < 3/4$, the axial velocity of the inviscid solution displays an infinite slope at the axis. This characteristic of the velocity distribution, which would disappear in the presence of viscous forces, is not revealed in the early results of Fiedler & Rotunno (Reference Fiedler and Rotunno1986) because the tabulated representations of $F_0$ and $G_0$ employed in their description, taken from Burggraf et al. (Reference Burggraf, Stewartson and Belcher1971), did not contain enough points to reproduce the near-wall behaviour (3.8a,b).

4.4. The boundary layer in the collision region

The inviscid flow described above is accompanied by a near-wall viscous boundary layer with characteristic thickness $z^*\sim \delta _v=\delta /\sqrt {{\textit {Re}}}$ at radial distances $r^* \sim \delta$. As seen in figure 7(b) this boundary layer develops under the action of a pressure gradient $-\partial \tilde {p}/\partial \tilde {r}=\tilde u_w\, \textrm {d}\tilde u_w/ \textrm {d} \tilde {r}$ that is negative (favourable) at large radial distances but becomes positive (adverse) on approaching the axis. Clearly, the validity of the inviscid solution as a representation of the flow for ${\textit {Re}} \gg 1$ requires that the boundary layer remains attached, that being the assumption underlying previous descriptions (Fiedler & Rotunno Reference Fiedler and Rotunno1986; Wilson & Rotunno Reference Wilson and Rotunno1986; Rotunno Reference Rotunno2013). Examination of this aspect of the problem requires introduction of the rescaled variables

(4.13a–d)\begin{equation} \hat r = \frac{r^*}{\delta}= \tilde r, \quad \hat z = \frac{z^*}{\delta_v}=\sqrt{{\textit{Re}}}\tilde z, \quad \hat u = \frac{u^*}{\Gamma/\delta}= \tilde u, \quad \hat w = \frac{\sqrt{{\textit{Re}}}w^*}{\Gamma/\delta}= \sqrt{{\textit{Re}}} \tilde w \end{equation}

to write the boundary-layer equations

(4.14)\begin{gather} \hat u{\frac{\partial \hat u}{\partial \hat r}} +\hat w{\frac{\partial \hat u}{\partial \hat z}} = \tilde u_w \frac{\textrm{d} \tilde u_w}{\textrm{d} \hat r} + {\frac{\partial ^2 \hat u}{\partial \hat z^2}}, \end{gather}

(4.15)\begin{gather}{\frac{\partial }{\partial \hat r}}(\hat r \hat u) + {\frac{\partial }{\partial \hat z}}(\hat r \hat w) =0, \end{gather}

and associated initial and boundary conditions

(4.16)\begin{equation} \left.\begin{array}{c@{}} \hat r \to \infty:\quad \hat r \hat u ={-}\psi'_0(\hat z/\hat r), \\ \hat z = 0: \quad \hat u = \hat w = 0,\\ \hat z \to \infty:\quad \hat u \to \tilde u_w(\hat r), \end{array}\right\} \end{equation}

involving the apparent slip velocity $\tilde u_w$ of the inviscid collision region, which carries the dependence on $\sigma$, and on the rescaled streamfunction $\psi _0$ across the viscous sublayer, determined from (3.1a,b).

Numerical integration of (4.14)–(4.16) for decreasing values of $\hat r$ reveals that for all $\sigma$ the boundary layer separates at a radial location $\hat r \sim 1$, where the velocity profile develops an inflection point at the wall, preventing integration beyond that point. Illustrative results are shown in figure 8 for $\sigma =1$, for which separation is predicted to occur at $\hat r \simeq 1.76$.

Figure 8. Boundary-layer profiles of radial velocity at selected locations obtained by integration of (4.14) for $\sigma = 1$.

4.5. The viscous structure of the collision region

The predicted separation of the boundary layer, questioning the validity of the inviscid description, was further investigated numerically by integrating the complete Navier–Stokes equations (4.1)–(4.4) with a Newton-Raphson method in combination with the finite-element solver FreeFem++ (Hecht Reference Hecht2012) for increasing values of ${\textit {Re}}$ and different values of $\sigma$. A cylindrical computational domain with outer radius $\tilde r_{max} \gg 1$ and height $\tilde z_{max} \gg 1$ was employed in the integrations. The three-level composite expansions of (3.18), written in terms of the collision-region variables, were used to provide the inlet boundary conditions at $\tilde r=\tilde r_{max}$. Additional boundary conditions include $\tilde u=\tilde v=\tilde w=0$ at $\tilde z=0$, $\tilde u=\tilde v=\partial \tilde w/\partial \tilde r=0$ at $\tilde r=0$, and the outflow condition $\partial p/\partial z = 0$ at $\tilde z=\tilde z_{max}$. The results for the flow in the collision region $\tilde r \sim \tilde z \sim 1$ were found to be independent of the size of the computational domain provided that the boundaries were selected in the ranges $5 \lesssim \tilde r_{max} \lesssim 10$ and $10 \lesssim \tilde z_{max} \lesssim 20$. Illustrative results using $\tilde r_{max}=8$ and $\tilde z_{max}=16$ are shown in figure 9(a–c) for $\sigma =1$ and three different values of the relevant Reynolds number ${\textit {Re}}$. It should be noted that previous experimental results (Phillips Reference Phillips1985) suggest that, for the two largest Reynolds numbers considered, namely, ${\textit {Re}} = 10^5$ and ${\textit {Re}} = 4\times 10^5$, the boundary layer of the steady-flow solutions considered here is probably unstable and would experience transition to a turbulent state, but that aspect of the problem is not investigated in our numerical computations, which are focused instead on the emergence of boundary-layer separation.

Figure 9. The viscous structure of the collision region for $\sigma =1$, including streamlines obtained by integration of (4.1)–(4.4) for (a) ${\textit {Re}} = 10^4$, (b) ${\textit {Re}} = 10^5$ and (c) ${\textit {Re}} = 4\times 10^5$, together with a comparison of radial profiles of (d) axial velocity $\tilde {w}$ and (e) circulation $\tilde {r}\tilde {v}$ at $\tilde {z} = 5$ with the inviscid results shown in figures 7(c) and 7(d). Streamlines in (a–c) are not equispaced. The arrows in (b,c) indicate the location where the boundary layer separates.

It can be seen in figure 9(a) that the structure of the flow for ${\textit {Re}}=10^4$ is very similar to the inviscid structure, in that the boundary layer remains attached and the resulting streamlines are similar to those shown in figure 7(a). By way of contrast, the flow structure found when the Reynolds number is increased to ${\textit {Re}}=10^5$ is markedly different (figure 9b). The streamline pattern reveals the presence of a slender recirculating bubble adjacent to the wall, generated by the separation of the boundary layer at $\tilde r \simeq 1.55$, and subsequent reattachment at $\tilde r \simeq 0.3$. Further increasing the Reynolds number to ${\textit {Re}} = 4\times 10^5$ causes the recirculation bubble to enlarge, and this also moves the location of the point at which the boundary layer separates to $\tilde {r} = 1.75$, approaching the value $\tilde {r} = 1.76$ predicted by the boundary-layer computations.

The closed recirculating bubble has a limited effect on the vertical jet issuing from the collision region. This is quantified in figures 9(d) and 9(e), where profiles of axial velocity and circulation at $\tilde {z} = 5$ for the three values of the Reynolds number considered before are compared with the inviscid results shown in figures 7(c) and 7(d). A noticeable difference is found near the axis, where the sharp peak of the inviscid axial velocity is smoothed by the viscous forces. The resulting near-axis boundary layer is thicker for smaller Reynolds numbers, resulting in a smaller peak velocity. The quantitative agreement everywhere else is quite satisfactory, with the viscous results approaching the inviscid profile for increasing values of the Reynolds number.

It is worth pointing out that closed recirculating bubbles, similar to those shown in figures 9(b) and 9(c), were observed near the end wall in early flow simulations of Ward-type vortex chambers with Reynolds number $10^3$ (based on the flow rate and on the vortex-chamber radius) when the swirl level was sufficiently low (Rotunno Reference Rotunno1979). The separation bubble disappears for increasing swirl level (Rotunno Reference Rotunno1979) and is not present in subsequent computations of the same flow at Reynolds number $10^4$ (Wilson & Rotunno Reference Wilson and Rotunno1986), for which the flow was shown to be fundamentally inviscid. No indication of boundary-layer separation was found in recent tornado simulations at much higher Reynolds numbers (Rotunno et al. Reference Rotunno, Bryan, Nolan and Dahl2016) employing a prescribed forcing term in the vertical momentum equation to generate the motion. The differences between the results of these previous simulations (Rotunno Reference Rotunno1979; Wilson & Rotunno Reference Wilson and Rotunno1986; Rotunno et al. Reference Rotunno, Bryan, Nolan and Dahl2016) and the predictions reported above are attributable to the differences in the associated flow field, suggesting that the detailed distribution of near-wall velocity plays a critical role in the occurrence of boundary-layer separation on approaching the axis.