1. Introduction
Investigations of wall-bounded cyclonic motions continue to serve as focus areas in fluid mechanics and propulsion research. This interest is largely driven by specific characteristics of swirl-assisted combustion chambers. These include intensified mixing, enhanced flame stability and internalised film cooling. Such features directly contribute to improved engine efficiency, extended operational longevity and reduced thermal stresses on components which, taken collectively, allow for more compact and cost-effective designs (Lilley Reference Lilley1977).
Beyond traditional combustors, cyclonic flow patterns, which are inspired by meteorological phenomena such as twisters, dust devils and waterspouts, have also been conceived and implemented into propulsive devices (Barber & Majdalani Reference Barber and Majdalani2017). One such example consists of the Vortex Combustion Cold-Wall (VCCW) thrust chamber, which is responsible for providing lifting, manoeuvring, docking and landing capabilities to the Dream Chaser spaceplane. As described by Chiaverini et al. (Reference Chiaverini, Malecki, Sauer, Knuth, Gramer and Majdalani2003), the VCCW chamber relies on a coaxial, bidirectional swirling motion that is generated by injecting an oxidiser tangentially to the inner circumference and directly upstream of the nozzle’s converging section, where the Mach number can be 0.3 or higher. The ensuing design concept has been shown to achieve high specific impulse across various chamber pressures while maintaining cool wall temperatures by specifically confining propellant mixing and combustion to the chamber core. The so-called ‘cold-wall’ attribute is ascribed in this case to the internal film cooling effect that accompanies the low-temperature oxidiser which, after being introduced tangentially into the chamber, continuously revolves around while sweeping the bounding walls.
In complementing these practical developments, a series of numerical studies may be seen to focus on providing insights into cyclonic flow patterns. For instance, investigations into both hydraulic and gaseous cyclones by Hsieh & Rajamani (Reference Hsieh and Rajamani1991) and Hoekstra, Derksen & Van Den Akker (Reference Hoekstra, Derksen and Van Den Akker1999) help to establish the limitations of conventional eddy-viscosity models at capturing complex swirl phenomena. Along similar lines, Knuth et al. (Reference Knuth, Gramer, Chiaverini, Sauer, Whitesides and Dill1998) undertake finite difference Navier–Stokes simulations to predict the helical flow features in a vortex-injection hybrid rocket engine. The regression-rate characteristics of the latter are also examined by Knuth et al. (Reference Knuth, Chiaverini, Sauer and Gramer2002). Then, using a large-eddy simulation scheme, Derksen & Van den Akker (Reference Derksen and Van den Akker2000) manage to predict vortex-core precession in a cyclonic chamber while closely matching experimental observations at similar Reynolds numbers. In the interim, Slack et al. (Reference Slack, Prasad, Bakker and Boysan2000) implement both Reynolds-stress and large-eddy simulation models on unstructured grids to achieve favourable agreement with experimental measurements across several axial positions in a high-efficiency Stairmand cyclone.
Building upon these computational advancements, Cervantes & Gustavsson (Reference Cervantes and Gustavsson2006) introduce a novel approach for estimating the radial velocity in a swirling field; they also showcase the role of swirl in prescribing the pressure distribution and recovery, particularly in regions comprising recirculation zones. Moreover, parametric trade analyses by Maicke & Talamantes (Reference Maicke and Talamantes2017) help to quantify the manner by which increasing the outlet diameter in a cyclonic chamber can shift the peak tangential velocity radially outward, thus reducing the swirl intensity and potentially producing multiple mantles. Lastly, several sequential studies may be seen to evaluate the effects of varying the geometric and injection parameters on the flow topology of VCCW configurations (Sharma & Majdalani Reference Sharma and Majdalani2021, Reference Sharma and Majdalani2022, Reference Sharma and Majdalani2024).
In seeking to reconcile numerical predictions with empirical observations, analytical models have also been devised for wall-bounded cyclonic motions in a variety of chamber configurations. A common starting point for such models consists of the Bragg–Hawthorne equation (BHE), a specialised reformulation of Euler’s momentum balance for axisymmetric flows. Also known as the Hicks equation (Hicks Reference Hicks1898), the BHE was developed by Bragg & Hawthorne (Reference Bragg and Hawthorne1950) to study actuator disc-induced flows. In actuality, the BHE’s explicit incorporation of specific angular momentum and stagnation properties, expressed in terms of the mean flow streamfunction, makes it particularly well suited for modelling cyclonic flow problems. Corresponding solutions extend over conical, cylindrical and hemispherical chamber geometries, with each entailing its unique challenges and flow characteristics.
In conical domains, for example, several theoretical contributions may be worth noting. These include the asymptotic formulation by Bloor & Ingham (Reference Bloor and Ingham1987), which focuses on industrial cyclones, and the Beltramian solution presented by Barber & Majdalani (Reference Barber and Majdalani2017). Due to its exact nature, the latter remains applicable irrespective of the cone’s divergence half-angle. In the meantime, Grimble & Agarwal (Reference Grimble and Agarwal2015) use dimensional analysis to develop a scaling framework for modelling hydrodynamic instabilities in conical cyclone separators. Their work unravels a simple empirical relation for the Strouhal number as a function of the swirl number and outlet-to-base-diameter ratio. It also connects the hydrodynamically unstable precessing vortex core to the often-reported audible hum in cyclone separators. As a follow-on, Grimble, Agarwal & Juniper (Reference Grimble, Agarwal and Juniper2017) apply linear stability analysis on Reynolds-averaged Navier–Stokes-derived mean flows to systematically characterise these oscillations and predict the onset of precessing vortex core type instabilities. In their essential sequel, they identify absolute instabilities in helical modes that match experimental frequencies and spatial patterns quite favourably, thus leading to deeper insights into the emergence of global oscillations in related experiments.
Unlike the limited number of confined conical profiles (Cortes & Gil Reference Cortes and Gil2007), a much broader spectrum of analytical solutions may be identified in cylindrical configurations. These are typically grouped into inviscid, viscous-corrected and compressibility-corrected categories. For instance, Vyas & Majdalani (Reference Vyas and Majdalani2006) obtain an exact Euler solution by assuming an azimuthal vorticity that varies linearly with the streamfunction. Their approach is further amended to produce Beltramian and Trkalian solutions through modifications of the stagnation head and specific angular momentum relations (Majdalani Reference Majdalani2012). The class of generalised Beltramian solutions is once again extended by Cecil & Majdalani (Reference Cecil and Majdalani2022) to right-cylindrical cyclonic chambers exhibiting hollow cores.
The effort to integrate viscous effects into cyclonic models of cylindrical chambers constitutes another fundamental aspect that may be worth noting. Inspired by the work of Stewartson & Hall (Reference Stewartson and Hall1963) and Mayer & Powell (Reference Mayer and Powell1992), Maicke & Majdalani (Reference Maicke and Majdalani2008) compare two models for the swirl velocity: a Rankine-like piecewise expression and a matched-asymptotic expansion that reproduces the problem’s inner and outer fields in their respective domains. In practice, the latter allows for a smooth transition to occur between the inviscid far field and the forced viscous core (Lewellen Reference Lewellen1962; Stewartson & Hall Reference Stewartson and Hall1963). Similarly, using matched-asymptotic expansions, Majdalani & Chiaverini (Reference Majdalani and Chiaverini2009) incorporate viscous effects into the inviscid solution developed by Vyas & Majdalani (Reference Vyas and Majdalani2006), thereby securing the velocity vanishing requirements at both the wall and centreline regions.
Given that many applications of cyclonic motions occur at relatively high speeds, the incorporation of dilatational effects into basic models of wall-bounded cylindrical chambers has also been attempted. As shown by Maicke, Cecil & Majdalani (Reference Maicke, Cecil and Majdalani2017), a compressibility correction to the inviscid Beltramian class of solutions can be obtained using a multidimensional Rayleigh–Janzen (RJ) expansion. The latter is further extended by Cecil, Little & Majdalani (Reference Cecil, Little and Majdalani2023) to a Trkalian profile. Their findings identify the axial coordinate, normalised by the distance from the chamber head-end to the self-choking sonic point, as a self-preserving, Mach number independent, scaling parameter.
In addition to conical and cylindrical formulations, several models have been developed in spherical domains. One may begin with Hill’s vortex (Hill Reference Hill1894), a spherical ring of constant vorticity translating steadily through an irrotational fluid due to its self-induced velocity. Being recognised as a special non-swirling case of the doubly infinite family of swirling vortices identified by Moffatt (Reference Moffatt1969), Hill’s vortex stability, uniqueness and physical extensions have been the subject of extensive analysis. For example, its compressible form is obtained asymptotically in the form of an RJ expansion by Moore & Pullin (Reference Moore and Pullin1998). Then, using a global spectral stability approach, Protas & Elcrat (Reference Protas and Elcrat2016) complement the work of Moffatt & Moore (Reference Moffatt and Moore1978) by quantifying the inherent instability of Hill’s vortex to axisymmetric deformations; in this process, two pairs of stable and unstable eigenmodes are identified that give rise to infinitely sharp peaks at the sphere’s front and rear stagnation points. Moreover, Scase & Terry (Reference Scase and Terry2018) manage to derive exact solutions for swirling spherical vortices that propagate steadily along the axis of a rotating ideal fluid. Their work shows that Hill’s vortex can be seamlessly matched onto a rotating outer fluid region, wherein concentric ‘sibling’ vortex rings can be engendered once a critical rotation rate (or inverse Rossby number) is exceeded. Along similar lines, Crowe, Kemp & Johnson (Reference Crowe, Kemp and Johnson2021) use numerical simulations and asymptotic expansions in the inverse Rossby number to predict the decay rate of the vortex’s translational speed and radius in the presence of a weakly rotating far field. As for the effects of small viscosity, Smith & Wang (Reference Smith and Wang2021) evaluate whether exact inviscid solutions such as Hill’s spherical vortex, Stuart cats’ eyes, Mallier–Maslowe vortices and the Lamb–Chaplygin dipole represent steady-state or quasi-steady-state approximations to the Navier–Stokes equations in the high Reynolds number limit. Using the Fredholm alternative, they show that some inviscid solutions, such as Hill’s vortex core and the Lamb–Chaplygin dipole, remain consistent at satisfying the problem’s solvability conditions under steady or quasi-steady viscous conditions, respectively. Conversely, Stuart cats’ eyes and Mallier–Maslowe vortices are shown to develop inconsistencies at satisfying the underlying solvability requirements by virtue of the cumulative effects of viscosity over time. The role of Hill’s spherical vortex as the swirl-free limit of the Hicks–Moffatt dipole is also confirmed by Viúdez (Reference Viúdez2022); therein, this classical vortex is recognised as the foundational member of a multipolar family of spherical solutions to the steady vorticity equation.
Besides the foregoing external solutions in spherical geometry, several incompressible models have been developed in confined hemispherical configurations, namely, for the class of vortex engines under development by Sierra Space Corporation. These include a series of studies leading to exact inviscid Beltramian, quasi-complex lamellar and irrotational profiles for hemispherically bounded cyclonic flows (Williams & Majdalani Reference Williams and Majdalani2021a , Reference Williams and Majdalanib ; Majdalani & Williams Reference Majdalani and Williams2021). In contrast, compressible models of confined cyclonic flows in hemispherical chambers, such as those depicted in figure 1, remain untreated. For this reason, the present study devotes itself to the development of a mathematical procedure for capturing the high-speed flow effects in confined cyclonic Beltramian models.
Graphical representations of a compact, hemispherically bounded cyclonic engine concept using (a) side, (b) top and (c) isometric views.
The paper is organised as follows. Initially in § 2, the conservation equations will be reduced to a dimensionless, compressible BHE in spherical coordinates. Appropriate relations between the tangential velocity and the streamfunction will then be established, and the primary flow variables, equations and boundary conditions will be expanded in terms of the squared reference Mach number. In § 3, this expansion will be shown to produce a coupled pair of equations exhibiting both incompressible and compressible scales. While the leading-order solution will be shown to recover the steady-state motion of an ideal incompressible profile, the first-order solution will systematically capture the effects of fluid compression. In § 4, several kinematic and thermodynamic variables will be evaluated and discussed. These include the velocity, pressure, density, temperature, Mach number, vorticity and helicity density distributions. After a brief discussion of the solution behaviour in a fully spherical domain, several conclusions and suggestions for future work will be provided in § 5.
2. Problem formulation
2.1. Fundamental assumptions and geometry
We begin by specifying the domain geometry and adopting a spherical coordinate system that is appropriate for the underlying symmetry. In this process, the position and angular coordinates,
$\bar {R}$
,
$\phi$
and
$\theta$
, are taken to denote the radial, polar and azimuthal directions, respectively. As shown in figure 2(a), the radial coordinate extends from the chamber origin to its outer boundary at
$\bar {R} = a$
, the polar angle sweeps counterclockwise from the vertical axis to the equatorial plane at
$\phi = \pi /2$
and the azimuthal variable,
$0\leqslant \theta \lt 2\pi$
, encircles the domain around the
$z$
axis.
Depiction of (a) right-handed spherical coordinate system that is anchored at the origin
$O$
, the epicentre of a hemispherically bounded chamber. Besides the main variables
$(\bar {R},\phi ,\theta )$
and geometric parameters
$(a,b)$
, we also show in (b) a planar rendering of the chamber dome with clearly identified inflow and outflow regions.
As depicted in figure 2(b), the fluid enters the chamber across the equatorial plane within the radial range of
$b \lt \bar {R} \lt a$
. As the flow traverses the chamber, it reverses its axial direction and exits through an inner vortex tunnel that recrosses the equatorial plane within
$0 \leqslant \bar {R} \lt b$
. Since the radial location
$b$
demarcates the inflow–outflow boundary, its value will depend on the resulting velocity formulation. To facilitate the attendant analysis, the working fluid is modelled as a steady, axisymmetric, reversible, adiabatic and calorically perfect gas. Additionally, gravitational forces and viscous effects are ignored.
2.2. Conservation equations and boundary conditions
2.2.1. Conservation equations
Our equations comprise five conservation statements: one for mass, three for momentum and one for energy. The system is closed using the ideal gas law. In alignment with our initial assumptions, these equations may be reduced to (Batchelor Reference Batchelor2000)
\begin{equation} \bar {{\boldsymbol{\nabla }}}\boldsymbol{\cdot }(\bar {\rho }\bar {\boldsymbol{u}}) = 0, \qquad \bar {\boldsymbol{u}} \boldsymbol{\cdot }\bar {{\boldsymbol{\nabla }}}\bar {\boldsymbol{u}} = -\bar {{\boldsymbol{\nabla }}}\bar {p}/\bar {\rho }, \qquad \bar {{\boldsymbol{\nabla }}}\boldsymbol{\cdot }(\bar {\vphantom {\boldsymbol{u}}\rho }\bar {H} \bar {\boldsymbol{u}}) = 0. \end{equation}
In the above,
$\bar {\rho }$
,
$\bar {\boldsymbol{u}}$
,
$\bar {p}$
and
$\bar {H}$
correspond to the fluid density, velocity, pressure and stagnation enthalpy, respectively. Under the stated assumptions, the stagnation enthalpy can be written as
\begin{equation} \bar {H} = \bar {e} + \bar {p}/\bar {\rho } + \tfrac {1}{2}\bar {\boldsymbol{u}}\boldsymbol{\cdot }\bar {\boldsymbol{u}} = c_{\!p}\bar {T} + \tfrac {1}{2}\bar {\boldsymbol{u}}\boldsymbol{\cdot }\bar {\boldsymbol{u}}, \end{equation}
where
$\bar {e}$
,
$c_{\!p}$
and
$\bar {T}$
stand for the specific internal energy, isobaric specific heat capacity and temperature.
2.2.2. Boundary conditions
Our boundary conditions correspond to an impervious chamber wall, a vanishing cross-flow along the axis of symmetry and a prescribed mass balance in the equatorial plane. These conditions translate into
\begin{align} & \bar {u}_{\!R}(a,\phi )= 0\qquad\!\! &&\text{(no cross-flow along rigid wall),}& \end{align}
\begin{align} & {\lim _{\phi \to 0^{+}}\bar {u}_\phi (\bar {R},\phi )}= 0\qquad\!\! & &\text{(no cross-flow along centreline),}& \end{align}
\begin{align} &{\int _0^b \!2\pi \bar {R}(\bar {\rho } \bar {u}_\phi )\big |_{\phi = \tfrac {1}{2}\pi } \, {\rm d}{\bar {R}}} = \dot {m}= \rho _{{i}} U A_{{i}}& \qquad\!\! &\text{(mass balance along exit plane).}& \end{align}
Here,
$\dot {m}$
represents the injection mass flow rate and
$\rho _{{i}}$
,
$U$
and
$A_{{i}}$
denote the injection density, average velocity and area, respectively. In this context, the incoming stream, which enters the vortex chamber tangentially at an average speed of
$U$
, is allowed to turn axially upon entry to produce the annular updraft. In this work, we follow Bloor & Ingham (Reference Bloor and Ingham1987) by assuming that the three-dimensional development of the tangentially injected stream into an upward-directed motion occurs instantaneously, namely, in a manner to fully occupy the space between
$b$
and the bounding wall. As for the value of
$b$
, which serves as the upper integration bound in (2.3c
), it is implicitly defined by
\begin{equation} \bar {u}_\phi (b,\pi /2) = 0. \end{equation}
Specifically,
$b$
corresponds to the radial position in the equatorial plane where the polar velocity vanishes. Because the spherical
$\bar {u}_\phi$
and cylindrical
$\bar {u}_z$
components coincide in this plane,
$b$
also denotes the radius at which the axial velocity reverses polarity, namely, at the mantle’s intersection with the equatorial plane. In this context, the mantle refers to the axisymmetric interfacial layer that separates the updraft and downdraft regions.
2.3. Normalisation
In what follows, dimensional quantities are denoted by an overbar, while their dimensionless counterparts are left unmarked. To make further headway, it is convenient to introduce a set of dimensionless variables that are consistent with those adopted in previous work (Williams & Majdalani Reference Williams and Majdalani2021a
). In this compressible analysis, however, we find it more beneficial to normalise the stagnation enthalpy,
$\bar {H}$
, by
$c_{\!p} T_0$
in lieu of
$U^2$
. The normalised variables can be established by rescaling the dimensional (barred) quantities using their respective reference values. This may be accomplished by setting
\begin{gather*} R = \bar {R}/a, \quad \beta = b/a, \!\quad {\boldsymbol{\nabla }} = a\bar {{\boldsymbol{\nabla }}}, \!\quad \boldsymbol{u} = \bar {\boldsymbol{u}}/U, \!\quad \boldsymbol{\omega } = a\bar {\boldsymbol{\omega }}/U, \!\quad \psi = \bar {\psi }/\bigl ( \rho _0 U a^2 \bigr ), \nonumber \\ H = \bar {H}/(c_{\!p} T_0), \!\quad B = \bar {B}/(Ua), \!\quad p = \bar {p}/p_0, \!\quad \rho = \bar {\rho }/\rho _0, \!\quad T= \bar {T}/T_0 \!\quad \text{and} \!\quad S = \bar {S}/c_{\!p}, \end{gather*}
where
$a$
,
$\bar {\boldsymbol{\omega }}$
,
$\bar {\psi }$
,
$\bar {B}$
and
$\bar {S}$
denote the chamber radius, vorticity, streamfunction, specific angular momentum and entropy, respectively. The pressure, density and temperature are normalised by their stagnation values that may be taken at the apex of the chamber, specifically by
$p_0$
,
$\rho _0$
and
$T_0$
. By implementing these variables, the dimensional equations of motion from (2.1) simplify into
\begin{equation} {\boldsymbol{\nabla }} \boldsymbol{\cdot }(\rho \boldsymbol{u}) = 0, \qquad \boldsymbol{u} \boldsymbol{\cdot }{\boldsymbol{\nabla }} \boldsymbol{u} = -{\boldsymbol{\nabla }}\! p / ( \gamma M_0^2 \rho ) \qquad \text{and} \qquad {\boldsymbol{\nabla }} \boldsymbol{\cdot }(\rho H \boldsymbol{u}) = 0, \end{equation}
where the pressure, density and temperature remain connected through the isentropic relations,
$p = \rho ^\gamma = T^{\gamma /(\gamma - 1)}$
, with
$\gamma \equiv c_{\!p}/c_v$
denoting the ratio of specific heats. Moreover,
$M_0 \equiv U/c_0$
specifies the reference Mach number, where
$c_0 = (\gamma \mathcal{R} T_0)^{1/2}$
determines the stagnation point acoustic speed, and
$\mathcal{R}$
stands for the specific gas constant. After substituting these transformations into (2.2), the normalised stagnation enthalpy becomes
\begin{equation} H = \rho ^{\gamma - 1} + \tfrac {1}{2}(\gamma - 1)M_0^2\,\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}. \end{equation}
In the same vein, the dimensional boundary conditions from (2.3a ), (2.3b ) and (2.3c ) may be readily converted into their dimensionless forms by writing
\begin{align} & u_{\!R}(1,\phi ) = 0 &\text{(no cross-flow along rigid wall),} \end{align}
\begin{align} &\lim _{\phi \rightarrow 0^{+}}u_\phi (R,\phi ) = 0 & \text{(no cross-flow along centreline),} \end{align}
\begin{align} & \int _0^\beta \!(\rho u_\phi )|_{\phi = \tfrac {1}{2}\pi }\, R\,{\rm d}{R} = \kappa & \text{(mass balance along exit plane),} \end{align}
where the Ekman-type scaling parameter,
$\kappa \equiv \rho _{{i}} A_{{i}}/(2\pi \rho _0 a^2)$
, is recovered. Often referred to as the ‘inflow parameter’,
$\kappa$
gauges the ratio of two characteristic mass-flux surfaces: the fluid injection surface,
$\rho _{{i}} A_{{i}}$
, and the chamber cross-sectional area in the equatorial plane,
$2\pi \rho _0 a^2$
.
Lastly, the condition in (2.4) may be restated in dimensionless form using
\begin{equation} \mathinner {u_\phi (\beta ,\pi /2)} = 0. \end{equation}
At this juncture, by leveraging the rotational symmetry about the
$z$
axis, we may proceed by reducing the system of equations given by (2.6) to the compressible form of the BHE in terms of the streamfunction,
$\psi$
. As usual, the latter may be related to the spherical radial and polar velocities through
\begin{equation} u_{\!R} = -\frac {1}{R^2\sin {\phi }\,\rho }\frac {\partial \psi }{\partial \phi } \qquad \text{and} \qquad u_\phi = \frac {1}{R\sin {\phi }\,\rho }\frac {\partial \psi }{\partial R}. \end{equation}
For the reader’s convenience, the underlying formulation is outlined in Appendix A.
2.4. Density–streamfunction formulation
2.4.1. Compressible Bragg–Hawthorne equation
Having specified the dimensionless equations of motion and associated boundary conditions, we may begin with the compressible BHE:
\begin{align} \textrm{L}^2 \psi + \rho ^2 \biggl (\! B\frac {{\rm d} B}{{\rm d} \psi } - R^2 \sin ^2{\!\phi }\,\,\frac {{\rm d} H}{{\rm d} \psi } \biggr ) = {\boldsymbol{\nabla }}\psi \boldsymbol{\cdot }{\boldsymbol{\nabla }}(\ln \rho ); \!\quad \textrm{L}^2 \equiv \frac {\partial ^2}{\partial R} + \frac {\sin {\phi }}{R^2}\frac {\partial }{\partial \phi }\biggl ( \frac {1}{\sin {\phi }}\frac {\partial }{\partial \phi } \biggr ), \end{align}
where
$B(\psi ) = R\sin {\phi } \mathinner {u_\theta }$
and
$\textrm{L}^2 $
stands for the Stokes operator. To achieve system closure, an additional relation between
$\rho$
and
$\psi$
is required. The stagnation enthalpy, given by (2.7), can be recast in terms of
$\psi$
to produce
\begin{equation} \mathinner {H{(\psi)}} = \rho ^{\gamma - 1} + \frac {\gamma - 1}{2}\biggl [ \frac {{\boldsymbol{\nabla }}\psi \boldsymbol{\cdot }{\boldsymbol{\nabla }}\psi + \rho ^2\,B^2(\psi )}{\rho ^2 R^2\sin ^2\!\phi } \biggr ]M_0^2. \end{equation}
By prescribing a functional dependence for
$B$
and
$H$
on
$\psi$
, the combination of (2.11) and (2.12) forms a complete set that can serve as the basis for analysis. In this effort, however, and irrespective of the particular functions
$B$
and
$H$
, one arrives at a nonlinear relation between
$\psi$
and
$\rho$
. The ensuing nonlinearity necessitates the use of approximation techniques as it precludes the extraction of exact analytical solutions.
2.4.2. Boundary condition reformulation
Having expressed (2.11) and (2.12) in terms of
$\psi$
and
$\rho$
, we now extend the treatment to the boundary conditions given by (2.8a
), (2.8b
) and (2.8c
). To facilitate this effort, it is helpful to first examine the possible domain and range of
$\psi$
and
$\rho$
vis-à-vis potential singularities. For example, the presence of either
$R$
or
$\sin \phi$
in the denominator of expressions such as (2.11) and (2.12) or (A2) gives rise to singularities at the centreline. In avoidance of such behaviour, the boundary conditions given by (2.8a
), (2.8b
) and (2.8c
) are evaluated using one-sided limits. To do so, we set
\begin{align} & \!\left ( \frac {1}{\rho }\frac {\partial \psi }{\partial \phi } \right )\Big |_{R = 1} = 0 & \text{(no cross-flow along rigid wall),} \end{align}
\begin{align} & \lim _{\phi \rightarrow 0^{+}} \,\!\left ( \frac {1}{\rho \sin {\phi }}\frac {\partial \psi }{\partial R} \right ) = 0 & \text{(no cross-flow along centreline),} \end{align}
\begin{align} & \int _0^\beta \!\frac {\partial \psi }{\partial R}\Big |_{\phi =\tfrac {1}{2}\pi } \, {\rm d}{R} = \kappa & \text{(mass balance along exit plane).} \end{align}
Furthermore, the requirement on the radial fraction
$\beta$
, given by (2.9), can be restated in terms of
$\psi$
by writing
\begin{equation} \frac {\partial \psi }{\partial R}\!\left ( \beta , \tfrac {1}{2}\pi \right ) = 0. \end{equation}
Turning once again to the boundary conditions at hand, it should be noted that, due to the physical constraint of a positive and finite fluid density, the limit property
\begin{equation} \lim _{\boldsymbol{x} \rightarrow \boldsymbol{x}_0} \frac {\mathinner {f(\boldsymbol{x})}}{\mathinner {\rho {\left ( \boldsymbol{x} \right )}}} = \frac {\lim _{\boldsymbol{x}\rightarrow \boldsymbol{x}_0} \mathinner {f{\left ( \boldsymbol{x} \right )}}}{\lim _{\boldsymbol{x}\rightarrow \boldsymbol{x}_0} \mathinner {\rho {\left ( \boldsymbol{x} \right )}}} \end{equation}
applies to any function
$\mathinner {f(\boldsymbol{x})}$
that is well defined at
$\boldsymbol{x}_0$
. As a result,
$\rho$
may be eliminated from (2.13a
) and (2.13b
) to the extent of reducing (2.13a
), (2.13b
) and (2.13c
) into
\begin{align} & \frac {\partial \psi }{\partial \phi }(1,\phi ) = 0 & \text{(no cross-flow along rigid wall),} \end{align}
\begin{align} & \lim _{\phi \rightarrow 0^{+}}\frac {1}{\sin \phi }\frac {\partial \psi }{\partial R} = 0 & \text{(no cross-flow along centreline),} \end{align}
\begin{align} & \psi (\beta ,\pi /2) = \kappa + \lim _{R \rightarrow 0^{+}}\mathinner {\psi (R,\pi /2)} & \text{(mass balance along exit plane).} \end{align}
Equations (2.16a ), (2.16b ) and (2.16c ) are used hereafter due to their equivalence to (2.13a ), (2.13b ) and (2.13c ). In the absence of singularities, these constraints hold for both compressible and incompressible formulations.
2.4.3. Specification of stagnation enthalpy and circulation
Besides assuming the flow to be isentropic, we take the stagnation enthalpy,
$H$
, to be uniform throughout the domain. It may be instructive to note that for a steady, inviscid flow, such a condition, termed homenergic, is not a consequence of isentropic motion alone. The relation between these quantities may be deduced from Crocco’s theorem (see § 3.5 in Batchelor Reference Batchelor2000). Accordingly, for steady, inviscid, calorically perfect gas, one can put
\begin{equation} \bar {\boldsymbol{u}}\times \bar {\boldsymbol{\omega }} =\bar {\boldsymbol{\nabla }} \kern1pt \bar {H} \kern1pt -\bar {T} \kern1pt \bar {\boldsymbol{\nabla }} \kern1pt \bar {S} \qquad \text{or} \qquad (\gamma -1)M_0^2 {\boldsymbol{u}}\times {\boldsymbol{\omega }} ={\boldsymbol{\nabla }} {H} -{T}{\boldsymbol{\nabla }} {S}. \end{equation}
When the isentropic flow originates from a region of uniform energy, as in the case here, it becomes homentropic, i.e. with both
$\boldsymbol{\nabla }{S} = {0}$
and
$\boldsymbol{\nabla }{H} = {0}$
. Inserting these into (2.17), one is left with a vanishing Lamb vector,
$\boldsymbol{u} \times \boldsymbol{\omega } = {0}.$
We hence recover the defining criterion for a Beltrami field, namely, one for which the velocity and vorticity vectors remain parallel everywhere. This condition enables us to further write
\begin{equation} k\boldsymbol{u} = {\boldsymbol{\nabla }}\times \boldsymbol{u}, \end{equation}
where
$k$
denotes a scalar function that may vary spatially. For the spatially invariant case of
$\boldsymbol{\nabla }k = {0}$
, the helical field reduces to a Trkalian motion.
As for the constant value of
$H$
, it may be evaluated at the reference stagnation state that is accompanied by a vanishing velocity and a temperature of
$T_0$
. This enables us to specify the dimensional stagnation enthalpy as
$\bar {H} = c_{\!p} T_0$
everywhere. Then, since
$H = \bar {H}/(c_{\!p} T_0)$
, we arrive at
$H=1$
.
Turning now to
$B(\psi )$
, and in an effort to recover analytically tractable solutions (Maicke et al. Reference Maicke, Cecil and Majdalani2017), we consider functional forms that are linear in
$\psi$
. By setting
$B(\psi ) {\rm d}B/{\rm d}\psi = C_{B}^2\psi$
, direct integration yields
$B^2 /2 = C_{B}^2 \psi ^2 /2 + K$
, where
$K$
denotes an arbitrary constant. For notational convenience, we simply multiply by 2 and set
$C_0^2 \equiv 2K$
; this enables us to proceed using
\begin{equation} B^2(\psi ) = C_{B}^2 \psi ^2 + C_{0}^2. \end{equation}
In the above, the addition of a non-negative constant,
$C_{0}^2$
, takes into account the physical requirement of pursuing strictly positive real values of
$B^2$
. Recalling the relation between
$B$
and
$u_\theta$
, the tangential velocity becomes
\begin{equation} u_{\theta } = \pm \frac {1}{R\sin {\phi }}\big ( C_{B{}}^2\psi ^2 + C_{0}^2 \big )^{1/2}. \end{equation}
Since the sign above determines the direction of swirl, the case of
$u_{\theta } \geqslant 0$
is selected. At this juncture, it may be helpful to briefly examine the physical significance of
$C_B$
and
$C_0$
, the two constants that appear in (2.20). In the particular case of
$C_B=0$
, the reduced form of
$C_0 = \pm R \sin {\phi } \mathinner {u_\theta } = \pm B$
prescribes the circulation of an irrotational vortex for which
$u_\theta \propto 1/r$
. In this instance, the term ‘irrotational’ refers to a swirling motion with a vanishing axial vorticity, namely, with
$\omega _z = r^{-1}{\partial }(ru_\theta )/{\partial } r = 0$
.
If, on the other hand, we set
$C_0=0$
, we obtain a rotational vortex. In the absence of viscosity, however, only the latter case leads to physically meaningful (non-negative) thermodynamic variables (Maicke et al. Reference Maicke, Cecil and Majdalani2017). Motivated by this condition, we take
$C_0=0$
in the remainder of this work. The full justification and implications of this choice are detailed in the velocity field analysis of § 4.3.4.
With
$B(\psi )$
and
$H(\psi )$
being fully specified, the BHE and stagnation enthalpy expressions in (2.11) and (2.12) can be consolidated into a pair of coupled relations. These are
\begin{align} &\textrm{L}^2\psi + C_B^2 \rho ^2 \psi = {\boldsymbol{\nabla }}\psi \boldsymbol{\cdot }{\boldsymbol{\nabla }}(\ln {\rho }) & \text{(compressible BHE),}& \end{align}
\begin{align} & \rho ^{\gamma - 1} = 1 - \frac {\gamma - 1}{2} M_0^2 \frac {{\boldsymbol{\nabla }}\psi \boldsymbol{\cdot }{\boldsymbol{\nabla }}\psi + C_B^2 \rho ^2 \psi ^2}{\rho ^2 R^2 \sin ^2{\phi }} \qquad & \text{(stagnation enthalpy).}& \end{align}
By inspection of (2.21b ), one may infer that viable solutions will require additional boundedness conditions to ensure that the velocity remains finite throughout the entire domain. Using asymptotic notation, one must have
\begin{equation} {\mathopen \| {\boldsymbol{\nabla }}\psi \mathclose \| = {O}(R \sin {\phi }) \quad \text{and} \quad \psi = o(R \sin {\phi }) \quad \text{as} \quad R\sin {\phi } \to 0.} \end{equation}
By adopting (2.22) as an auxiliary constraint, one ensures that the velocity as well as the main thermodynamic variables remain bounded.
3. Solution procedure
3.1. Asymptotic expansion approach
3.1.1. Rayleigh–Janzen technique
In accordance with the RJ technique, all dependent variables can be expanded as a power series in the squared reference Mach number,
$M_0^2 \ll 1$
. This procedure transforms the nonlinear, coupled density–streamfunction system into a sequence of tractable linear equations that can be solved consecutively. At leading order, this expansion serves a confirmatory role by recovering the same streamfunction of the known incompressible solution to this problem (Williams & Majdalani Reference Williams and Majdalani2021a
). As for the effects of fluid compressibility, they are systematically captured by subsequent terms. In this manner, the RJ expansions for the velocity components,
$(u_{\!R},u_{\phi }, u_{\theta })$
, thermodynamic variables,
$(\rho , p, T)$
, streamfunction,
$\psi$
, vorticity,
$\boldsymbol{\omega }$
, and specific angular momentum,
$B$
, become
\begin{align} \begin{alignedat}{5} &u_{\!R} \sim u_{\!R}^{\scriptscriptstyle {(0)}} + M_0^2 u_{\!R}^{\scriptscriptstyle {(1)}}, \qquad & & u_\phi \sim u_\phi ^{\scriptscriptstyle {(0)}} + M_0^2 u_\phi ^{\scriptscriptstyle {(1)}}, \qquad &&u_{\theta } \sim u_\theta ^{\scriptscriptstyle {(0)}} + M_0^2 u_\theta ^{\scriptscriptstyle {(1)}}, \\ & \rho \sim 1 + M_0^2\rho _1 + M_0^4\rho _2, \qquad & & p \sim 1 + M_0^2p_1 + M_0^4 p_2, \qquad & & T \sim 1 + M_0^2T_1 + M_0^4 T_2, \\ & \psi \sim \psi _0 + M_0^2\psi _1, \qquad & & \boldsymbol{\omega } \sim \boldsymbol{\omega }_0 + M_0^2\boldsymbol{\omega }_1, \qquad & & B \sim B_0 + M_0^2B_1. \end{alignedat} \end{align}
In (3.1), the original variables depend on
$(R,\phi ,M_{0}^2)$
, whereas the expanded terms depend only on
$(R,\phi )$
. It may be helpful to remark that the zeroth-order unit terms that lead the thermodynamic expansions are normalised by their uniform stagnation values. For this reason, the fundamental variations that accompany
$\psi _0$
are reflected in the seemingly first-order corrections,
$\rho _1$
,
$p_1$
and
$T_1$
. As such, one must wait for the second-order terms,
$\rho _2$
,
$p_2$
and
$T_2$
, to actually capture the dilatational effects stemming from
$\psi _1$
. When these expansions are substituted into (2.21), terms of like powers in
$M_0$
may be grouped and segregated into leading- and first-order equations. For composite functions such as
$\mathinner {B(\psi )}$
, the expansion can be reproduced using single-variable differentiation:
\begin{align} \mathinner {B(\psi )} & = \mathinner {B(\psi _0)} + M_{0}^2 \psi _1 \!\mathinner {B'\!(\psi _0)} + M_{0}^4\!\left [ \psi _2 B'\! (\psi _0) + \tfrac {1}{2}\psi _1^2 B''\! (\psi _0) \right ] + \cdots \end{align}
\begin{align} &= B_0 + M_{0}^2 B_1 + M_{0}^4 B_2 + \cdots , \end{align}
where primes denote differentiation with respect to
$\psi _0$
. Comparing terms of like order in
$M_0$
, the first three coefficients,
$B_0$
,
$B_1$
and
$B_2$
, may be readily retrieved. One gets
\begin{equation} B_0 = \mathinner {B(\psi _0)}, \qquad B_1 = \psi _1\!\mathinner {B'\!(\psi _0)} \qquad \text{and}\qquad B_2 = \psi _2\!\mathinner {B'\!(\psi _0)} + \tfrac {1}{2} \psi _1^2\!\mathinner {B''\!(\psi _0)}. \end{equation}
This procedure may be systematically extended to evaluate higher-order terms to any targeted degree of accuracy.
It may be instructive to note that the RJ methodology applies equally well to special parameters that appear in the foregoing equations. In particular,
$C_B$
, first introduced through (2.19), will be shown to play the notable role of an eigenvalue in the forthcoming analysis. In fact, the corresponding linearisation and variable expansions confirm that
$C_B$
must be uniquely determined at each order. To properly allocate its individual contributions across successive perturbative orders, we expand
$C_B$
using
\begin{equation} \mathinner {C_B\bigl (M_{0}^2\bigr )} = C_B^{\scriptscriptstyle {(0)}} + M_{0}^2 C_B^{\scriptscriptstyle {(1)}} + \cdots .\end{equation}
Along similar lines, it may be readily shown that the asymptotic expansion of
$\psi$
, which modifies (2.14), requires that
$\beta$
, the open fraction of the outlet radius, be incrementally reconstructed at each order. Although
$\beta$
is routinely treated as a geometry-dependent constant under incompressible conditions, it becomes dependent on the Mach number in the presence of flow dilatation. To adequately account for the sensitivity of
$\beta$
to
$M_0$
, we let
\begin{equation} \mathinner {\beta \bigl (M_{0}^2\bigr )} = \beta _0 + M_{0}^2 \, \beta _1 + \cdots .\end{equation}
In view of this expansion, the effect of compressibility on the ideal outlet radius,
$\beta$
, can be accurately quantified. By substituting the expanded forms of
$\psi$
and
$\beta$
into (2.16c
) and (2.14), one obtains, for
$R=\beta _0$
and
$\phi =\pi /2$
,
\begin{align} \frac {\partial }{\partial R}\!\!\left [ \!\psi _0 + M_{0}^2\bigl ( \beta _1 {\mathop {\partial _{R}}}\psi _0 + \psi _1 \!\bigr ) +\! M_{0}^4\!\left ( \!\tfrac {1}{2} \beta _1^2 {\mathop {\partial _{RR}}}\psi _0 + \beta _2 {\mathop {\partial _{R}}}\psi _0 + \beta _1{\mathop {\partial _{R}}}\psi _1 + \psi _2 \!\right )\! + \cdots \right ] = 0,& \end{align}
\begin{align} {\psi _0 - \kappa + M_{0}^2\,\bigl ( \beta _1 {\mathop {\partial _{R}}}\psi _0 + \psi _1 \!\bigr ) + M_{0}^4\!\left ( \!\tfrac {1}{2} \beta _1^2 \,{\mathop {\partial _{RR}}}\psi _0 + \beta _2 {\mathop {\partial _{R}}}\psi _0 + \beta _1{\mathop {\partial _{R}}}\psi _1 + \psi _2 \!\right )\! + \cdots } = 0,& \end{align}
where
${\mathop {\partial _{R}}}f \equiv \partial f/\partial R$
. From these expressions, one may recover the leading-order value of the streamfunction at
$\beta _0$
as well as the Mach number correction to
$\beta _0$
. The latter represents the shift in the ideal outlet radius or mantle location in the exit plane. Forthwith, by evaluating (3.6) at each order, one systematically retrieves at
$R = \beta _0$
and
$\phi = \pi /2$
\begin{equation} \begin{aligned} {O}(1)\colon & \quad\psi _0 = \kappa , \quad {\mathop {\partial _{R}}}\psi _0 = 0, \\ {O}(\textstyle M_0^2)\colon & \quad\psi _1 = 0, \quad \beta _1= -{\mathop {\partial _{R}}}\psi _1 \big / {\mathop {\partial _{RR}}}\psi _0. \end{aligned} \end{equation}
To the authors’ knowledge, the effort leading to (3.7) constitutes the first attempt to account for the shifting of the cyclonic mantle due to flow dilatation (Maicke et al. Reference Maicke, Cecil and Majdalani2017). Moreover, although
$\beta _1 \neq 0$
, its effect on the streamfunction only appears at
${O}(M_0^4)$
. Nonetheless, radial derivatives of
$\psi$
remain non-zero at all subsequent orders. Consequently, only the velocity components,
$u_{\!R}$
and
$u_\phi$
, exhibit first-order variations with
$M_0$
on the ring defined by
$(R,\phi ) = (\,\beta _0,\pi /2)$
. This behaviour may be attributed to the dependence of
$u_{\!R}$
and
$u_\phi$
on derivatives of both
$\psi$
and
$\rho$
, unlike
$u_{\theta }$
, which depends solely on
$\psi$
. At the outset, compression-induced variations in
$u_\theta$
along the outlet ring appear as higher-order contributions.
3.1.2. Two-term expansion
Using (3.1), the expansions of (2.21a ) and (2.21b ) enable us to separate the leading- and first-order equations. We get
\begin{align} \hspace {-6pt} &{O}{(1)}\colon \hspace {-4pt}& &\begin{cases} \textrm{L}^2 \psi _0 + \bigl ( C_B^{\scriptscriptstyle {(0)}} \bigr ){}^2 \psi _0 = 0, \\[5pt] \rho _1 = -\dfrac {1}{2}\dfrac {{\boldsymbol{\nabla }}\psi _0 \boldsymbol{\cdot }{\boldsymbol{\nabla }} \psi _0 + \bigl ( C_B^{\scriptscriptstyle {(0)}} \bigr ){}^2 \psi _0^2}{{\vphantom {\sin }R}^2 \sin ^2{\phi }}, \end{cases} \end{align}
\begin{align} \hspace {-6pt} &{O}\bigl ( M_{0}^2 \bigr )\colon \hspace {-4pt}& &\begin{cases} \textrm{L}^2 \psi _1 + \bigl ( C_B^{\scriptscriptstyle {(0)}} \bigr ){}^2 \psi _1 = {\boldsymbol{\nabla }}\! \rho _1\boldsymbol{\cdot }{\boldsymbol{\nabla }} \psi _0 - 2 {C_B^{\scriptscriptstyle {(0)}}}\bigl ( C_B^{\scriptscriptstyle {(1)}} + C_B^{\scriptscriptstyle {(0)}} \rho _1 \bigr )\psi _0, \\[5pt] \rho _2 = \dfrac {2 - \gamma }{2}\rho _1^2 - \dfrac {{\boldsymbol{\nabla }} \psi _0\boldsymbol{\cdot }({\boldsymbol{\nabla }} \psi _1 - \rho _1 {\boldsymbol{\nabla }} \psi _0) + C_B^{\scriptscriptstyle {(0)}} \psi _0\,\bigl(C_B^{\scriptscriptstyle {(0)}}\psi _1 + C_B^{\scriptscriptstyle {(1)}} \psi _0\bigr)}{{\vphantom {\sin }R}^2 \sin ^2{\phi }}. \end{cases} \end{align}
In (3.8a
), the equation for the density correction,
$\rho _1$
, appears as part of the
${O}(1)$
system because it remains solely prescribed by the leading-order streamfunction,
$\psi _0$
. Conversely,
$\rho _2$
may be seen to depend on
$\rho _1$
as well as on the interactions between
$\psi _0$
and
$\psi _1$
. A fortuitous outcome of this expansion may be the decoupling of
$\psi$
and
$\rho$
at each perturbative order. Because the density components
$\rho _1$
and
$\rho _2$
rely respectively on
$\psi _0$
and
$\psi _1$
, the characteristic equations for
$\psi _1$
and
$\psi _2$
may be directly expressed in terms of lower-order streamfunctions. As such, the density components can be conveniently eliminated from (3.8a
) and (3.8b
) to retrieve
\begin{align} {O} (1)\colon \quad & \textrm{L}^2 \psi _0 + \bigl ( C_B^{\scriptscriptstyle {(0)}} \bigr ){}^2\psi _0 = 0, \end{align}
\begin{align} {O}(M_{0}^2)\colon \quad & \textrm{L}^2 \psi _1 + \bigl ( C_B^{\scriptscriptstyle {(0)}} \bigr ){}^2\psi _1 = -2 C_B^{\scriptscriptstyle {(0)}} C_B^{\scriptscriptstyle {(1)}} \psi _0 + \biggl [ \bigl ( C_B^{\scriptscriptstyle {(0)}} \bigr ){}^2 \psi _0 - \frac {1}{2}{\boldsymbol{\nabla }}\psi _0\boldsymbol{\cdot }{\boldsymbol{\nabla }} \biggr ]\notag \\&\hspace {93pt}\times \biggl [ \frac {{\boldsymbol{\nabla }}\psi _0\boldsymbol{\cdot }{\boldsymbol{\nabla }}\psi _0 + \bigl ( C_B^{\scriptscriptstyle {(0)}} \bigr ){}^2 \psi _0^2}{{R\vphantom {\sin }}^2\sin ^2\phi } \biggr ]. \end{align}
The expansion of the remaining thermodynamic properties,
$p$
and
$T$
, may be recovered by combining the definition of stagnation enthalpy with the adiabatic flow relations. One obtains
\begin{align} &p =\rho ^\gamma \!: \; M_{0}^2 (p_1 - \gamma \rho _1) + M_{0}^4\!\left \{ p_2 - \gamma \!\left [ \rho _2 - \tfrac {1}{2}(\gamma - 1)\rho _1^2 \right ] \right \} + \cdots = 0, \end{align}
\begin{align} &T =\rho ^{\gamma -1}\!:\; M_{0}^2 \bigl [ T_1 - (\gamma -1)\,\rho _1 \bigr ] + M_{0}^4 \!\left \{ T_2 - (\gamma - 1)\!\left [ \rho _2 + \tfrac {1}{2}(\gamma - 2)\rho _1^2 \right ] \right \} + \cdots = 0. \hspace {26pt} \\[9pt] \nonumber \end{align}
These enable us to extract the first- and second-order pressure and temperature corrections, specifically
\begin{align}& p_1 = \gamma \rho _1, \quad T_1 = (\gamma - 1)\rho _1, \quad p_2 =\textstyle \gamma \!\left [ \rho _2 + \tfrac {1}{2}(2 - \gamma )\rho _1^2 \right ]\!, \nonumber\\& T_2 = (\gamma - 1) \!\left [ \rho _2 - \tfrac {1}{2} (2 - \gamma )\rho _1^2 \right ]\!. \end{align}
3.1.3. Two-term expansion of boundary conditions
Having linearised the principal equations, we may shift our focus to the boundary conditions in (2.16a ), (2.16b ) and (2.16c ). The underlying requirements may be treated analogously to (2.11) and (2.12). Assuming a smooth solution across the domain, the RJ expansions outlined in (3.1) may be substituted into (2.16a ) , (2.16b ) and (2.16c ) to produce
\begin{align} & \mathinner {\frac {\partial \psi }{\partial \phi }(1,\phi )} = 0 & \text{or}\qquad & && \begin{cases} {O}{(1)}\colon \quad & \mathinner {{\mathop {\partial _{\phi }}}\psi _0(1,\phi )} = 0, \\[5pt] {O}{(M_{0}^2)}\colon \quad & \mathinner {{\mathop {\partial _{\phi }}}\psi _1(1,\phi )} = 0, \\ \end{cases} \end{align}
\begin{align} & \lim _{\phi \to 0^+}\frac {1}{\sin \phi }\frac {\partial \psi }{\partial R} = 0 & \text{or}\qquad & & & \begin{cases} \displaystyle {O}{(1)}\colon \quad & {\lim _{\phi \to 0^+}}\csc {\phi }\,{\mathop {\partial _{R}}}\psi _0 = 0, \\[5pt] \displaystyle {O}{(M_{0}^2)}\colon \quad & {\lim _{\phi \to 0^+}}\csc {\phi }\,{\mathop {\partial _{R}}}\psi _1 = 0, \\ \end{cases} \end{align}
\begin{align} & \psi (\beta ,\pi /2) - \kappa = 0 & \text{or}\qquad & & & \begin{cases} {O}{(1)}\colon \quad & \psi _0(\beta _0,\pi /2) = \kappa , \\[5pt] {O}{(M_{0}^2)}\colon \quad & \mathinner {\psi _1(\beta _0,\pi /2)} = 0. \\ \end{cases} \end{align}
In (3.12c
), the RJ expansion of
$\beta$
is effectively applied to rewrite the condition entirely in terms of
$\psi _m$
and
$\beta _0$
, where
$m=\{0,1,\ldots \}$
.
3.2. Compressible Bragg–Hawthorne equation treatment
With the availability of the compressible flow equations and boundary conditions at each asymptotic order, the solution can be constructed sequentially. First, the zeroth-order variables can be determined. Then, the procedure leading to the first-order compressibility corrections can be strategically implemented.
3.2.1. Leading-order solution
At leading order, (3.9a
) along with the requirements encapsulated through (3.12a
), (3.12b
) and (3.12c
) must be resolved. The problem reduces to the specification of a function
$\psi _0$
that satisfies the partial differential equation (PDE) given by
\begin{equation} \frac {\partial ^2 \psi _0}{\partial R^2} + \frac {\sin \phi }{R^2}\frac {\partial }{\partial \phi }\biggl ( \frac {1}{\sin {\phi }} \frac {\partial \psi _0}{\partial \phi } \biggr ) + \bigl ( C_B^{\scriptscriptstyle {(0)}} \bigr ){}^2\psi _0 = 0. \end{equation}
Equation (3.13) must be solved in conjunction with
\begin{align} & \mathinner {\frac {\partial \psi _0}{\partial \phi }{ \left ( 1,\phi \right )}} = 0 & \hspace {1in} & \text{(no cross-flow along rigid wall),}& \end{align}
\begin{align} &{\lim _{\phi \to 0^+}\frac {1}{\sin {\phi }}\frac {\partial \psi _0}{\partial R}} {}= 0 & & \text{(no cross-flow along centreline),}& \end{align}
\begin{align} & \psi _0(\,\beta _0,\pi /2) = \kappa & & \text{(mass balance along exit plane),}& \end{align}
\begin{align} &{\lim _{\substack {R \to 0^{+} \\ \phi \to 0^{+}} }\frac {\psi _0}{R\sin {\phi }}} = 0 & & \text{(bounded flow profile at origin).}& \end{align}
It may be instructive to note that (3.14d
) promotes the development of a physically meaningful solution by suppressing possible singularities at the origin, namely, at the epicentre of the hemispherical enclosure in figure 2. Upon multiplying (3.13) by
$R^2$
, one may infer that the solution,
$\psi _0(R, \phi )$
, can be expressed as the product of two independent functions of the form
\begin{equation} \mathinner {\psi _0(R,\phi )} = f_0(R)g_0(\phi ). \end{equation}
Then, by introducing a positive semidefinite separation constant,
$\upsilon ^2 \geqslant 0$
, a pair of ordinary differential equations (ODEs) and auxiliary conditions may be identified. These are
\begin{align} &f_0'' + \bigl [ \bigl ( C_B^{\scriptscriptstyle {(0)}} \bigr ){}^2 - \upsilon ^2R^{-2} \bigr ]f_0 = 0; & f_0(1) = 0, \quad \lim _{R \to 0^{+}} \bigl [ f_0(R)/R \bigr ] = 0, \end{align}
\begin{align} & g_0'' - \cot {\phi }\,g_0' + \upsilon ^2g_0 = 0; & \lim _{\phi \rightarrow 0^+} \bigl [ g_0(\phi )/\sin {\phi } \bigr ] = 0. \end{align}
In the above, the stated conditions on
$f_0$
and
$g_0$
ensure that the solution remains regular at the origin as
$R \to 0^{+}$
irrespective of
$\phi$
. Additionally, the inflow–outflow mass balance enables us to connect
$f_0$
and
$g_0$
through
\begin{equation} \mathinner {f_0(\beta _0)} \mathinner {g_0(\pi /2)} = \kappa . \end{equation}
Through these manipulations, (3.16) and (3.17) enable us to fully resolve the separated eigenvalue problem. While the algebraic constraints on
$f_0$
correspond to a rigid hemispherical wall, those on
$g_0$
prevent flow crossing along the centreline. In practice, the solution for
$f_0$
may be expressed using Bessel functions and written as
\begin{equation} f_0(R) = R^{1/2} \,\bigl [ C_1 J_\nu \bigl ( C_B^{\scriptscriptstyle {(0)}} R \bigr ) + C_2 Y_{\nu }\bigl ( C_B^{\scriptscriptstyle {(0)}} R \bigr ) \bigr ]; \quad \nu = \sqrt {\tfrac {1}{4} + \upsilon ^2}. \end{equation}
Here,
$J_\nu$
and
$Y_\nu$
denote Bessel functions of the first and second kind, respectively, with
$\nu$
designating their order. In view of
$\upsilon \geq 0$
, one realises that
$\nu \geqslant 1/2$
. Then, given the direct bearing of
$\nu$
on the solution character, it is adopted hereafter as the primary separation parameter instead of
$\upsilon$
.
Next, the boundary condition on
$f_0$
may be evaluated at the origin using the asymptotic behaviour of
$f_0/R$
at
$R \rightarrow 0^+$
. One has
\begin{align} \frac {f_0(R)}{R} \sim C_1 \frac {\bigl ( C_B^{\scriptscriptstyle {(0)}}/2 \bigr )^\nu }{\nu {\Gamma }(\nu )} R^{\nu -1/2} - C_2 \frac {\Gamma (\nu )}{\pi \bigl ( C_B^{\scriptscriptstyle {(0)}}/2 \bigr )^\nu } \frac {1}{R^{\nu +1/2}}, \end{align}
where we take
$C_B^{(0)} \gt 0$
without any loss of generality. To avoid secular behaviour at the origin, it may be seen that two conditions must be satisfied. These consist of
$C_2 = 0$
and
$\nu \gt 1/2$
. The first equality eliminates the divergent contribution from
$Y_\nu$
, whereas the second inequality prevents divergence at
$R \to 0^+$
. As for the rigid-wall condition at
$R = 1$
, it translates into a vital solvability condition, specifically
\begin{equation} J_\nu \bigl ( C_B^{\scriptscriptstyle {(0)}} \bigr ) \equiv J_\nu (j_{\nu ,m}) = 0, \end{equation}
where
$j_{\nu ,m}$
stands for the
$m$
th non-zero root of
$\mathinner {J_\nu (z)}$
.
Moving on to the ODE that controls polar variations, which is given by (3.17), an exact solution for
$g_0(\phi )$
may be found in terms of the Gaussian hypergeometric function
${}_2F_1$
(Polyanin & Zaitsev Reference Polyanin and Zaitsev2018):
\begin{align} g_0(\phi ) = C_3 \mathinner {{}_2F_1\!\left ( \substack {\tfrac {1}{2}\nu - \tfrac {1}{4}, -\tfrac {1}{2}\nu - \tfrac {1}{4}\\[5pt]\tfrac {1}{2}}; \cos ^2{\!\phi }\right )} + C_4 \cos (\phi ) \mathinner {{}_2F_1\!\left ( \substack {\tfrac {1}{2}\nu + \tfrac {1}{4}, -\tfrac {1}{2}\nu + \tfrac {1}{4}\\[5pt]\tfrac {3}{2}}; \cos ^2{\!\phi }\right )}\!. \end{align}
Furthermore, the centreline boundary condition requires that
$g_0'(\phi )/\sin \phi$
remain bounded as
$\phi \to 0^+$
. For the general solution given by (3.22), this requirement is met when the function itself vanishes along the axis of rotation. In fact, by enforcing
$\lim _{\phi \to 0^+} g_0(\phi )=0$
, one gets
\begin{align} \lim _{\phi \rightarrow 0^+} \mathinner {g_0{\left ( \phi \right )}} = C_3\dfrac {\pi ^{1/2}}{{\Gamma }\!\left ( \tfrac {3}{4} - \tfrac {1}{2}\nu \right ) \Gamma \!\left ( \tfrac {3}{4} + \tfrac {1}{2}\nu \right )} + C_4\dfrac {\tfrac {1}{2} \pi ^{1/2}}{ {\Gamma }\!\left ( \tfrac {5}{4} - \tfrac {1}{2}\nu \right ) {\Gamma }\!\left ( \tfrac {5}{4} + \tfrac {1}{2}\nu \right )} = 0. \end{align}
Here too, the resulting equality will be secured when two conditions are met. Firstly, one of the two integration constants,
$C_3$
or
$C_4$
, must be suppressed. Secondly, the permissible values of
$\nu$
must be constrained in such a manner as to eliminate the remaining term. In this work, we find it convenient to set
$C_4=0$
and then proceed to evaluate
$\nu$
while leveraging the infinite product representation of the Gamma function. Specifically, one can substitute the identity
\begin{equation} \frac {1}{{\Gamma }(z)} = z \,\textrm{e}^{\gamma z} \prod _{k=1}^{\infty }\Bigl ( 1 + \frac {z}{k} \Bigr )\,\textrm{e}^{-z/k}, \end{equation}
where
$\gamma \approx 0.577216$
stands for Euler’s constant. By inserting (3.24) into (3.23), suitable values of
$\nu$
may be identified that will observe the necessary constraint. One finds
\begin{equation} \nu = 2k-\tfrac {1}{2}; \quad k = 1,2,3,\ldots \qquad \text{or}\qquad \nu = \!\left \{ \tfrac {3}{2},\tfrac {7}{2},\tfrac {11}{2},\ldots \right \}. \end{equation}
To avoid solutions exhibiting multiple mantles or flow passes, we proceed with the simple partial solution corresponding to the first root at
$k = 1$
and a fixed value of
$\nu = 3/2$
. Given this choice,
$f_0$
and
$g_0$
collapse into
\begin{equation} f_0(R) = C_1 R\, \,\textrm{j}_{1}(\lambda _m R) \qquad \text{and}\qquad g_0(\phi ) = C_3\sin ^2{\phi }, \end{equation}
where
$\mathinner {\textrm{j}_{n}(z)} \equiv \sqrt {\pi /(2z)}\mathinner {J_{n+1/2}(z)}$
refers to the spherical Bessel function of the first kind and
$\lambda _m \equiv j_{3/2,m}$
designates the
$m$
th positive root of the first-order spherical Bessel function of the first kind,
$\mathinner {\textrm{j}_{1}(z)}$
. For the general expressions of
$\textrm{j}_{n}$
and
$\textrm{y}_{n}$
in terms of trigonometric functions, see Appendix B. Note that in (3.26),
$\sqrt {\pi /2}$
is absorbed, with no loss of generality, into the yet to be determined integration constant
$C_1$
. Then, for
$k = 1$
, (3.21) leads to the classical transcendental equation
\begin{equation} \tan {\lambda _m} = \lambda _m \qquad \text{(no cross-flow along rigid wall).} \end{equation}
After some algebra, the asymptotic solution for positive roots of (3.27) may be expressed as
\begin{align} \begin{aligned} & \lambda _m = m \pi + \textstyle \tfrac {1}{2}\pi - \!\left(m\pi + \tfrac {1}{2}\pi \right)^{-1} - \tfrac {2}{3} \!\left(m \pi + \tfrac {1}{2}\pi \right)^{-3} + \cdots ;\quad m = 1, 2, 3, \ldots \\ & \quad \;\approx \mathopen \{ 4.49341, 7.72525, 10.9041, \ldots \mathclose \}. \end{aligned} \end{align}
Hereafter, only the smallest positive root,
$\lambda _1 \approx 4.49341$
, is considered. The significance of these roots lies in their correspondence to the number of mantles evolving in a hemispherical enclosure. For instance,
$m = 1$
will induce a profile exhibiting a single mantle, whereas
$m=2, 3, 4$
will give rise to cyclonic solutions with two, three or four mantles, respectively. In keeping the present article to a manageable length, the compressible analysis of solutions exhibiting multiple mantles is omitted.
For a single mantle, the product of the remaining constants,
$C_1 C_3$
, can be determined from the mass-flux boundary condition given by (3.18). At the outset, the leading-order streamfunction can be fully determined and written as
\begin{equation} \psi _0(R,\phi ) = \tilde {\kappa }R\mathinner {\textrm{j}_{1}(\lambda _1 R)}\sin ^2{\phi }; \qquad \kappa = \tilde {\kappa }\beta _0\mathinner {\textrm{j}_{1}(\lambda _1\beta _0)}. \end{equation}
Next, the leading-order velocity components may be extracted using straightforward differentiation of (3.29). One collects
\begin{align} u_{\!R}^{\scriptscriptstyle {(0)}}(R,\phi ) &= -2 \tilde {\kappa } \textstyle R^{-1}\,\,\textrm{j}_{1}(\lambda _1 R)\cos \phi , \end{align}
\begin{align} u_{\phi }^{\scriptscriptstyle {(0)}}(R,\phi ) &= \tilde {\kappa }\,\bigl [ \textstyle \lambda _1\mathinner {\textrm{j}_{0}(\lambda _1 R)} - R^{-1}\mathinner {\textrm{j}_{1}(\lambda _1 R)} \bigr ] \sin {\phi }, \end{align}
\begin{align} u_\theta ^{\scriptscriptstyle {(0)}}(R,\phi ) &= \tilde {\kappa }\,\,\lambda _1\mathinner {\textrm{j}_{1}(\lambda _1 R)}\sin \phi . \end{align}
3.2.2. Mantle analysis
With
$\psi _0$
in hand, attention can be turned to the recovery of
$\beta _0$
, the radial transition point separating the updraft and downdraft motions in the equatorial plane. The latter can be specified by setting
$u_\phi ^{\scriptscriptstyle {(0)}}(\beta _0,\pi /2)=0$
. Letting
$R=\beta _0$
in (3.30b
), we get
\begin{align} \lambda _1 \beta _0 &= \frac {\mathinner {\textrm{j}_{1}(\lambda _1 \beta _0)}}{\mathinner {\textrm{j}_{0}(\lambda _1 \beta _0)}} = \frac { ( \lambda _1 \beta _0 )^{-2} \sin (\lambda _1\beta _0) - ( \lambda _1 \beta _0 )^{-1} \cos (\lambda _1 \beta _0) }{ (\lambda _1\beta _0)^{-1} \sin (\lambda _1 \beta _0) } \nonumber \\ &=\frac { \sin (\lambda _1\beta _0) - { \lambda _1 \beta _0 } \cos (\lambda _1 \beta _0) }{ \lambda _1\beta _0\sin (\lambda _1 \beta _0) }, \end{align}
where the expanded fraction represents an identical form to the
$\textrm{j}_{1}/\textrm{j}_{0}$
ratio. In fact, by letting
$x\equiv \lambda _1 \beta _0$
, (3.31) can be collapsed into
\begin{equation} \textstyle \tan {x} = x/\bigl ( 1 - x^2 \bigr ), \quad x \in (0,\lambda _1]. \end{equation}
Similarly to the treatment of
$\lambda _1$
,
$\beta _0$
may be reconstructed from a series expansion. Pursuant to the expansion technique used to derive (3.28), we obtain a power series for
$\beta _0$
in terms of
$\pi$
. We find
\begin{equation} \beta _0 = \tfrac {2}{3} \!\left \{ 1 - \tfrac {5}{9}\pi ^{-2} - \tfrac {20}{81} \pi ^{-4} - \!\left ( \tfrac {80}{729} + \tfrac {1119}{729}\pi ^2 \right )\pi ^{-6} + \cdots \right \} \approx 0.610607. \end{equation}
This leading-order approximation for
$\beta$
confirms that, in the absence of compressibility, the ideal outlet radius reduces to approximately
$61\,\%$
of the hemispherical radius. It is reassuring that this analytically calculated
$\beta _0$
coincides with the numerically extracted root for the incompressible Beltramian mantle reported in the same geometric configuration (Williams & Majdalani Reference Williams and Majdalani2021a
).
3.2.3. First-order thermodynamic corrections
Being a parent function,
$\psi _0$
can be relied upon to derive the first-order corrections to the thermodynamic variables,
$\rho _1$
,
$T_1$
and
$p_1$
. Through appropriate manipulations of (3.8a
)–(3.10b
), we arrive at
\begin{align} \mathinner {\rho _1{\left ( R,\phi \right )}} &= -\frac {1}{2} \biggl [ \frac {\mathinner {\psi _0(R,\phi )}}{R^2 \sin {\phi }} \biggr ]^2\,\biggl \{ \biggl [ \lambda _1 R \,\frac { \mathinner {\textrm{j}_{0}(\lambda _1 R)}}{\mathinner {\textrm{j}_{1}(\lambda _1 R)}} - 1 \biggr ]^2 + 4 \cot ^2{\phi } + \lambda _1^2 \biggr \}, \end{align}
\begin{align} \mathinner {p_1{\left ( R,\phi \right )}} & = \gamma \rho _1(R,\phi ), \end{align}
\begin{align} \mathinner {T_1{\left ( R,\phi \right )}} & = (\gamma -1)\rho _1(R,\phi ). \end{align}
Since
$\psi _0 = {O}(\kappa )$
, the resulting expression for
$\rho _1$
in (3.34a
) reveals that the first-order thermodynamic corrections scale with
${O}(\kappa ^2 M_0^2)$
. The relevance of this composite scaling parameter is discussed in § 4.1.
3.2.4. First-order solution
With both
$\psi _0$
and
$\rho _1$
in hand, we may proceed to evaluate the first-order streamfunction,
$\psi _1$
. Prior to its recovery, however, we find it useful to implement a change of independent variables. Guided by the general solution of
$\psi _0$
in (3.26), a dual coordinate transformation may be introduced, specifically
\begin{equation} \xi = \lambda _1 R \qquad \text{and} \qquad \eta = \cos {\phi }, \end{equation}
with the corresponding differentials
\begin{equation} {\rm d}{R} = {\rm d}{\xi }/\lambda _1 \qquad \text{and} \qquad {\rm d}{\phi } = -{\rm d}{\eta }/\textstyle (1 - \eta ^2)^{1/2}. \end{equation}
Here,
$0\lt \xi \leqslant \lambda _1$
and
$0 \lt \eta \leqslant 1$
. When expressed in terms of these coordinates, the leading-order streamfunction, Stokes operator and gradient become
\begin{align} & \varPsi _0(\xi ,\eta ) = \psi _0(\xi /\lambda _1,\arccos {\eta }) = \textstyle \lambda _1^{-1} \tilde {\kappa } \,(1 - \eta ^2) \,\xi \,\,\textrm{j}_{1}(\xi ), \end{align}
\begin{align} & \textrm{L}^2\!{f} = \lambda _1^2\,\biggl ( \frac {\partial ^2 {f}}{\partial {\xi }^2 }+ \frac {1-\eta ^2}{\xi ^2}\frac {\partial ^2{f}} {\partial {\eta }^2}\biggr ), \end{align}
\begin{align} & {\boldsymbol{\nabla }} f = \lambda _1\,\biggl [ \hat {\boldsymbol{e}}_{\!R} \frac {\partial f}{\partial \xi } - \hat {\boldsymbol{e}}_\phi \frac {(1 - \eta ^2){}^{1/2}}{\xi }\frac {\partial f}{\partial \eta } \biggr ], \end{align}
where
$\hat {\boldsymbol{e}}_{R}$
and
$\hat {\boldsymbol{e}}_{\phi }$
stand for the spherical unit vectors in the meridional plane. Strictly for convenience, we also relabel the streamfunction as
$\varPsi (\xi ,\eta )$
in lieu of
$\psi (R,\phi )$
in this transformed-variable domain. Defining
$\mathinner {\varPsi _1{ ( \xi ,\eta )}} \equiv \mathinner {\psi _1{ ( \xi /\lambda _1,\arccos {\eta } )}}$
, the first-order streamfunction equation turns into
\begin{align} \frac {\partial ^2 \varPsi _1}{\partial \xi ^2} + \frac {1 - \eta ^2}{\xi ^2} \frac {\partial ^2 \varPsi _1}{\partial {\eta }^2} + \varPsi _1 & = -2\frac {{C_B^{\scriptscriptstyle {(1)}}}}{\lambda _1}\varPsi _0 + \biggl [ \varPsi _0 - \frac {1}{2}\biggl ( \frac {\partial \varPsi _0}{\partial \xi }\frac {\partial }{\partial \xi } + \frac {1 - \eta ^2}{\xi ^2} \frac {\partial \varPsi _0}{\partial \eta } \frac {\partial }{\partial \eta } \biggr ) \biggr ] \nonumber \\ & \quad \times \biggl \{ \frac {1}{\xi ^2 (1 - \eta ^2)}\biggl [ {\varPsi ^2_0} + \biggl ( \frac {\partial \varPsi _0}{\partial \xi } \biggr )^2 + \frac {1 - \eta ^2}{\xi ^2} \biggl ( \frac {\partial \varPsi _0}{\partial \eta } \biggr )^2 \biggr ] \biggr \}, \end{align}
with
\begin{align} & \mathinner {\frac {\partial \varPsi _1}{\partial \eta }(\lambda _1,\eta )} = 0 & \text{(no cross-flow along rigid wall),} \end{align}
\begin{align} & \lim _{\eta \rightarrow 1^{-}} \frac {1}{(1 - \eta ^2)^{1/2}}\frac {\partial \varPsi _1}{\partial \xi } = 0 & \text{(no cross-flow along centreline),} \end{align}
\begin{align} & \varPsi _1(\lambda _1\beta _0,0) = 0 & \text{(mass balance along exit plane),} \end{align}
\begin{align} & \lim _{\substack {\xi \to 0^{+} \\ \eta \to 1^{-}} } \frac {\varPsi _1}{\xi \, (1 - \eta ^2)^{1/2}} = 0 & \text{(bounded flow profile at origin).} \end{align}
Following the leading-order analysis, the fourth boundary condition given by (3.39d ) may be applied to suppress any attendant singularities. The primary challenge in solving (3.38) lies in determining the appropriate separation of variables that will decompose the PDE into a linear set of coupled ODEs. Naturally, a general product solution to a linear PDE may be achieved through superposition, i.e. by juxtaposing multiple product solutions of the form
\begin{equation} \varPsi _1(\xi ,\eta ) = F_1(\xi )\,G_1(\eta ) + F_2(\xi ) \, G_2(\eta ) + \cdots + F_n(\xi ) \, G_n(\eta ). \end{equation}
In (3.40), the individual functions
$\{\mathinner {F_i{ ( \xi )}},\mathinner {G_i{ ( \eta )}}; i=1,2,\ldots ,n \}$
, which are indexed by
$i$
, only need to collectively satisfy the right-hand side of (3.38). Each pair may only satisfy one or more source terms.
Upon close examination of the boundary conditions, particularly (3.39b
), we find it essential for
$\varPsi _1$
to decay as rapidly as
$(1-\eta ^2)^{1/2}$
when
$\eta \to 1^{-}$
. This prompts us to set
$G_i(\eta )=(1-\eta ^2)^{i}$
, where the suitability of different exponents may be systematically investigated through backward substitution into (3.39b
), (3.39d
) and (3.38). To further streamline the analysis, we divide
$\varPsi _1$
by
$\lambda _1\tilde {\kappa }^3$
and eliminate the constant in
$\varPsi _0^3$
that emerges on the right-hand side of (3.38). After some algebra, we arrive at a viable
$\varPsi _1$
expansion, particularly
\begin{equation} \textstyle \varPsi _1(\xi ,\eta ) = \lambda _1 \tilde {\kappa }^{3} \!\left [(1 - \eta ^2) \, F_1(\xi ) + (1 - \eta ^2)^2\, F_2(\xi ) \right ]. \end{equation}
In the above, terms comprising
$(1 - \eta ^2)^n$
for
$n\gt 2$
are not included because their corresponding coefficients end up vanishing upon substitution into (3.38). This proposed ansatz may be shown to satisfy (3.39b
) identically, thereby reducing the problem to that of securing the remaining conditions in (3.39a
), (3.39c
) and (3.39d
). Verification of this ansatz may be readily carried out by evaluating the right-hand side of (3.38) while grouping terms in powers of
$(1-\eta ^2)$
. This operation yields a pair of ODEs for the functions
$F_1$
and
$F_2$
that may be gathered into
\begin{align} \xi ^2\,\frac {{\rm d}^2 {F_2}}{{\rm d} \xi ^2} + \bigl ( \xi ^2 - 12 \bigr )F_2 &= \xi \bigl \{ \xi ^2 \mathinner {\textrm{j}_{1}^3(\xi )} - \xi \mathinner {\textrm{j}_{1}^2(\xi )} \mathinner {\textrm{j}_{2}(\xi )} - 10 \mathinner {\textrm{j}_{1}(\xi )}\mathinner {\textrm{j}_{2}^2(\xi )} + 4\xi \mathinner {\textrm{j}_{2}^3(\xi )} \nonumber\\ &\quad - \bigl [ \xi \mathinner {\textrm{j}_{2}(\xi )} - 2\mathinner {\textrm{j}_{1}(\xi )} \bigr ]^2 \,\textrm{j}_{3}(\xi ) \bigr \}\equiv \mathinner {g_2(\xi )}, \end{align}
\begin{align} \xi ^2 \frac {{\rm d}^2{F_1}}{{\rm d} \xi ^2 } + \bigl ( \xi ^2-2 \bigr )\,F_1 &= -8F_2 - 2\bigl [ \tilde {C}_B^{(1)} \xi ^3 - \xi \mathinner {\textrm{j}_{1}^2(\xi )} - 8\mathinner {\textrm{j}_{1}(\xi )} \mathinner {\textrm{j}_{2}(\xi )} + 3\xi \mathinner {\textrm{j}_{2}^2(\xi )} \bigr ]\mathinner {\textrm{j}_{1}(\xi )} \nonumber\\ & \equiv \mathinner {g_1(\xi )}, \\[9pt] \nonumber \end{align}
where
$\tilde {C}_B^{\scriptscriptstyle {(1)}} \equiv {C_B^{\scriptscriptstyle {(1)}}}\lambda ^{-3}_1 \tilde {\kappa }^{-2}$
. In terms of
$F_i$
, the boundary conditions in (3.39a
), (3.39b
) and (3.39c
) turn into
\begin{align} & \mathinner {F_i{\left ( \lambda _1 \right )}} = 0 & \text{(no cross-flow along rigid wall),} \end{align}
\begin{align} & \lim _{\eta \to 1^{-}}(\textstyle 1-\eta ^2)^{i-1/2} F_i'(\xi ) = 0 & \text{(no cross-flow along centreline),} \end{align}
\begin{align} & F_1(\lambda _1\beta _0) + F_2(\lambda _1 \beta _0)= 0 & \text{(mass balance along exit plane),} \end{align}
\begin{align} & \lim _{\xi \to 0^{+}} \textstyle \xi ^{-1}F_i(\xi ) = 0 & \text{(bounded flow profile at origin),} \end{align}
where, let us recall,
$i=\{1,2\}$
. As expected, the solution corresponding to the homogeneous forms of (3.42a
) and (3.42b
) gives rise to spherical Bessel functions of the first and second kind,
$\textrm{j}_n$
and
$\textrm{y}_n$
, respectively. A detailed description of these is given by Olver et al. (Reference Olver, Lozier, Boisvert and Clark2010). For both
$i=\{1,2\}$
, the homogeneous solution,
$F_{i,{h}}$
, takes the form
\begin{equation} \mathinner {F_{i,\text{h}}{\left ( \xi \right )}} = C_{2i-1} \xi \mathinner {\textrm{j}_{2i-1}(\xi )} + C_{2i}\xi \mathinner {\textrm{y}_{2i - 1}(\xi )}, \end{equation}
where
$C_{2i}$
and
$C_{2i-1}$
refer to the integration constants.
At this juncture, as we seek to resolve the tightly coupled set in (3.42a ) and (3.42b ) , it may be helpful to recall some basic theory. Considering a generally inhomogeneous, second-order, linear ODE of the form
\begin{equation} f_2(x)\,y'' + f_1(x)\,y' + f_0(x) \, y = g(x) \quad \text{with} \quad f_2(x) \neq 0, \end{equation}
with two homogeneous solutions,
$y_1$
and
$y_2$
, a direct solution may be achieved using
\begin{align} y(x) = C_1 y_1(x) + C_2 y_2(x) + y_2(x) \int \! \frac {y_1(x) g(x)}{f_2(x) \mathinner {\mathscr{W}(y_1,y_2)}} {\rm d}{x} - \mathinner {y_1{\!\left ( x \right )}}\int \! \frac {y_2(x) g(x)}{f_2(x) \mathinner {\mathscr{W}(y_1,y_2)}}{\rm d}{x}, \end{align}
where
$\mathinner {\mathscr{W}(y_1,y_2)} \equiv y_1y_2' - y_2y_1'$
defines the Wronskian of
$y_1$
and
$y_2$
. Presently, one can employ the corresponding
$F_{i,{h}}$
in (3.44) to determine the Wronskian using
\begin{equation} \mathscr{W}\,\mathopen [\xi \mathinner {\textrm{j}_{2i-1}{(\xi)}},\xi \mathinner {\textrm{y}_{2i-1}{(\xi)}}\mathclose ] = 1. \end{equation}
Then, using
$f_2(\xi ) = \xi ^2$
, the general solutions to (3.42a
) and (3.42b
) may be recovered. For each index
$i$
, one gets
\begin{align} \mathinner {F_i(\xi )} &= \xi \mathinner {\textrm{j}_{2i-1}{\left ( \xi \right )}} \biggl [ C_{2i-1} - \int \frac {\mathinner {\textrm{y}_{2i-1}{(\xi)}}\mathinner {g_i{(\xi)}}}{\xi }{\rm d}{\xi } \biggr ]\nonumber\\&\quad + \xi \mathinner {\textrm{y}_{2i-1}{\left ( \xi \right )}} \biggl [ C_{2i} + \int \frac {\mathinner {\textrm{j}_{2i-1}{(\xi)}}\mathinner {g_i{(\xi)}}}{\xi } {\rm d}{\xi } \biggr ]. \end{align}
In (3.48),
$g_i$
corresponds to the right-hand sides of (3.42a
) and (3.42b
). Using spherical Bessel function identities, these integrals can be evaluated to secure
$F_i(\xi )$
analytically. While the resulting expressions tend to be lengthy, they can be expressed in terms of well-known functions to the extent of obviating the need for numerical integration. For the sake of completeness, the corresponding detail is provided in Appendix C.
Finally, to complete the derivation of
$\varPsi _1$
, the integration constants
$C_1$
to
$C_4$
must be determined by satisfying (3.43a
), (3.43b
), (3.43c
) and (3.43d
). Although the explicit forms of
$F_i$
are not presented here, an expansion of
$\mathinner {F_i{ ( \xi )}}/\xi$
about
$\xi = 0$
confirms that
\begin{equation} \frac {\mathinner {F_i{(\xi)}}}{\xi }= C_{2i-1}\frac {\xi ^{2i-1}}{{\!\left ( 4i-1 \right )}!!} - C_{2i}\frac {{\!\left ( 4i-3 \right )}!!}{\xi ^{2i}} + o(\xi ), \quad \text{for} \quad i= \{1,2\}, \end{equation}
where the double factorial symbol
$n!!=n(n-2)(n-4)\ldots$
denotes a factorial that skips every other integer. Next, to satisfy the boundedness condition at the origin, (3.43d
), we find it necessary to set
$C_{2i}=0$
. The remaining constants,
$C_1, C_3,$
and the eigenvalue correction,
$\tilde {C}{}_B^{\scriptscriptstyle {(1)}}$
, may be evaluated by applying the remaining boundary conditions in (3.43a
) and (3.43c
). These enable us to retrieve
\begin{align} C_3 &= \!\left . {\!\left [ \int \frac {\textrm{y}_{3}(\xi ) g_2(\xi )}{\xi } {\rm d}{\xi } - \frac {\mathinner {\textrm{y}_{3}(\xi )}}{\mathinner {\textrm{j}_{3}{(\xi)}}} \int \frac {\textrm{j}_{3}(\xi ) g_2(\xi )}{\xi } {\rm d}{\xi } \right ]}\right |_{\xi = \lambda _1}, \end{align}
\begin{align} C_1 &=\!\left \{\int \frac {\textrm{y}_{1}(\xi ) g_1(\xi )}{\xi } {\rm d}{\xi } - \!\left [\frac {\textrm{y}_{1}\bigl ( \xi \bigr )}{\textrm{j}_{1}\bigl ( \xi \bigr )} \int \frac {\textrm{j}_{1}(\xi ) g_1(\xi )}{\xi }{\rm d}{\xi } + \frac {\textrm{j}_{3}\bigl ( \xi \bigr )}{\textrm{j}_{1}(\xi )}{\!\left ( C_3 - \int \frac {\textrm{y}_{3}(\xi ) g_2(\xi )}{\xi }{\rm d}{\xi } \right )} \right . \right . \nonumber \\ &\quad \!\left . \!\left .\!\left . + \frac {\textrm{y}_{3}\bigl ( \xi \bigr )}{\textrm{j}_{1}\bigl ( \xi \bigr )} \int \frac {\textrm{j}_{3}(\xi ) g_2(\xi )}{\xi } {\rm d}{\xi }\right ]\right \}\right |_{\xi = \lambda _1\,\beta _0}, \end{align}
\begin{align} \tilde {C}_B^{(1)} &= \bigg \{ {\!\left [ \int \xi ^2 \mathinner {\textrm{j}_{1}^2(\xi )} {\rm d}{\xi } \right ]}^{-1} \bigg [\int \biggl ( \mathinner {\textrm{j}_{1}^2(\xi )} + \frac {8\mathinner {\textrm{j}_{1}(\xi )}\textrm{j}_{2}(\xi )}{\xi } - 3\mathinner {\textrm{j}_{2}^2(\xi )} \biggr )\mathinner {\textrm{j}_{1}^2(\xi )}{\rm d}{\xi } \notag \\ & \quad - 4\int \frac {\textrm{j}_{1}(\xi )F_2(\xi )}{\xi }{\rm d}{\xi } \bigg ] \bigg \} \bigg |_{\xi = \lambda _1}. \end{align}
With
$\psi _1$
being fully determined, the first-order velocity components may be deduced from the expansion of (A2) and (2.20). Reverting back to the original variables, we collect
\begin{align} u_{\!R}^{\scriptscriptstyle {(1)}} & = -\frac {1}{R^2\sin {\phi }}\Bigl ( \frac {\partial \psi _1}{\partial \phi } - \rho _1 \frac {\partial \psi _0}{\partial \phi } \Bigr ), \end{align}
\begin{align} u_\phi ^{\scriptscriptstyle {(1)}} & = \frac {1}{R\sin {\phi }} \Bigl ( \frac {\partial \psi _1}{\partial R} - \rho _1 \frac {\partial \psi _0}{\partial R} \Bigr ), \end{align}
\begin{align} u_\theta ^{\scriptscriptstyle {(1)}} & = \biggl ( \frac {\psi _1}{\psi _0} + \frac {C_B^{\scriptscriptstyle {(1)}}}{\lambda _1} \biggr )\,u_\theta ^{\scriptscriptstyle {(0)}}. \end{align}
Additionally, the first-order correction to the outlet radius,
$\beta _1$
, may be extracted from (3.7) to get
\begin{equation} \beta _1\approx 8.61201\kappa ^2. \end{equation}
Lastly, by combining (3.52) with (3.33), the two-term expansion of
$\beta$
may be expressed as
\begin{equation} \beta = 0.610607 + 8.61201\kappa ^2 M_0^2 + {O}\bigl ( M_0^4 \bigr ). \end{equation}
The formal expressions for the second-order thermodynamic corrections,
$\rho _2$
,
$p_2$
and
$T_2$
, are given by (3.8b
) and (3.11); these quantities are evaluated and discussed in § 4.4.
4. Results and discussion
4.1. A composite group parameter
Prior to detailing the specific flow characteristics, it may be instructive to examine the role of the inflow parameter,
$\kappa =\rho _{{i}} A_{{i}}/(2\pi \rho _0 a^2)$
. From the boundary condition in (2.16c
), it may be seen that the streamfunction remains linearly proportional to the mass influx, with
$\psi = {O}(\kappa )$
and, in turn,
$\boldsymbol{u}_0 = {O}(\kappa )$
.
This particular scaling has a direct bearing on the compressibility corrections. The first-order density correction, for instance, is given by (3.34a
) as
$\rho _1 \propto \mathopen \| \boldsymbol{u}_0 \mathclose \|^2$
. Since
$\boldsymbol{u}_0 = {O}(\kappa )$
, it follows that
$\rho _1 = {O}(\kappa ^2)$
. The first-order density correction,
$M_0^2 \rho _1$
, becomes of order
$M_0^2 \kappa ^2$
. This behaviour proves to be characteristic of all first-order quantities, which simply depend on the
$(\kappa M_0)^2$
product. We are thus led to define a single composite parameter,
$M_\kappa \equiv \kappa M_0$
, that naturally combines the effects of the reference Mach number with those of the inflow geometry and mass flux into the hemispherical domain.
In the forthcoming analysis, all results are presented in a manner that showcases this characteristic group parameter. In so doing, leading-order quantities are normalised by the appropriate power of
$\kappa$
(e.g.
$\mathopen \| \boldsymbol{u} \mathclose \|_0/\kappa$
), while first-order corrections are presented such that their magnitudes remain independent of both
$\kappa$
and
$M_0$
(e.g.
$\mathopen \| \boldsymbol{u} \mathclose \|_1/\kappa ^3$
). This allows the ensuing graphical results to capture fundamental solution shapes, whose magnitudes can be determined through straightforward multiplication by the relevant scaling parameters.
4.2. Streamline and mantle displacements
We begin in figure 3 by evaluating the streamline displacements that are induced by increasing dilatational effects.
Vector streamlines and loci of
$u_z=0$
(mantle, in blue) and
$(u_\theta )_{max}$
(vortex core radius, in red) for (a) incompressible and (b) compressible motions in the
$r$
–
$z$
plane. Using a
$\kappa$
-normalised streamfunction, streamlines, mantle and
$(u_\theta )_{max}$
lines in part (b) correspond to
$M_\kappa ^2 = 0$
(
),
$0.003$
(
) and
$0.006$
(
).
4.2.1. Compressible streamlines
For the incompressible motion at
$M_0 = 0$
, the streamlines in figure 3(a) may serve as a visual benchmark for tracing particle trajectories that are unaffected by Mach number variations. In contrast, figure 3(b) illustrates the evolution of streamline curvatures under increasing compressibility levels for three characteristic values of
$M_\kappa ^2 = 0$
,
$0.003$
and
$0.006$
. To enable direct comparisons, contours of the normalised streamfunction,
$\psi /\kappa = 0.00826,\ 0.0744,\ 0.207,\ 0.669$
and
$1$
, are plotted for each
$M_\kappa$
. Therein, one may ascertain that compressibility effects cause streamlines of constant
$\psi /\kappa$
to deviate from their incompressible counterparts. Graphically, it may be seen that as
$M_\kappa$
is raised, higher tangential forces push the streamlines radially outward from the centreline, whereas an increase in axial forces causes them to overshoot relative to their incompressible trajectories as the chamber apex is approached. As a result, the compressible flow turning becomes sharper, with higher-velocity streamlines bending more abruptly and exiting the enclosure across the equatorial plane closer to the chamber wall than their slower zeroth-order counterparts. Interestingly, the locus of points prescribed by
$\psi _1(R,\phi ) = 0$
delineates a surface where the compressible and incompressible forms of the streamfunction become identical. In fact, the multiple
$\psi$
-intersections that appear in figure 3(b) occur among streamlines taken at different Mach numbers. Downstream of this surface, the separation between incompressible and compressible trajectories becomes gradually more noticeable due to the outward displacement of the higher-momentum outflow, particularly as it exits the domain at a larger radius than in the case of a slower solenoidal motion.
4.2.2. Axial mantle
Another manifestation of dilatational effects stands in the radially outward expansion of the axisymmetric mantle surface. As first introduced in § 2.2.2, the latter refers to the interfacial layer separating the inner downdraft from the outer updraft regions and where, quite specifically, the axial velocity vanishes. This behaviour may be readily inferred from the streamlines in figure 3(b). The induced outward expansion can be quantified by relating the normalised outlet mantle radius,
$\beta$
, to the compressibility parameter,
$M_\kappa$
, as captured in the equatorial plane through (3.53). More globally, the axial mantle can be pinpointed by setting
$u_z=0$
, namely,
\begin{equation} {\partial } \psi /{\partial } r= 0. \end{equation}
Based on (4.1), analytical expressions can be derived for the leading- and first-order contributions to the mantle radius,
$r_{{m}}$
, as functions of the axial elevation,
$z$
. To this end, the mantle radius,
$r_{{m}}(z)$
, can be expanded using
\begin{equation} r_{{m}}(z) = r_{{m}}^{\scriptscriptstyle {(0)}}(z) + M_0^2 r_{{m}}^{\scriptscriptstyle {(1)}}(z) + {O}\bigl ( M_0^4 \bigr ), \end{equation}
where
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
denotes the incompressible mantle radius and
$r_{{m}}^{\scriptscriptstyle {(1)}}(z)$
captures its first-order correction. To make further headway, one may insert (4.2) and the expansion of
$\psi$
into (4.1), collect terms in like powers of
$M_0$
and return with two equations for the mantle surface at successive orders of
$M_0^2$
. As shown in Appendix D, one simply recovers
\begin{equation} {\partial }\psi _0/{\partial } r = 0 \quad \text{and} \quad {\partial } \psi _1/{\partial } r + r_{{m}}^{\scriptscriptstyle {(1)}} {\partial }^2\psi _0/{\partial } r^2 = 0, \end{equation}
where all derivatives of
$\psi _0$
and
$\psi _1$
can be evaluated at the leading-order mantle radius,
$r=r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
.
Next, by substituting
$\psi _0$
and
$\psi _1$
from (3.29) and (3.41) into (4.3), a transcendental relation for
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
can be identified along with a linear algebraic equation for
$r_{{m}}^{\scriptscriptstyle {(1)}}(z)$
. An additional transformation from spherical to cylindrical coordinates leads to
\begin{equation} \bigl [ 2 \bigl ( z / r_{{m}}^{\scriptscriptstyle {(0)}} \bigr ){}^2 - 1 \bigr ] \bigl [ \sin \bigl ( \lambda _1 R_{{m}}^{\scriptscriptstyle {(0)}} \bigr ) / \bigl ( \lambda _1 R_{{m}}^{\scriptscriptstyle {(0)}} \bigr ) - \cos \bigl ( \lambda _1 R_{{m}}^{\scriptscriptstyle {(0)}} \bigr ) \bigr ] + \sin \bigl ( \lambda _1 R_{{m}}^{\scriptscriptstyle {(0)}} \bigr ) = 0, \end{equation}
where
$R_{{m}}^{\scriptscriptstyle {(0)}} = \sqrt {(r_{{m}}^{\scriptscriptstyle {(0)}}){}^2 + z^2}$
. Although an exact solution to (4.4) does not seem to be tractable from a finite combination of elementary functions, a relatively straightforward and closed-form approximation for
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
may be obtained:
\begin{equation} \bigl [ r_{{m}}^{\scriptscriptstyle {(0)}}(z) \bigr ]^2 \approx \beta _0^2 \bigl ( 1-z^2 \bigr )\tfrac {1 + \tfrac {1}{5} \bigl ( 8\beta _0{}^{-2} - 21 \bigr )z^2 }{1 - \tfrac {1}{5}\bigl ( 32\beta _0{}^2 - 11 \bigr )z^2}. \end{equation}
In practice, this fractional approximation identically satisfies (4.4) at the vertical endpoints,
$z = \{0,1\}$
, while entailing a relative error of under
$0.17\,\%$
over the entire mantle surface. Indeed, a cursory inspection of figure 4(a) confirms that the approximation in (4.5) remains visually indiscernible from the numerically obtained solution at the given resolution. For the reader’s convenience, this compact expression for the mantle surface is depicted in figure 3(a). As for the Mach-number-dependent mantle lines in figure 3(b), they are deduced directly from the second member of (4.3). In constructing (4.5), it may be helpful to recognise that the geometric configuration of the mantle surface may be interpreted as the upper segment of a prolate spheroid that can be nearly prescribed by
$r^2 = \beta _0^2(1-z^2)$
, i.e. a function that satisfies the mantle’s endpoint conditions. The surface curvature of this shape can be subsequently perturbed in accordance with the rational expression
$[1+{1}/{5}(8\beta _0^{-2} - 21)z^2]/[1 - {1}/{5}(32\beta _0^2 - 11)z^2]$
. Guided by Hermite interpolation concepts, the perturbing function can be formulated such that the resulting approximation for
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
matches the numerical solution and several of its derivatives at the domain boundaries, which are bracketed by
$z = \{0,1\}$
. More specifically, we find it optimal for the analytical approximation to match the first derivative of the numerical solution at
$z = 0$
as well as the derivatives up to the fourth order at
$z = 1$
. Then, once
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
is determined,
$r_{{m}}^{\scriptscriptstyle {(1)}}(z)$
can be readily retrieved from the second relation in (4.3).
Comparison between analytical solutions (
) and the numerical solution (
$\circ$
) for the mantle radius,
$r_{{m}}(z)$
. The numerical roots are obtained from (4.3) using a Newton–Raphson solver. These plots display (a) the leading-order distribution,
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
, and (b) the first-order compressible correction,
$r_{{m}}^{\scriptscriptstyle {(1)}}(z)$
.
The accuracy of the resulting variations in
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
and
$r_{{m}}^{\scriptscriptstyle {(1)}}(z)$
may be ascertained from figures 4(a) and 4(b). Therein, both leading- and first-order mantle variations can be seen to compare rather favourably with their numerical solutions. The latter are extracted from (4.4) using a Newton–Raphson root-finding solver. One may also infer that as the axial elevation is reduced from the apex to the outlet plane,
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
continues to increase monotonically. Particularly at
$z = 0$
, the incompressible mantle radius reaches a maximum value of approximately
$0.6106$
. In the same vein, the compressibility-induced contributions vanish at the apex, where a stagnation region seems to develop. Here too, as the equatorial plane is approached, a monotonic increase in
$r_{{m}}^{\scriptscriptstyle {(1)}}(z)$
is observed, namely, where dilatational forces become so appreciable that they produce a maximum correction of
$r_{{m}}^{\scriptscriptstyle {(1)}}(0) = (r_{{m}}^{\scriptscriptstyle {(1)}})_{max} \approx 8.612\kappa ^2$
.
4.3. Velocity field
In the interest of clarity, the analysis is presented using
$\boldsymbol{u}$
and its cylindrical components (
$u_r$
,
$u_{\theta }$
,
$u_z$
). These components relate to their spherical counterparts
$(u_{\!R}, u_\phi , u_{\theta })$
through the coordinate transformations
\begin{equation} u_r = u_{\!R} \sin {\phi } + u_\phi \cos {\phi } \quad \text{and}\quad u_z = u_{\!R}\cos {\phi } - u_\phi \sin {\phi }. \end{equation}
Before proceeding further, it may be helpful to remark that, starting with figures 5(a) and 5(c), all forthcoming contour plots display a standardised colour scheme. Accordingly, warm hues designate positive scalar quantities, whereas cool hues denote negative quantities. As for the colour intensity, it scales directly with the magnitude of a given function. For enhanced contrast, darker shades are used to reflect higher absolute values, while lighter shades mark regions where the function becomes insignificant.
Total velocity magnitudes in (a) leading-order and (c) first-order isocontours. Radial variations are shown using (b) incompressible and (d) compressible velocity components at equispaced axial stations of
$z = 0$
(
),
$0.25$
(
),
$0.50$
(
) and
$0.75$
(
). Everywhere,
$x\equiv r \cos \theta$
denotes the polarised cylindrical radius.
4.3.1. Total velocity
Attention is first directed towards the magnitude of the velocity field,
$\mathopen \| \boldsymbol{u} \mathclose \|$
. Using the decomposition technique described previously,
$\mathopen \| \boldsymbol{u} \mathclose \|$
can be expanded asymptotically for small Mach numbers. Retaining the initial two terms, we have
\begin{equation} \begin{aligned} \mathopen \| \boldsymbol{u} \mathclose \| & = \bigl ( \boldsymbol{u}_0\boldsymbol{\cdot }\boldsymbol{u}_0 + 2 M_{0}^2\,\boldsymbol{u}_1\boldsymbol{\cdot }\boldsymbol{u}_0+\cdots \bigr )^{1/2} = {\left ( \boldsymbol{u}_0\boldsymbol{\cdot }\boldsymbol{u}_0 \right )}^{1/2} + M_{0}^2 \,\boldsymbol{u}_1 \boldsymbol{\cdot }\hat {\boldsymbol{u}}_0 + \mathinner {{O}\bigl (M_{0}^4\bigr )} \\ & = \mathopen \| \boldsymbol{u} \mathclose \|_0 + M_{0}^2\,\mathopen \| \boldsymbol{u} \mathclose \|_1 + \mathinner {{O}\bigl (M_{0}^4\bigr )}; \quad \mathopen \| \boldsymbol{u} \mathclose \|_1 \equiv \boldsymbol{u}_1 \boldsymbol{\cdot }\hat {\boldsymbol{u}}_0, \end{aligned} \end{equation}
where
$\hat {\boldsymbol{u}}_0 \equiv \boldsymbol{u}_0/\mathopen \| \boldsymbol{u}_0 \mathclose \|$
defines the unit vector that remains aligned with the leading-order component,
$\boldsymbol{u}_0$
. In the above, the scalar quantity,
$\mathopen \| \boldsymbol{u} \mathclose \|_1 \equiv \boldsymbol{u}_1 \boldsymbol{\cdot }\hat {\boldsymbol{u}}_0$
, denotes the first-order dilatational correction to the velocity magnitude along the base streamlines. Note that
$\mathopen \| \boldsymbol{u} \mathclose \|_1$
constitutes a scalar that can assume both positive and negative values, unlike a conventional Euclidean norm.
In what follows, the spatial distribution of the leading-order velocity magnitude, normalised by the inflow parameter,
$\kappa$
, and the corresponding first-order correction term,
$\mathopen \| \boldsymbol{u} \mathclose \|_1$
, normalised by
$\kappa ^3$
, are illustrated in figure 5. Specifically, figure 5(a) displays isocontours of
$\mathopen \| \boldsymbol{u} \mathclose \|_0/\kappa$
, whereas figure 5(b) presents radial profiles of this quantity at select
$z$
stations. In a similar manner, figure 5(c) displays isocontours of the normalised correction term,
$\mathopen \| \boldsymbol{u} \mathclose \|_1/\kappa ^3$
, with its radial variations detailed in figure 5(d).
The spatial distribution of the leading-order velocity magnitude,
$\mathopen \| \boldsymbol{u} \mathclose \|_0$
, is illustrated in figures 5(a) and 5(b). As one would expect, it attains its maximum value at the chamber origin and decays monotonically towards the boundaries. In contrast, the first-order correction,
$\mathopen \| \boldsymbol{u} \mathclose \|_1$
, exhibits a more complex structure (figures 5
c and 5
d); its maximum is displaced radially outward from the origin and is located in the outlet plane. This displacement is a direct consequence of dilatational forces redistributing momentum as the flow accelerates towards the exit. The specific values and locations of these extrema are summarised in table 1.
Locations and magnitudes of the kinematic and thermodynamic extrema.
4.3.2. Radial velocity profile
Next, attention is turned to the spatial characteristics of the cylindrical radial velocity,
$u_r$
, whose distribution is illustrated in figure 6. In cyclonic flows, irrespective of the chamber configuration, the radial velocity remains invariably smaller in magnitude compared with its axial and azimuthal counterparts (Sheen et al. Reference Sheen, Chen, Jeng and Huang1996; Cervantes & Gustavsson Reference Cervantes and Gustavsson2006). These relative magnitudes are influenced by both inlet conditions and geometric considerations. In fact, a small radial velocity may be deemed characteristic of most swirl-driven flows, particularly where the majority of fluid momentum remains concentrated in the azimuthal component (Majdalani Reference Majdalani2012, Reference Majdalani2022).
Cylindrical radial velocity profiles in (a) leading-order and (c) first-order isocontours. Radial variations are shown using (b) incompressible and (d) compressible velocity components at
$z = 0$
(
),
$0.25$
(
),
$0.50$
(
) and
$0.75$
(
).
As specified by the centreline boundary conditions, it can be readily confirmed that both the incompressible
$u_r^{\scriptscriptstyle {(0)}}$
and compressible
$u_r^{\scriptscriptstyle {(1)}}$
components vanish along the centreline. Notwithstanding this common property, the spatial distribution and polarity of both
$u_r^{\scriptscriptstyle {(0)}}$
and
$u_r^{\scriptscriptstyle {(1)}}$
differ significantly across the domain. For instance, figure 6(a) provides isocontours of
$u_r^{\scriptscriptstyle {(0)}}/\kappa$
, whereas figure 6(b) displays data at evenly spaced axial elevations of
$z = 0$
,
$0.25$
,
$0.50$
and
$0.75$
. Through visual inspection, these results indicate quite clearly that
$u_r^{\scriptscriptstyle {(0)}}$
remains negative throughout the chamber, as prescribed by the requirement for a radially inward flow (Colonius, Lele & Moin Reference Colonius, Lele and Moin1991). Indeed, this observation aligns well with the classical theory of vortex dynamics, where pressure gradients induced by swirling motion direct the flow radially inward (Wu, Ma & Zhou Reference Wu, Ma and Zhou2006).
Conversely, the compressible contribution
$u_r^{\scriptscriptstyle {(1)}}$
, which is illustrated in figures 6(c) and 6(d), leads to a polarity reversal across a quasi-ellipsoidal boundary corresponding to
$u_r^{\scriptscriptstyle {(1)}} = 0$
. This boundary separates regions of positive and negative contributions in a manner that is reminiscent of that achieved by the mantle surface described in § 4.2. In fact, such behaviour proves to be consistent with the spatial variations of the radial velocity field encountered in the treatment of compressible Beltramian fields in right-cylindrical chambers (Cecil et al. Reference Cecil, Little and Majdalani2023). These outward shifts, which also appear in the streamline plots of figure 3(b), are likely driven by centrifugal forces associated with specific angular momentum redistributions; the latter are prescribed by
$u_{\theta }$
and its relation to
$\mathinner {B{ ( \psi )}}$
.
Another noteworthy distinction between
$u_r^{\scriptscriptstyle {(0)}}$
and
$u_r^{\scriptscriptstyle {(1)}}$
concerns the number and spatial manifestations of extrema within the domain. The incompressible component
$u_r^{\scriptscriptstyle {(0)}}$
exhibits a single global minimum, whereas the compressible correction,
$u_r^{\scriptscriptstyle {(1)}}$
, gives rise to both a minimum and a maximum whose values and locations are evaluated and catalogued in table 1. This behaviour, which is also inferable from the streamline plots of figure 3, can be attributed to the balance between dilatational and centrifugal forces in prescribing the ensuing motion.
4.3.3. Axial velocity profile
The axial velocity component,
$u_z$
, defines the vertically bidirectional nature of this swirl-dominated motion by separating the inner (core) downdraft from the outer (toroidal) updraft. The leading-order component,
$u_z^{\scriptscriptstyle {(0)}}$
, depicted in figures 7(a) and 7(b), establishes this fundamental structure. In this context, the interface given by
$u_z^{\scriptscriptstyle {(0)}}=0$
forms the mantle surface, whose properties are already described in § 4.2. The peak downdraft velocity for the incompressible base flow occurs at the origin, with a value that is presented in table 1.
Axial velocity profiles in (a) leading-order and (c) first-order isocontours. Radial variations are shown using (b) incompressible and (d) compressible velocity components at
$z = 0$
(
),
$0.25$
(
),
$0.50$
(
) and
$0.75$
(
).
Turning now to the first-order correction,
$u_z^{\scriptscriptstyle {(1)}}$
, compressibility appears to have a visible bearing on the axial flow field, as reflected in figures 7(c) and 7(d). This term introduces appreciable spatial variations as it features both negative and positive contributions. More specifically, negative values of
$u_z^{\scriptscriptstyle {(1)}}$
tend to produce a compressibility-induced acceleration of the downdraft and, conversely, a deceleration of the updraft. For example, figure 7(c) contains a centrally located region of negative
$u_z^{\scriptscriptstyle {(1)}}$
in the vicinity of the outlet plane; it also displays a blanketing zone of progressively increasing
$u_z^{\scriptscriptstyle {(1)}}$
in the radially outward direction, particularly near the wall. While extrema for both components can be observed near the outlet, the locations of their most negative values differ. As documented in table 1, the minimum for
$u_z^{\scriptscriptstyle {(1)}}$
is displaced radially outward from the origin. This outward shift may be caused by swirl-induced compression near the outlet; the latter tends to redistribute the axial momentum flux in a manner that displaces the peak downdraft speed slightly away from the centreline.
Perhaps equally revealing is the effect of
$u_z^{\scriptscriptstyle {(1)}}$
on the overall axial profile at different vertical stations, as detailed in figure 7(d). Near the outlet,
$u_z^{\scriptscriptstyle {(1)}}$
remains predominantly negative across the core region, thus leading to an intensification of the downdraft velocity. Radially outward from the core and still near the outlet,
$u_z^{\scriptscriptstyle {(1)}}$
becomes positive and hence additive to the upward speed as the wall is approached. Further into the chamber and towards the apex, this correction switches polarity:
$u_z^{\scriptscriptstyle {(1)}}$
becomes largely positive across most of the radius, suggesting that compressibility primarily acts to quicken the updraft while delaying the downdraft in the upper portions of the chamber. It may therefore be seen that the overall effect of compressibility is not limited to a simple flattening but rather entails a complex spatial redistribution of axial momentum. The visible intensification of the downdraft near the outlet core, where
$u_z^{\scriptscriptstyle {(1)}}$
reaches its largest negative values, may prevent conditions that are conducive to flow reversal and the local formation of recirculation zones (Lilley Reference Lilley1977; Wang & Yang Reference Wang and Yang2018).
4.3.4. Tangential, azimuthal or swirl velocity profile
We now consider the tangential velocity,
$u_{\theta }$
, whose spatial distribution is illustrated in figure 8. The properties of this component are closely tied to the conservation of specific angular momentum,
$B(\psi ) \equiv R \sin {\phi } \, u_{\theta }$
. The radial profiles for the leading-order component,
$u_\theta ^{\scriptscriptstyle {(0)}}$
, exhibit a local maximum at each axial elevation,
$z$
. The locus of these maxima defines the vortex core radius,
$r_{{c}}(z)$
, which delineates the boundary of the so-called forced vortex region. The surface of the latter can be determined by setting
${\partial } u_\theta /{\partial } r = 0$
, as detailed for the mantle in § 4.2.2 and Appendix D. The prescribing equations for the leading-order core radius,
$r_c^{\scriptscriptstyle {(0)}}$
, and its first-order correction,
$r_c^{\scriptscriptstyle {(1)}}$
, are
\begin{equation} {\partial } u_\theta ^{\scriptscriptstyle {(0)}}/{\partial } r = 0 \qquad \text{and} \qquad {\partial } u_\theta ^{\scriptscriptstyle {(1)}}/{\partial } r + r_{{c}}^{\scriptscriptstyle {(1)}} {\partial }^2 u_\theta ^{\scriptscriptstyle {(0)}}/{\partial } r^2 = 0, \end{equation}
where both expressions are evaluated at
$r = r_{{c}}^{\scriptscriptstyle {(0)}}(z)$
. Forthwith, one finds that the maximum angular speed for the incompressible motion occurs in the outlet plane, with its magnitude and location being provided in table 1.
Azimuthal (swirling) velocity profiles in (a) leading-order and (c) first-order isocontours. Radial variations are shown using (b) incompressible and (d) compressible velocity components at
$z = 0$
(
),
$0.25$
(
),
$0.50$
(
) and
$0.75$
(
).
Naturally, an asymptotic expression for
$r_{{c}}^{\scriptscriptstyle {(0)}}(z)$
may be constructed following the same methodology used to obtain (4.5) for
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
. Specifically, the condition
${\partial } u_\theta ^{\scriptscriptstyle {(0)}} / {\partial } r = 0$
, when written explicitly, leads to a transcendental relation of the form
\begin{equation} \bigl [ \bigl ( z / r_{{c}}^{\scriptscriptstyle {(0)}} \bigr ){}^2 - 2 \bigr ]\bigl [ \sin \bigl ( \lambda _1 R_{{c}}^{\scriptscriptstyle {(0)}} \bigr ) / \bigl ( \lambda _1 R_{{c}}^{\scriptscriptstyle {(0)}} \bigr ) - \cos \bigl ( \lambda _1 R_{{c}}^{\scriptscriptstyle {(0)}} \bigr ) \bigr ] + \lambda _1 R_{{c}}^{\scriptscriptstyle {(0)}}\sin \bigl ( \lambda _1 R_{{c}}^{\scriptscriptstyle {(0)}} \bigr ) = 0, \end{equation}
where
$R_{{c}}^{\scriptscriptstyle {(0)}} = \sqrt {\bigl ( r_{{c}}^{\scriptscriptstyle {(0)}} \bigr ){}^2 + z^2}$
. In the absence of an exact solution to (4.9), an approximation for the leading-order vortex core radius,
$r_{{c}}^{\scriptscriptstyle {(0)}}(z)$
, may be constructed using
\begin{equation} \bigl [ r_{{c}}^{\scriptscriptstyle {(0)}}(z) \bigr ]^2 \approx \alpha _0^2\bigl(1-z^2\bigr)\dfrac { 1 - \tfrac {1}{5}\bigl(14 - 3 \alpha _0^{-2}\bigr) z^2 }{ 1 - \tfrac {1}{5}\bigl(27 \alpha _0^2 - 4\bigr) z^2 }, \end{equation}
where
$\alpha _0 \equiv r_{{c}}^{\scriptscriptstyle {(0)}}(0) \approx 0.463251$
. This expression satisfies the same conditions on the residual derivatives at the domain endpoints,
$z=\{ 0,1\}$
, as does the approximation for
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
.
To ascertain its accuracy, the maximum relative error incurred by (4.10) is calculated and found to be
$0.22\,\,\%$
, a comparable value to the
$0.17\,\%$
relative error associated with
$r_{{m}}^{\scriptscriptstyle {(0)}}(z)$
. As further corroborated through figure 9, the analytical expressions remain uniformly valid and practically imperceptible from their numerical solutions over the entire range of interest.
Comparison between analytical solutions (
) and the numerical solution (
$\circ$
) for the vortex core radius,
$r_{{c}}(z)$
. The numerical roots are extracted from (4.9) using a Newton–Raphson solver. These plots display (a) the leading-order distribution,
$r_{{c}}^{\scriptscriptstyle {(0)}}(z)$
, and (b) the first-order compressible correction,
$r_{{c}}^{\scriptscriptstyle {(1)}}(z)$
.
For the reader’s convenience,
$r_{{c}}^{\scriptscriptstyle {(0)}}(z)$
is superimposed onto the streamline plot of figure 3(a), where it appears to be entirely enclosed by the axial mantle. This confirms that the peak swirl velocity will always occur in the downdraft, before crossing the mantle into the toroidal updraft. Separately,
$r_{{c}}^{\scriptscriptstyle {(1)}}(z)$
is recovered through straightforward algebraic manipulation of the second relation in (4.8). Its dependence on
$M_\kappa$
is also illustrated in figure 3(b), where successive expansions of the compressible vortex core may be seen to follow incremental shifts in
$M_\kappa$
. Therein, it may be seen that the expansions of the vortex core radius tend to be consistently more appreciable than those of the mantle for the same excursions in
$M_\kappa$
. This may be attributed to
$r_{{c}}^{\scriptscriptstyle {(1)}}(z)$
being slightly more sensitive to
$M_\kappa$
than
$r_{{m}}^{\scriptscriptstyle {(1)}}(z)$
.
In the outlet plane, the first-order correction,
$\alpha _1 \equiv r_{{c}}^{\scriptscriptstyle {(1)}}(0)$
, can be determined analogously to the mantle radius,
$\beta _1$
. In the present case, we recover
$\alpha _1 \approx 13.7169\kappa ^2$
. As such, the two-term expansion for the normalised vortex core radius at the outlet,
$\alpha \equiv r_{{c}}(0)$
, turns into
\begin{equation} \alpha = 0.463251 + 13.7169\kappa ^2 M_0^2 + {O}\bigl ( M_0^4 \bigr ). \end{equation}
Another hallmark of the present solution for
$u_{\theta }$
consists of its compliance with the physical requirement for a vanishing
$u_\theta$
along the axis of symmetry. This particular characteristic precludes the emergence of a singularity at
$r=0$
, a deficiency that often accompanies inviscid vortex models. Unlike the strictly incompressible model of Williams & Majdalani (Reference Williams and Majdalani2021a
), where imposing a finite
$u_{\theta }$
at the wall engenders an infinite swirl velocity at the axis, our formulation circumvents such behaviour. Presently, enforcing a non-zero circulation constant, such as setting
$C_0 = 1$
in
$u_{\theta }R\sin {\phi } = \pm (C_{B{}}^2\psi ^2 + C_{0}^2)^{1/2}$
, is deliberately avoided. Given that both centreline and impervious wall boundaries belong to the same streamline, maintaining free-vortex behaviour of the
$u_\theta \sim 1/r$
form at the wall inevitably precipitates an undesirable singularity at
$r=0$
. By enforcing
$u_\theta =0$
at the centreline through a judicious choice of
$B(\psi )$
, the resulting model gives rise to physically meaningful thermodynamic variables throughout the domain. In fact, Maicke et al. (Reference Maicke, Cecil and Majdalani2017) reach similar conclusions and adopt comparable assumptions in their compressible treatment of a cylindrically bounded cyclonic motion. While free-vortex behaviour is generally admissible in the context of BHE analysis away from regions of non-uniformity, its presence in the compressible RJ expansion can be problematic in the absence of viscous corrections. Unless the latter are incorporated upfront, setting
$C_0 = 0$
proves to be essential lest negative densities are produced.
Proceeding to the first-order correction,
$u_\theta ^{\scriptscriptstyle {(1)}}$
may be examined through its isocontours and radial profiles in figures 8(c) and 8(d). Here too, we identify two distinct regions that remain separated by a quasi-ellipsoidal surface corresponding to
$u_\theta ^{\scriptscriptstyle {(1)}} = 0$
. This surface, which should not be confused with the axial mantle, delineates zones where compressibility either amplifies or attenuates the local swirl velocity relative to the base flow, thereby altering the swirl polarity.
Based on a cursory graphical examination of figure 8(d), it may be seen that compressibility induces a negative correction near the outlet plane and a corresponding reduction in swirl within the core region; conversely, a positive correction serves to speed up the motion while moving radially outward. Further into the chamber,
$u_\theta ^{\scriptscriptstyle {(1)}}$
becomes predominantly positive across most radii. For the reader’s convenience, the locations and magnitudes of the extrema for
$u_\theta ^{\scriptscriptstyle {(1)}}$
are listed in table 1.
Before leaving this section, a systematic recapitulation of the observed trends in
$ u_\theta ^{\scriptscriptstyle {(1)}}$
may be useful to undertake. In the core region and near the outlet plane, the negative values of
$u_\theta ^{\scriptscriptstyle {(1)}}$
show that flow dilatation leads to a deceleration of the local swirl velocity, particularly in areas characterised by intense axial outflow. Such a reduction can be attributed to fluid expansion during axial acceleration in the vicinity of the outlet, thereby reducing the local density. To maintain local conservation of specific angular momentum, the tangential velocity component must diminish. Conversely, the prevalence of positive
$u_\theta ^{\scriptscriptstyle {(1)}}$
values in the outer radial zones near the outlet and across a substantial portion of the upper chamber volume confirms that compressibility can also enhance the swirl velocity relative to the incompressible case. In these regions, flow deceleration or reduced radial inflow can lead to an increase in density or centrifugal forces, with both factors contributing to an enhanced
$u_{\theta }$
. Collectively, these trends show that compressibility induces a spatial redistribution of the swirling motion. At the axis of rotation, where both velocity components vanish, swirl remains unaffected.
4.4. Pressure, density and temperature fields
Having fully characterised the velocity field, we now turn to the pressure, density and temperature distributions. As illustrated in figure 10(a), isocontours of
$p_1$
decrease monotonically from the outer wall to the chamber origin. The trend in these variations appears to inversely mirror those of
$\mathopen \| \boldsymbol{u} \mathclose \|_0$
, which are captured in figure 5(a). The inverse relation between pressure and velocity can be inferred from energy conservation for an isentropic motion, namely, where an increase in kinetic energy necessitates a decrease in enthalpy. This overarching trend is actually supported by the line distributions in figure 10(b); therein, it can be seen that
$p_1$
reaches its minimum at the origin.
Pressure profiles in (a) first-order (
$p_1/\kappa ^2$
) and (c) second-order (
$p_2/\kappa ^4$
) isocontours. Radial variations are shown using (b) first-order and (d) second-order pressure approximations at
$z = 0$
(
),
$0.25$
(
),
$0.50$
(
) and
$0.75$
(
).
In contrast, as one may surmise from figures 10(c) and 10(d), the second-order pressure component,
$p_2$
, exhibits a distinctly dissimilar trend. Here, a region of positive pressure forms at the origin, where
$p_2$
reaches its maximum value, and then remains surrounded by a spherical region of negative pressure corrections. The minimum value of
$p_2$
also occurs in the outlet plane, at a point that lies approximately halfway between the origin and outer wall. This distribution is driven by dilatational forces, wherein centrifugal acceleration in the outlet plane induces localised fluid expansions. To be brief, precise values of the pressure extrema are summarised in table 1.
A closer examination of the solution character suggests that, for typical values of
$\kappa$
and
$M_0$
, the origin-centred peak of
$p_2$
remains dominated by
$p_1$
. The second-order correction seems to counteract the steep pressure gradients introduced by
$p_1$
, thereby diminishing the sharp variations that are characteristic of incompressible cyclonic flow models. The resulting effect of
$p_2$
on the total pressure may be described as a broadening of the radial gradient associated with the incompressible motion and a redistribution of the thermodynamic energy within the chamber. In actuality, the strong coupling among thermodynamic variables produces similar first-order trends in the density and temperature fields. Specifically, the linearised isentropic relations between pressure, density and temperature enforce direct proportionality relations among the first-order contributions:
$\rho _1$
,
$p_1$
and
$T_1$
, in conformance with (3.11).
Particularly in figure 11(a), the first-order density correction,
$\rho _1$
, exhibits a monotonic decrease that may be attributed to fluid expansion as the centreline is approached. This effect is manifested in a radially outward increase in density. The line plots in figure 11(b) further confirm this trend by showing a consistently increasing density differential as the equatorial plane is approached. This differential diminishes, as one would expect, while approaching the stagnation point at the apex.
Density profiles in (a) first-order (
${\rho }_1/\kappa ^2$
) and (c) second-order (
${\rho }_2/\kappa ^4$
) isocontours. Radial variations are shown using (b) first-order and (d) second-order density approximations at
$z = 0$
(
),
$0.25$
(
),
$0.50$
(
) and
$0.75$
(
).
At second order,
$\rho _2$
in figure 11(c) introduces more elaborate variations. Unlike the first-order field, the second-order correction accounts for nonlinear velocity and thermodynamic interactions. The corresponding line plots in figure 11(d) capture localised density deviations that closely resemble those of the second-order pressure field in figure 10(d).
The first-order temperature field,
$T_1$
, shown in figure 12(a), exhibits a predominant decrease due to expansion cooling, a direct consequence of the pressure drop. The radial profiles in figure 12(b) confirm this temperature reduction towards the core. At second order,
$T_2$
in figure 12(c) introduces additional refinements, and these may be ascribed to compressibility-induced redistributions. The line plots in figure 12(d) showcase localised variations that deviate from first-order trends although they consistently follow those of the second-order pressure and density in figures 10(d) and 11(d).
Temperature profiles in (a) first-order (
${T}_1/\kappa ^2$
) and (c) second-order (
${T}_2/\kappa ^4$
) isocontours. Radial variations are shown using (b) first-order and (d) second-order temperature approximations at
$z = 0$
(
),
$0.25$
(
),
$0.50$
(
) and
$0.75$
(
).
4.5. Range of applicability of the asymptotic solution
The foregoing approximations remain contingent upon the asymptotic correction at each order staying small relative to its preceding order. Moreover, one must ensure that the thermodynamic quantities, such as
$p$
,
$\rho$
and
$T$
, remain positive throughout the flow domain. For example, the pressure expansion,
$p \sim 1 + M_0^2 p_1$
, becomes unphysical when the first-order correction drops to sufficiently negative values to the extent of producing
$p \leqslant 0$
; in fact, it may be easily shown that the most negative and therefore critical value of
$p$
occurs at the origin, i.e. where
$(p_1)_{min } \approx -112.2\kappa ^2$
. Ensuring that
$1 + M_0^2 (p_1)_{min } \gt 0$
translates into
$1 - 112.2 \kappa ^2 M_0^2 \gt 0$
or
\begin{equation} M_\kappa \lt \sqrt {1/112.2} \approx 0.0944. \end{equation}
This upper bound delimits the operating range for which the present analysis continues to hold.
4.6. Total Mach number maps
Maps of the total Mach number,
$\mathopen \| \boldsymbol{M} \mathclose \|$
, remain among the most effective visual tools to assess the significance of flow acceleration. To this end, we use figure 13 to illustrate the evolution of
$\mathopen \| \boldsymbol{M} \mathclose \| = M_0 \mathopen \| \boldsymbol{u} \mathclose \|/\sqrt {T}$
, which gauges the relative importance of inertial and dilatational forces, throughout the hemispherical domain. By using regular expansions for the flow speed,
$\mathopen \| \boldsymbol{u} \mathclose \|$
, and temperature,
$T$
, one may recover a two-term asymptotic representation of the scalar Mach number field:
\begin{equation} \begin{aligned} \mathopen \| \boldsymbol{M} \mathclose \| & \sim M_0 \mathopen \| \boldsymbol{u} \mathclose \|_0 + M_0^3\!\left ( \mathopen \| \boldsymbol{u} \mathclose \|_1 - \tfrac {1}{2} \mathopen \| \boldsymbol{u} \mathclose \|_0 \, T_1 \right )\\ & = \mathopen \| \boldsymbol{M} \mathclose \|_0 + M_0^2\mathopen \| \boldsymbol{M} \mathclose \|_1. \end{aligned} \end{equation}
In the above, the leading- and first-order terms can be readily identified as
$\mathopen \| \boldsymbol{M} \mathclose \|_0 = M_0\mathopen \| \boldsymbol{u} \mathclose \|_0$
and
$\mathopen \| \boldsymbol{M} \mathclose \|_1 = M_0\,(\mathopen \| \boldsymbol{u} \mathclose \|_1 - \mathopen \| \boldsymbol{u} \mathclose \|_0 T_1/2 )$
, respectively. At leading order,
$\mathopen \| \boldsymbol{M} \mathclose \|_0$
provides a spatial map of the fluid inertia scaled by
$\kappa M_0$
. Its distribution, shown in figure 13(a), remains directly proportional to the leading-order velocity magnitude,
$\mathopen \| \boldsymbol{u} \mathclose \|_0$
. Indeed, this proportionality arises because, at this order, the effects of variations in the local speed of sound, dominated by the first-order temperature correction,
$T_1$
, prove to be inconsequential. At the outset, the spatial structure of
$\mathopen \| \boldsymbol{M} \mathclose \|_0/(\kappa M_0)$
may be seen to precisely mirror that of
$\mathopen \| \boldsymbol{u} \mathclose \|_0/\kappa$
in figure 5(a). Distinct features, such as stagnation points and peak locations, become directly reflected within the
$\mathopen \| \boldsymbol{M} \mathclose \|_0$
field. For example, the maximum Mach number may be captured at the origin, where
$(\mathopen \| \boldsymbol{M} \mathclose \|_0)_{max} \approx 12.66\kappa M_0$
; this numerical coefficient is determined by the maximum value of the leading-order velocity, as given in table 1.
Total Mach number profiles in (a)
$\mathopen \| \boldsymbol{M} \mathclose \|_0/(\kappa M_0)$
and (c)
$\mathopen \| \boldsymbol{M} \mathclose \|_0/(\kappa ^3 M_0)$
isocontours. Radial variations are shown using (b) first-order and (d) second-order Mach number approximations at
$z = 0$
(
),
$0.25$
(
),
$0.50$
(
) and
$0.75$
(
).
In contrast, the first-order Mach number correction,
$\mathopen \| \boldsymbol{M} \mathclose \|_1$
, exhibits a spatial structure that differs from both
$\mathopen \| \boldsymbol{M} \mathclose \|_0$
and the first-order velocity correction,
$\mathopen \| \boldsymbol{u} \mathclose \|_1$
. As shown in figure 13(c), the largest
$\mathopen \| \boldsymbol{M} \mathclose \|_1$
remains concentrated in the outflow plane along the axis rather than at the radially shifted peak of
$\mathopen \| \boldsymbol{u} \mathclose \|_1$
. This dissimilarity may be attributed to the dependence of the local speed of sound on the temperature field, which directly affects the Mach number. The first-order pressure drop along the centreline,
$p_1$
, which is depicted in figure 10(a), induces a corresponding reduction in the local temperature,
$T_1$
. Since the speed of sound varies with
$T^{1/2}$
, this localised cooling reduces the acoustic speed while amplifying the local Mach number for a given velocity magnitude.
In fact, one may compare
$\mathopen \| \boldsymbol{M} \mathclose \|_1$
with
$\mathopen \| \boldsymbol{u} \mathclose \|_1$
in order to identify the distinct role of compressibility in controlling the Mach number distribution. While
$\mathopen \| \boldsymbol{u} \mathclose \|_1$
captures modifications to the velocity magnitude due to compressibility,
$\mathopen \| \boldsymbol{M} \mathclose \|_1$
incorporates both velocity and thermodynamic effects in such a seamless way to produce a more complete spatial representation.
4.7. Vorticity field
In this section, the spatial variations of the total vorticity as well as its most dominant component,
$\omega _z$
, are evaluated using both leading- and first-order approximations. In view of the vector parallelism between vorticity and velocity fields, the behaviour of
$\omega _r$
and
$\omega _\theta$
is omitted.
To illustrate the vorticity’s global structure and local variations within the hemispherical domain, figures 14 and 15 are used to present isocontours side-by-side with radial line plots at the four equidistant elevations of
$z = 0$
,
$0.25$
,
$0.50$
and
$0.75$
.
Total vorticity magnitudes in (a) leading-order and (c) first-order isocontours. Radial variations are shown using (b) incompressible and (d) compressible vorticity components at
$z = 0$
(
),
$0.25$
(
),
$0.50$
(
) and
$0.75$
(
).
Axial vorticity in (a) leading-order (
$\omega _{z}^{(0)}$
) and (c) first-order (
$\omega _{z}^{(1)}$
) isocontours. Radial variations are shown using (b) incompressible and (d) compressible vorticity components at
$z = 0$
(
),
$0.25$
(
),
$0.50$
(
) and
$0.75$
(
).
In accordance with the defining Beltramian criterion (
$k\boldsymbol{u} = {\boldsymbol{\nabla }}\times \boldsymbol{u}$
),
$\boldsymbol{\omega }$
and
$\boldsymbol{u}$
remain linearly related through
\begin{equation} \boldsymbol{\omega } = {\boldsymbol{\nabla }} \times \boldsymbol{u} = - C_B \rho \boldsymbol{u}. \end{equation}
Because the scalar factor
$- C_B \rho$
varies spatially with the density, the flow differs fundamentally from a Trkalian field, wherein a constant scalar factor is required (Maicke et al. Reference Maicke, Cecil and Majdalani2017). Expanding (4.14) in orders of
$M_0^2$
confirms that the leading-order vorticity vector,
$\boldsymbol{\omega }_0$
, remains linearly proportional to the leading-order velocity through
$\boldsymbol{\omega }_0 = -\lambda _1 \boldsymbol{u}_0$
. Consequently, the spatial structure of each leading-order vorticity component may be seen to directly mirror that of its velocity counterpart. Since the velocity field is already detailed in § 4.3, the discussion of the leading-order vorticity will be abbreviated. In contrast, first-order density variations induce a structure in
$\boldsymbol{\omega }_1$
that remains wholly dissimilar from that of
$\boldsymbol{u}_1$
. As such, our forthcoming analysis will focus on the distinct characteristics of
$\boldsymbol{\omega }_1$
, which embody the vast majority of dilatational effects.
4.7.1. Total vorticity
The vorticity magnitude,
$\mathopen \| \boldsymbol{\omega } \mathclose \|$
, can be synthesised using an analogous decomposition to that of
$\mathopen \| \boldsymbol{u} \mathclose \|$
in (4.7). The corresponding leading- and first-order contributions yield
\begin{equation} \mathopen \| \boldsymbol{\omega } \mathclose \|_0 = \lambda _1 \mathopen \| \boldsymbol{u} \mathclose \|_0 \qquad \text{and} \qquad \mathopen \| \boldsymbol{\omega } \mathclose \|_1 = \lambda _1 \mathopen \| \boldsymbol{u} \mathclose \|_1 + \bigl ( C_B^{\scriptscriptstyle {(1)}} + \lambda _1 \rho _1 \bigr )\mathopen \| \boldsymbol{u} \mathclose \|_0, \end{equation}
where
$\mathopen \| \boldsymbol{u} \mathclose \|_1 = \boldsymbol{u}_1 \boldsymbol{\cdot }\hat {\boldsymbol{u}}_0$
denotes the first-order correction to the flow speed along leading-order streamlines. In what follows, the spatial distributions of both
$\mathopen \| \boldsymbol{\omega } \mathclose \|_0$
and
$\mathopen \| \boldsymbol{\omega } \mathclose \|_1$
are characterised by describing their spatial distributions in figure 14.
The leading-order term,
$\mathopen \| \boldsymbol{\omega } \mathclose \|_0$
, reflects the vorticity of the base flow. As per (4.15), since it remains directly proportional to
$\mathopen \| \boldsymbol{u} \mathclose \|_0$
, its spatial distributions in figures 14(a) and 14(b) mirror those of the velocity magnitude in figure 5(a). Here, we find
$\mathopen \| \boldsymbol{\omega } \mathclose \|_0$
to be strictly positive with a peak intensity at the origin that diminishes monotonically outward. Its maximum value may be evaluated to be
\begin{equation} (\mathopen \| \boldsymbol{\omega } \mathclose \|_0)_{max} = \frac {2}{3}\frac {\kappa \lambda _1^2}{\beta _0 \, \mathinner {{j}_{1}(\lambda _1\beta _0)}} \approx 56.89\kappa . \end{equation}
On the one hand, and following the behaviour of
$\mathopen \| \boldsymbol{u} \mathclose \|_0$
, the base vorticity can be seen to decrease away from the origin while approaching the hemispherical boundary. The underlying spatial decay may be further inferred from the radial line plots of figure 14(b); therein, the
$z = 0$
profile captures the highest magnitudes as they depreciate radially outward, while profiles at higher elevations exhibit progressively reduced magnitudes and flatter distributions with respect to
$r$
.
On the other hand, the first-order correction,
$\mathopen \| \boldsymbol{\omega } \mathclose \|_1$
, which is normalised by
$\kappa ^3$
in figures 14(c) and 14(d), quantifies the change in vorticity due to first-order dilatational effects. Unlike
$\mathopen \| \boldsymbol{\omega } \mathclose \|_0$
,
$\mathopen \| \boldsymbol{\omega } \mathclose \|_1$
displays both positive and negative contributions that can locally augment or suppress the vorticity magnitude. In fact,
$\mathopen \| \boldsymbol{\omega } \mathclose \|_1$
may be viewed as the superposition of two separate contributions according to (4.15).
The first contribution,
$\lambda _1 \mathopen \| \boldsymbol{u} \mathclose \|_1$
, stems from the compressibility-induced shift in flow speed along the base streamlines,
$\mathopen \| \boldsymbol{u} \mathclose \|_1$
, multiplied by the leading-order eigenvalue,
$\lambda _1$
. As seen in figures 5(c) and 5(d),
$\mathopen \| \boldsymbol{u} \mathclose \|_1$
remains predominantly positive, especially around the centreline in the outlet plane, where it serves to speed up the motion. As such, it appears to contribute favourably to
$\mathopen \| \boldsymbol{\omega } \mathclose \|_1$
in regions where streamline acceleration is prominent. The second contribution,
$(C_B^{\scriptscriptstyle {(1)}} + \lambda _1 \rho _1)\mathopen \| \boldsymbol{u} \mathclose \|_0$
, captures the combined effects of the first-order correction to the eigenvalue,
$C_B^{\scriptscriptstyle {(1)}}$
, and the first-order density perturbation,
$\lambda _1 \rho _1$
, which jointly act upon
$\mathopen \| \boldsymbol{u} \mathclose \|_0$
.
Based on figures 11(a) and 11(b),
$\rho _1$
remains negative throughout the domain, thus leading to fluid expansion that is particularly appreciable near the origin. This expansion tends to produce negative contributions while the role of
$C_B^{\scriptscriptstyle {(1)}}$
, which depends on the first-order thermodynamic field, only serves to further enhance this term. The overall sign and magnitude of the second term depend on the balance between the density perturbation and
$C_B^{\scriptscriptstyle {(1)}}$
, all multiplied by the local base flow speed,
$\mathopen \| \boldsymbol{u} \mathclose \|_0$
. The resulting spatial distribution of
$\mathopen \| \boldsymbol{\omega } \mathclose \|_1$
, shown in figure 14(c), illustrates the outcome of this coupling.
At first glance, a large region of negative
$\mathopen \| \boldsymbol{\omega } \mathclose \|_1$
may be seen to dominate the lower hemisphere, where a net reduction in vorticity may be caused by dilatational effects. Therein, the lowest value,
$(\mathopen \| \boldsymbol{\omega } \mathclose \|_1)_{min} \approx -2100 \kappa ^3$
, appears at the origin. The observed depreciation occurs because the second term,
$(C_B^{\scriptscriptstyle {(1)}} + \lambda _1 \rho _1)\mathopen \| \boldsymbol{u} \mathclose \|_0$
, yields a substantially negative value near the origin; since
$\rho _1$
reaches its lowest negative value locally, its product with
$\mathopen \| \boldsymbol{u} \mathclose \|_0$
, which reaches its peak value at
$r=0$
, tends to eclipse the relatively small contribution from
$\lambda _1 \mathopen \| \boldsymbol{u} \mathclose \|_1$
in the same neighbourhood.
In contrast, a region of positive
$\mathopen \| \boldsymbol{\omega } \mathclose \|_1$
may be observed in the toroidal segment of the domain, which is bounded by the wall and lying outside of the elliptical surface defined by
$(37r/14)^2 + (86z/29)^2 = 1$
. This region of positive
$\mathopen \| \boldsymbol{\omega } \mathclose \|_1$
corresponds approximately to the area of maximum positive streamline acceleration in figure 5(c). Within this toroidal segment, the positive contribution from the first term,
$\lambda _1 \mathopen \| \boldsymbol{u} \mathclose \|_1$
, exceeds the opposing effects stemming from the second term,
$(C_B^{\scriptscriptstyle {(1)}} + \lambda _1 \rho _1)\mathopen \| \boldsymbol{u} \mathclose \|_0$
. The radial line plots presented in figure 14(d) serve as a quantitative confirmation of this structure.
4.7.2. Axial vorticity profile
Next, we examine the axial component of vorticity,
$\omega _z$
, which represents the local rate of fluid rotation about the
$z$
axis. Its leading- and first-order approximations are characterised in figure 15. As expected, the leading-order axial vorticity stays proportional to the negative of the axial velocity,
$\omega _z^{\scriptscriptstyle {(0)}} = -\lambda _1 u_z^{\scriptscriptstyle {(0)}}$
. Since
$u_z^{\scriptscriptstyle {(0)}}$
switches polarity across the mantle surface, as per figure 7(a),
$\omega _z^{\scriptscriptstyle {(0)}}$
undergoes a sign reversal. Specifically,
$\omega _z^{\scriptscriptstyle {(0)}}$
remains positive in the inner downdraft region, where
$u_z^{\scriptscriptstyle {(0)}} \lt 0$
, and negative in the outer updraft region, where
$u_z^{\scriptscriptstyle {(0)}} \gt 0$
. The corresponding isocontour and line graphs in figures 15(a) and 15(b) clearly depict these patterns. Therein, the magnitude of
$\omega _z^{\scriptscriptstyle {(0)}}$
can be seen to reach its peak value at the origin, where
$(\omega _z^{\scriptscriptstyle {(0)}})_{max } \approx 56.89 \kappa$
, and then decreases radially outward.
As for the spatial distribution of the first-order correction,
$\omega _z^{\scriptscriptstyle {(1)}}$
, it may be inferred from figures 15(c) and 15(d). Here too,
$\omega _z^{\scriptscriptstyle {(1)}}$
can be reconstructed from two additive terms that form
\begin{equation} \omega _z^{\scriptscriptstyle {(1)}} = -\bigl ( C_B^{\scriptscriptstyle {(1)}} + \lambda _1 \rho _1 \bigr ) u_z^{\scriptscriptstyle {(0)}} - \lambda _1 u_z^{\scriptscriptstyle {(1)}}. \end{equation}
The first member,
$-(C_B^{\scriptscriptstyle {(1)}} + \lambda _1 \rho _1) u_z^{\scriptscriptstyle {(0)}}$
, superimposes the constant first-order eigenvalue correction,
$C_B^{\scriptscriptstyle {(1)}}$
, and the spatially varying first-order density perturbation,
$\lambda _1 \rho _1$
, both acting upon
$u_z^{\scriptscriptstyle {(0)}}$
. In the presence of a negative
$\rho _1$
, the term
$-\lambda _1 \rho _1 u_z^{\scriptscriptstyle {(0)}}$
contributes negatively in the inner downdraft region, where
$u_z^{\scriptscriptstyle {(0)}} \lt 0$
, and positively in the outer updraft, where
$u_z^{\scriptscriptstyle {(0)}} \gt 0$
. The term
$-C_B^{\scriptscriptstyle {(1)}} u_z^{\scriptscriptstyle {(0)}}$
adds to this effect, with its sign depending on the value of
$C_B^{\scriptscriptstyle {(1)}}$
. From a practical perspective, the addition of these two terms helps to clarify how fluid expansion, paired with the base helical motion and eigenvalue correction, can alter the axial vorticity. As for the second member,
$-\lambda _1 u_z^{\scriptscriptstyle {(1)}}$
, it contributes negatively in the outer toroidal vortex, where
$u_z^{\scriptscriptstyle {(1)}}\gt 0$
, and positively across the inner core vortex, where
$u_z^{\scriptscriptstyle {(1)}}\lt 0$
.
The superposition of these two principal contributions determines the final structure of
$\omega _z^{\scriptscriptstyle {(1)}}$
. In practice, their sum leads to the complex multi-lobed pattern observed in figure 15(c). Therein, a region of intense negative axial vorticity surrounds the origin, where the lowest value of
$\omega _z^{\scriptscriptstyle {(1)}} \approx -2097\kappa ^3$
is realised. Such a feature may be attributed to the large negative contribution from the first term in the downdraft, which outweighs the positive contribution from the second term,
$-\lambda _1 u_z^{\scriptscriptstyle {(1)}}$
. Radially outward, particularly near the mid-radius and mid-elevation, distinct positive lobes emerge. These develop where the first term provides a net positive contribution, especially in the updraft region, and combines constructively with the second term, depending on the local sign of
$u_z^{\scriptscriptstyle {(1)}}$
. The radial variations of the latter, which are provided in figure 15(d), help to discern the transitional points separating positive and negative net contributions at different elevations.
4.8. Helicity density and global helicity
Another characteristic quantity of interest consists of the total helicity,
$\mathcal{H} = \int _{\mathcal{V}} h \,{\rm d}{V}$
, where
$h \equiv \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\omega }$
denotes the helicity density (Moffatt & Tsinober Reference Moffatt and Tsinober1992; Sharma & Sameen Reference Sharma and Sameen2020). While the volume integral serves as a measure of the global knottedness of vortex lines throughout the domain
$\mathcal{V}$
,
$h$
quantifies the degree of local alignment between the velocity and vorticity vectors; particularly,
$h$
reaches its peak value when the direction of a fluid particle’s rotational axis coincides precisely with that of its translational motion. As elegantly described by Moffatt (Reference Moffatt1969), the sign of
$h$
determines the handedness of the local helical structures. For the Beltramian motion under consideration, where
$\boldsymbol{\omega } = -C_B \rho \boldsymbol{u}$
, the helicity density becomes directly proportional to the local kinetic energy density, given that it reduces to
$h = -C_B \rho (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{u})$
.
In conformance with other variable expansions,
$h$
can be asymptotically approximated as a series in
$M_0^2$
; this may be accomplished using
$h = h_0 + M_0^2 h_1 + {O}(M_0^4)$
, where the leading- and first-order terms can be gathered into
\begin{equation} h_0 = -C_B^{\scriptscriptstyle {(0)}}\boldsymbol{u}_0 \boldsymbol{\cdot }\boldsymbol{u}_0 \quad \text{and} \quad h_1 = - C_B^{\scriptscriptstyle {(0)}}\bigl ( \rho _1 \boldsymbol{u}_0 \boldsymbol{\cdot }\boldsymbol{u}_0 + 2 \boldsymbol{u}_0 \boldsymbol{\cdot }\boldsymbol{u}_1 \bigr ) - C_B^{\scriptscriptstyle {(1)}} \boldsymbol{u}_0 \boldsymbol{\cdot }\boldsymbol{u}_0. \end{equation}
These relations enable us to evaluate both the leading-order helicity density and its first-order compressible correction. They also permit a global assessment of the total helicity.
4.8.1. Leading-order helicity
Recalling that
$C_B^{\scriptscriptstyle {(0)}} = \lambda _1$
, the leading-order helicity density can be written as
$h_0 = -\lambda _1 \mathopen \| \boldsymbol{u}_0 \mathclose \|^2$
. Moreover, since
$\lambda _1 \gt 0$
and
$\mathopen \| \boldsymbol{u}_0 \mathclose \|^2 \geqslant 0$
,
$h_0$
becomes necessarily non-positive throughout the domain. This consistent sign leads to a uniformly handed incompressible vortex. Specifically for
$h_0 \leqslant 0$
, the incompressible base flow will mainly exhibit left-handed helicity.
Interestingly, the spatial distribution of
$h_0$
, which is showcased in figure 16(a), directly mirrors that of
$-\mathopen \| \boldsymbol{u}_0 \mathclose \|^2$
, captured earlier through figure 5(a). Therein, the leading-order helicity density reaches its most negative value at the origin, where the peak of
$\mathopen \| \boldsymbol{u}_0 \mathclose \|$
occurs. As for its absolute magnitude, it diminishes rather monotonically while approaching the curved wall; it then vanishes at the apex due to flow bending away from the centreline.
Spatial distribution of the helicity density using (a) leading-order (
$h_0$
) and (b) first-order (
$h_1$
) isocontours.
Next, we may proceed to evaluate
$\mathcal{H}_0$
, the total helicity at the leading order. This may be accomplished by integrating
$h_0$
over the chamber volume, namely, by taking
\begin{equation} \mathcal{H}_0 = \int _{\mathcal{V}} h_0\, {\rm d}{V} = -\lambda _1 \int _{\mathcal{V}} \mathopen \| \boldsymbol{u}_0 \mathclose \|^2 \,{\rm d}{V}. \end{equation}
Using the Beltramian solution for
$\boldsymbol{u}_0$
in (3.30), this integral can be numerically evaluated to obtain
$\mathcal{H}_0 \approx -320.4\kappa ^2$
. This expression enables us to reaffirm the strong quadratic dependence of the total helicity on the inflow parameter. In fact, the same outcome may be achieved using an alternative transformation that is outlined in Appendix E.
4.8.2. First-order helicity correction
The first-order correction to the helicity density,
$h_1$
, which is prescribed by (4.18) at
${O}(M_0^2)$
, helps to quantify deviations from the local flow topology due to dilatational effects. This correction may be viewed as a byproduct of the first-order density perturbations, flow velocity fluctuations and corrections to the Beltramian eigenvalues, which remain controlled in large part by the base flow velocity,
$\boldsymbol{u}_0$
. Specifically,
$h_1$
incorporates a
$-C_B^{\scriptscriptstyle {(0)}}\rho _1 \boldsymbol{u}_0\boldsymbol{\cdot }\boldsymbol{u}_0$
product that captures the effect of fluid dilatation on the helicity associated with the base flow kinetic energy. As such, regions of expansion where negative
$\rho _1$
values yield a positive contribution tend to locally disperse the magnitude of the left-handed baseline helicity. Another contributor,
$-2 C_B^{\scriptscriptstyle {(0)}}\boldsymbol{u}_0 \boldsymbol{\cdot }\boldsymbol{u}_1$
, represents the change in helicity due to the projection of the first-order velocity onto
$\boldsymbol{u}_0$
. Since the first-order component,
$\boldsymbol{u}_1$
, stays parallel to the leading-order velocity during flow acceleration along the base streamlines, this second term generates a negative contribution. The third and final term,
$-C_B^{\scriptscriptstyle {(1)}}\boldsymbol{u}_0\boldsymbol{\cdot }\boldsymbol{u}_0$
, may be seen to account for the first-order shift in the proportionality constant,
$C_B^{\scriptscriptstyle {(1)}}$
, which multiplies the Beltramian flow kinetic energy.
By way of illustration, the volumetric distribution of the normalised helicity correction,
$h_1/\kappa ^4$
, is showcased in figure 16(b). Therein, it may be seen that the first-order correction undergoes appreciable spatial variations that span regions of both positive and negative helicity. The strongly positive region that forms around the origin may be viewed as being primarily modulated by
$\rho _1$
, whose spatial expansion in the core region remains commensurate with the locally accelerating motion. Conversely, the stacked layers of negative
$h_1/\kappa ^4$
that surround the core region may be attributed to the interplay causing the velocity and eigenvalue correction terms to more than offset the local density effect.
4.8.3. Total helicity correction
The first-order correction to the total helicity can be readily ascertained by integrating
$h_1$
over the hemispherical volume
$\mathcal{V}$
. Using
$\boldsymbol{u}_0,\ \boldsymbol{u}_1$
and
$\rho _1$
from (3.30), (3.51) and (3.34a
), as well as the available constants,
$C_B^{\scriptscriptstyle {(0)}}$
and
$C_B^{\scriptscriptstyle {(1)}}$
, one may compute
$\mathcal{H}_1 = \int _{\mathcal{V}} h_1 \,{\rm d}{V}$
to obtain
$\mathcal{H}_1 = -19041\kappa ^4$
. In fact, the quartic dependence of
$\mathcal{H}_1$
on the inflow parameter may be expected given the quadratic dependence of
$\mathcal{H}_0$
on
$\kappa$
. This confirms that the first-order effects, encapsulated by the
$M_0^2$
terms in the RJ expansion, serve to modify the helicity density within the hemispherical chamber, depending on the size of
$\kappa$
. In practice, the positive and negative contributions in figure 16(b) can be shown to prescribe the integrated helicity over the entire domain.
In fact, this behaviour can be further justified by recognising that
$h_1$
represents the first-order term in the expansion of
$h = \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\omega }$
. As such,
$h_1=\boldsymbol{u}_0 \boldsymbol{\cdot }\boldsymbol{\omega }_1 + \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\omega }_0$
can be decomposed into two parts that may be readily integrated. In this effort, the integral
$\mathcal{H}_1 = \int _{\mathcal{V}} (\boldsymbol{u}_0 \boldsymbol{\cdot }\boldsymbol{\omega }_1 + \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\omega }_0) \,{\rm d}{V}$
denotes the first-order variation in total helicity. Physically, one may infer that compressibility at
${O}(M_0^2)$
helps to reorganise the internal structure and twisting of vortex lines in a manner to slightly distort the topological signature of the baseline helicity,
$\mathcal{H}_0$
. For added clarity, an alternative method to evaluate
$\mathcal{H}_1$
is provided in Appendix E.
4.8.4. Implications of helicity distribution
The spatial distribution of helicity density provides useful insights into the flow’s structural organisation. Regions of large helicity density correspond to areas where flow twisting exhibits strong, coherent vortical structures that retain a firm alignment of velocity and vorticity. The consistently negative
$h_0$
corresponds to a predominantly left-handed helicity for the incompressible motion. The first-order correction,
$h_1$
, may be seen to prescribe how compressibility perturbs this inherent handedness and linkage of vortex lines. The emergence of positive
$h_1$
regions marks areas where compressibility acts to locally reduce the intensity of this left-handedness or, alternatively, induce a right-handed helical motion for cases where
$|h_1|$
sufficiently exceeds
$|h_0|$
to alter the sign of the net helicity,
$h_0 +M_0^2h_1$
.
4.9. Implications for cyclonically driven chambers
From a broader perspective, and based on the foregoing analysis of cyclonic flow features, it may be helpful to evaluate some of the characteristics that may be conveyable to geometrically and kinematically similar vortex engines.
Firstly, this work establishes that increasing the compressibility parameter,
$M_\kappa$
, results in an outward displacement of the mantle surface with a volumetric expansion of the central downdraft region relative to that of the annular updraft. The corresponding change in volumetric distributions affects the residence time of the core reactants and, in turn, the combustion efficiency. When considering the class of ‘cold-wall’ vortex engines, the hot reactants are deliberately confined to the central downdraft region, while the outer updraft is relied upon to insulate the walls by virtue of a rapidly convecting stream of a low-temperature oxidiser (Chiaverini et al. Reference Chiaverini, Malecki, Sauer, Knuth, Gramer and Majdalani2003; Maicke et al. Reference Maicke, Cecil and Majdalani2017). The expansion of the core region will, in this context, shift the volumetric ratio between the centrally located reaction zone and the wall-insulating outer zone. Such a shift can affect the mixing rate, mixture residence time, thermal loading on the wall and overall combustion efficiency (Lilley Reference Lilley1977). It is hoped that these various factors, including the accurate prediction of the mantle interface, can be used to inform the development of more compact and cost-effective designs.
Secondly, alongside the outward mantle excursions, the position of the peak swirl velocity, which determines the size of the core vortex, is also displaced radially outward (as per table 1). In traditional swirl-stabilised combustors, the diameter of the core vortex delimits the primary flame-anchoring zone. In a vortex chamber, the successive expansions of the core at increasing Mach numbers serve to displace the peak swirl velocity along with the principal reaction zone away from the chamber centreline and closer to the wall. Such volumetric distortions can have an appreciable bearing on thermal stress calculations as well as the proper orientation of the headwall injectors that are typically employed in vortex engines.
Thirdly, we have seen that increasing
$M_\kappa$
leads to the spatial steepening of the axial velocity profile. While approaching the exit plane, for instance, the central downdraft becomes less elongated. This may be attributed to the peak negative axial velocity being displaced radially outward from the chamber axis (see figure 7). Since the onset of vortex breakdown remains highly sensitive to the radial distribution of the axial and swirl velocity components (Syred & Beér Reference Syred and Beér1974), it can be speculated that the flattening of the axial velocity profile may have an adverse effect on flow stability and the onset of oscillatory motions in actual vortex engines. A dedicated study of the stability characteristics of this helical motion may be warranted to systematically investigate its response to small perturbations.
Fourthly, the outer expansion of the mantle with increasing dilatational levels may be helpful to consider at the conceptual design stage of vortex engines, namely, to allow for a seamless and unobstructed outflow into the nozzle entrance region. In practice, the suppression of stagnation or reversed-flow pockets will be essential for achieving stable engine performance.
Lastly, from a procedural standpoint, variations in inlet parameters, including the nozzle shape and size, may be systematically tailored to suppress the emergence of acoustic instabilities and multiple flow reversals. In this process, the non-dimensional grouping recovered here, which combines the injection Mach number and inflow parameter, may perhaps be leveraged as a practical guide for assessing dilatational contributions across a broad range of chamber configurations and mass inflow rates.
4.10. A spherically cyclonic Beltramian vortex
Before closing, it may be instructive to note that the leading-order streamfunction,
$\psi _0$
, given by (3.29), can be conceptually extended to describe the flow evolution either within or around a wholly spherical domain. The corresponding velocity components, given by (3.30), comprise spherical Bessel functions and powers of
$\sin \phi$
. As such, they remain mathematically consistent with a domain where the polar angle
$\phi$
ranges from
$0$
to
$\pi$
. In this section, we briefly explore the application of the present solution to a spherical configuration.
The flow evolving inside a sphere, which is illustrated in figure 17(a), exhibits a streamline pattern that is akin to that of Hill’s spherical vortex (Hill Reference Hill1894). It consists of two counter-rotating vortices that remain symmetrically positioned about the vertical axis. These vortices comprise closed-loop streamlines that spiral inward towards their centres. Each vortex shows circulatory motion with opposing rotational directions that lead to a balanced dipolar structure. Inside this spherical domain, an axial mantle that separates regions of upward and downward axial flow is depicted using a broken line. A chained line is also used to track the evolution of the peak tangential speed,
$(u_\theta )_{max}$
. The latter remains engulfed by the mantle surface.
Display of (a) internal streamlines and loci of
$u_z=0$
(axial mantle, 
) and
$(u_\theta )_{max}=0$
(vortex core radius, 
) for a spherically cyclonic Beltramian vortex; we also show in (b) the external streamlines from the same solution around a sphere.
As one may infer from the graph, the streamline patterns combine poloidal circulation in the meridional plane with a toroidal rotation in the azimuthal plane. They remain smooth and continuous with high levels of vorticity near the vortex cores. This dual rotation gives rise to a self-sustaining, spherically confined cyclonic motion that satisfies the problem’s fundamental boundary conditions. These include the absence of radial cross-flow at the wall, with
$u_{\!R} = 0$
at
$R = 1$
, as well as symmetry at the poles, with
$u_\phi = 0$
at
$\phi = \{0,\pi \}$
. Once established, the flow appears to retain its stability in the absence of viscosity while behaving as a spherically confined ‘perpetual fluid flywheel’.
Kinematically, the dipolar vortex structure that evolves may be compared with Hill’s spherical vortex, which exhibits a single vortex core (Hill Reference Hill1894), as well as the Lamb–Chaplygin dipole, which is accompanied by two-dimensional counter-rotating vortices (Meleshko & van Heijst Reference Meleshko and van Heijst1994). Unlike the Taylor–Proudman motion (Jacobs Reference Jacobs1964) or the rotating spherical Couette flow (Munson & Joseph Reference Munson and Joseph1971), which entail columnar structures or viscous effects, the present solution does not require the forced rotation of the entire sphere because of its allowance for wall slippage.
As for the external motion around the sphere, which is showcased in figure 17(b), it leads to an axisymmetric Beltramian vortex enveloping a spherical void and featuring a stagnation ring in the equatorial plane. In contrast to the irrotational potential motion described by Williams & Majdalani (Reference Williams and Majdalani2021b ), the present streamlines consist of a collection of nested, closed-loop stream surfaces wrapping around the sphere, namely, with recirculation lobes in the sphere’s lateral flanks. The resulting motion may represent vortex structures forming around spherical obstacles, such as those arising in geophysical or astrophysical research, including the motion of confined planetary cores (Ardes, Busse & Wicht Reference Ardes, Busse and Wicht1997; Dormy et al. Reference Dormy, Soward, Jones, Jault and Cardin2004).
5. Conclusion
This work presents a comprehensive analytical framework for describing the steady, compressible, cyclonic motion of a Beltramian profile in a hemispherically bounded chamber. This particular flow configuration is motivated by a vortex engine concept that is developed by Sierra Space Corporation for the Dream Chaser spaceplane. The analysis relies on an asymptotic RJ expansion of the compressible BHE in conjunction with basic conservation principles. This is accomplished using spherical coordinates and expansions in terms of the squared injection Mach number,
$M_{0}^{2}$
.
Our main findings include a complete set of closed-form expressions describing the compressible Beltramian flow field, with each quantity being resolved to the first or second order in
$M_0^2$
. The leading-order solution identically recovers a known incompressible Beltramian profile as a limiting process verification for
$M_0 \to 0$
. The first-order corrections, which are numerically verified, yield explicit approximations for all velocity components, thermodynamic properties, vorticity and helicity. They also provide compact analytical relations for the spatial deformations of the mantle and vortex core surfaces, namely, by capturing their outward dilatational displacements. In this vein, the fractional radii of the mantle,
$\beta$
, and the vortex core,
$\alpha$
, are found to expand according to
$\beta \approx 0.6106 + 8.612\kappa ^2 M_0^2$
and
$\alpha \approx 0.4633 + 13.717\kappa ^2 M_0^2$
, respectively.
Another helpful outcome of the asymptotic expansion is, perhaps, the identification of a single compressibility parameter,
$M_\kappa \equiv \kappa M_0$
, which consolidates the inlet Mach number,
$M_0$
, and the Ekman-type inflow parameter,
$\kappa$
. This grouping may be seen to encapsulate all compressibility sources including the combined effects of varying the mass influx, geometry and speed of sound. So long as
$M_\kappa \leqslant 0.0944$
, which serves as a conservative upper limit that can be difficult to exceed,
$M_\kappa$
may be used as a suitably small perturbation parameter.
Additional results comprise a detailed vorticity analysis to characterise the spatial distribution and orientation of the axial component of vorticity,
$\omega _z$
, as well as the total vorticity contribution. These confirm the Beltramian alignment of the velocity and vorticity vectors at leading order. In this process, the leading-order axial vorticity,
$\omega _z^{\scriptscriptstyle {(0)}}$
, is seen to undergo a polarity switch across the mantle surface. Our findings also suggest that compressibility introduces appreciable first-order contributions to the vorticity distribution. These corrections help to identify areas where vorticity is either amplified or suppressed, particularly near the outlet plane, i.e. where fluid expansion is most prominent. In this context,
$\omega _z^{\scriptscriptstyle {(1)}}$
is seen to exhibit multi-lobed patterns and, in some regions, achieve higher magnitudes than its zeroth-order counterpart.
Moreover, by evaluating the helicity density and total helicity in closed form, a quantification of the topological complexity of the flow is undertaken. For this Beltramian motion, the helicity density,
$h = \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\omega }$
, is found to be directly proportional to the local kinetic energy density,
$h = -C_B \rho \mathopen \|\boldsymbol{u} \mathclose \|^2$
, with the leading-order helicity corresponding to a uniformly left-handed vortex (
$h_0\leq 0$
). On the one hand, the helicity density is seen to peak in the mantle region, where the interaction between the radial and azimuthal velocities is most pronounced. On the other hand, integrating this quantity across the entire volume enables us to evaluate the global helicity. The latter serves as a scalar measure of the flow’s knottedness and dynamic coherence. Through systematic reductions of the helicity integrals, the leading- and first-order corrections to the total helicity are found to be proportional to
$\kappa ^2$
and
$\kappa ^4$
, respectively. As such, dilatational effects appear to reorganise and redistribute the internal vortical structures by altering the global topological knottedness of the baseline motion depending on the size of
$\kappa ^4$
.
Acknowledgements
We are deeply indebted to R.J. Hartfield, Jr. and V. Ahuja at Altair Engineering, and to M.J. Chiaverini, D. Benner, B. Pomeroy and A. Sauer at Sierra Space Corporation. We greatly appreciate their continued support of our cyclonic flow investigations.
Funding
This work was supported partly by AFWERX, through contract no. AF-FA8649-24-P-0021 with Sierra Space Corporation and Altair Engineering, and partly by the Hugh and Loeda Francis Chair of Excellence, Department of Aerospace Engineering, Auburn University.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available upon request.
Appendix A. Compressible Bragg–Hawthorne formulation
In what follows, the assumptions that accompany the compressible BHE are revisited and discussed. Given the flow’s rotational symmetry about
$z$
, mass conservation, given by (2.6), can be written as
\begin{equation} \frac {\partial }{\partial R}(R^2 \sin {\phi }\,\rho u_{\!R}) + \frac {\partial }{\partial \phi }(R \sin {\phi } \,\,\rho u_\phi ) = 0. \end{equation}
Equation (A1) can be inherently satisfied by expressing the radial and polar velocity components in terms of streamfunction derivatives. The velocity components become connected to
$\psi$
through
\begin{equation} u_{\!R} = -\frac {1}{R^2\sin {\phi }\,\rho }\frac {\partial \psi }{\partial \phi } \qquad \text{and} \qquad u_\phi = \frac {1}{R\sin {\phi }\,\rho }\frac {\partial \psi }{\partial R}. \end{equation}
Within (A2), the scaling of
$\psi$
follows the convention used to normalise
$\boldsymbol{u}$
,
$\rho$
and
$R$
in (2.5). Proceeding to the treatment of the energy equation, the latter may be combined with continuity to recover
\begin{equation} \boldsymbol{u} \boldsymbol{\cdot }{\boldsymbol{\nabla }} {{H}} = 0. \end{equation}
Equation (A3) establishes an orthogonal relation between
$\boldsymbol{u}$
and the gradient of the stagnation enthalpy,
$\boldsymbol{\nabla }H$
. Physically, this condition implies that
$H$
must remain invariant along streamlines. In practice, (A3) can be identically satisfied by expressing the stagnation enthalpy as a functional composition of
$\psi$
. This may be accomplished by first setting
\begin{equation} {\boldsymbol{\nabla }}(H\circ \psi ) = \frac {{\rm d} H}{{\rm d} \psi }{\boldsymbol{\nabla }}\psi . \end{equation}
One may then substitute (A4) into (A3) to obtain
\begin{equation} \boldsymbol{u}\boldsymbol{\cdot }{\!\left ( \frac {{\rm d} H}{{\rm d} \psi } {\boldsymbol{\nabla }} \psi \right )} =\frac {{\rm d} H}{{\rm d} \psi } \,\boldsymbol{u}\boldsymbol{\cdot }{ {\boldsymbol{\nabla }} \psi } = 0. \end{equation}
Recalling that
$\boldsymbol{u}\boldsymbol{\cdot }{\boldsymbol{\nabla }} \psi = 0$
in any coordinate system, it may be seen that both (A3) and (A5) become identically satisfied. Physically, the underlying orthogonality implies that
$H$
varies only in directions that are perpendicular to the velocity vector at every point in the domain. To further refine our set of governing equations, a re-examination of the azimuthal component of the momentum conservation equation may be warranted. Using
\begin{equation} \boldsymbol{u}\boldsymbol{\cdot }{\boldsymbol{\nabla }} \boldsymbol{u} \equiv \tfrac {1}{2}{\boldsymbol{\nabla }}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}) + \boldsymbol{\omega }\times \boldsymbol{u}, \end{equation}
the azimuthal component of the momentum equation can be expressed as
\begin{equation} u_\phi \omega _{\!R} - u_{\!R}\omega _\phi = 0. \end{equation}
Then, based on the definition of the mean flow vorticity,
$\boldsymbol{\omega } \equiv {\boldsymbol{\nabla }}\times \boldsymbol{u}$
, (A7) may be transformed algebraically into
\begin{equation} u_{\!R} \frac {\partial }{\partial R}(R \sin {\phi } \mathop {u_{\theta }}) + \frac {u_\phi }{R} \frac {\partial }{\partial \phi }(R \sin {\phi } \, u_{\theta }) = \boldsymbol{u}\boldsymbol{\cdot }{\boldsymbol{\nabla }}(R\sin {\phi }\,u_{\theta }) = 0. \end{equation}
Furthermore, in mirroring the treatment of (A3), we may examine the specific angular momentum,
$B$
, which is given by
\begin{equation} \mathinner {B(\psi )}= R \sin \phi \,u_{\theta } = \varGamma _{C}/(2\pi ), \end{equation}
where
$\varGamma _{C}$
denotes the circulation function. Equation (A9) can be interpreted as the dimensionless circulation per unit radian along a circular path that is centred around the
$z$
axis. Extending the rationale used to confirm that
$\mathinner {H(\psi )}$
satisfies (A3), it can be similarly shown that
$\mathinner {B(\psi )}$
inherently satisfies (A8).
At this juncture, the momentum equation can be simplified by leveraging the adiabatic relations between
$p$
and
$\rho$
. In view of these, the pressure-gradient term in the momentum equation can be reformulated using the identity
\begin{equation} \frac {{\boldsymbol{\nabla }}\! p}{\rho } = \frac {\gamma }{\gamma -1}\,{\boldsymbol{\nabla }}\!\left ( \frac {p}{\rho } \right )\!. \end{equation}
From (A10), the momentum equation becomes
\begin{equation} \frac {1}{2}\mathinner {{\boldsymbol{\nabla }}(\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{u})} + \frac {1}{M_{0}^2(\gamma -1)}\,{\boldsymbol{\nabla }}\!\left ( \frac {p}{\rho } \right ) = \boldsymbol{u}\times \boldsymbol{\omega }. \end{equation}
Additionally, by recognising that the first two terms in (A11) correspond to the gradient of
$\mathinner {H(\psi )}$
, this expression may be further reduced to
\begin{equation} \frac {{\rm d} H}{{\rm d} \psi }{\boldsymbol{\nabla }}\psi =\boldsymbol{u}\times \boldsymbol{\omega }. \end{equation}
In practice, the compressible BHE emerges from either the radial or polar components of (A12). One gets
\begin{gather} \textrm{L}^2 \psi + \rho ^2 \biggl ( B\frac {{\rm d} B}{{\rm d} \psi } - R^2 \sin ^2{\!\phi }\,\,\frac {{\rm d} H}{{\rm d} \psi } \biggr ) = \frac {{\boldsymbol{\nabla }}\psi \boldsymbol{\cdot }{\boldsymbol{\nabla }}\! \rho }{\rho } = {\boldsymbol{\nabla }}\psi \boldsymbol{\cdot }{\boldsymbol{\nabla }}(\,\ln \rho ), \end{gather}
where
${\boldsymbol{\nabla }}\! \rho /\rho \equiv {\boldsymbol{\nabla }}(\ln {\rho })$
is substituted, and where the Stokes operator,
$\textrm{L}^2\!$
, corresponds to
\begin{equation} \textrm{L}^2{} \equiv \frac {\partial ^2}{\partial R^2} + \frac {\sin {\phi }}{R^2}\frac {\partial }{\partial \phi }\!\left ( \frac {1}{\sin {\phi }}\frac {\partial }{\partial \phi } \right )\!. \end{equation}
Appendix B. Expansions of spherical Bessel functions
In what follows, we specify the spherical Bessel functions that appear in the streamfunction formulation associated with (3.26) for a hemispherically bounded cyclonic motion. In particular, we focus on the trigonometric representations of
$\mathinner {{j}_n{ ( z )}}$
and
$\mathinner {{y}_n{ ( z )}}$
:
\begin{align} \mathinner {\textrm{j}_n}\!{\left ( z \right )} &= \!\left ( \frac {\pi }{2z} \right )^{1/2}\mathinner {J_{n+\tfrac {1}{2}}{\!\left ( z \right )}} \nonumber \\&= \sin \!\left ( z - \tfrac {1}{2}n\pi \right )\sum _{k=0}^{{\!\left \lfloor n/2 \right \rfloor }}(-1)^k\tfrac {a_{2k}\bigl ( n+\tfrac {1}{2} \bigr )}{z^{2k+1}} \nonumber \\&\quad + \cos \!\left ( z - \tfrac {1}{2}n\pi \right )\sum _{k=0}^{{\!\left \lfloor (n-1)/2 \right \rfloor }}\!\!\!\!(-1)^k\tfrac {a_{2k+1}\bigl ( n+\tfrac {1}{2} \bigr )}{z^{2k+2}}, \end{align}
\begin{align} \mathinner {\textrm{y}_n}\!{\left ( z \right )} &= \!\left ( \frac {\pi }{2z} \right )^{1/2}\mathinner {Y_{n+\tfrac {1}{2}}{\!\left ( z \right )}} \nonumber \\&= -\cos \!\left ( z - \tfrac {1}{2}n\pi \right )\sum _{k=0}^{{\!\left \lfloor n/2 \right \rfloor }}(-1)^k\tfrac {a_{2k}\bigl ( n+\tfrac {1}{2} \bigr )}{z^{2k+1}} \nonumber \\&\quad + \sin \!\left ( z - \tfrac {1}{2}n\pi \right )\sum _{k=0}^{{\!\left \lfloor (n-1)/2 \right \rfloor }}\!\!\!\!(-1)^k\tfrac {a_{2k+1}\bigl ( n+\tfrac {1}{2} \bigr )}{z^{2k+2}}, \end{align}
where the floor function
$\!\left \lfloor x \right \rfloor$
returns the integer bracketed within
$x - 1 \lt {\!\left \lfloor x \right \rfloor } \leqslant x$
. As for the
$a_{k}$
coefficients, they are given by
\begin{equation} a_{k}\!\left ( n+\tfrac {1}{2} \right ) = \begin{cases} \dfrac {(n+k)!}{2^k k!\,(n-k)!}, & k = 0,1,\ldots ,n,\\ 0, & k = n+1,n+2,\ldots \, . \end{cases} \end{equation}
Appendix C. Closed-form expressions of first-order integrals
The functions
$F_1$
and
$F_2$
, which represent solutions to the system of ODEs given by (3.42), can be expressed using the general solution in (3.48). These describe the radial variations of the first-order streamfunction,
$\varPsi _1$
, within the chosen functional decomposition. From (3.48),
$F_k$
may be defined as
\begin{equation} \mathinner {F_k(x)} = x\mathinner {\textrm{j}_{2k-1}(x)}\mathopen [C_{2k-1} - I_{1,k}(x)\mathclose ] + x\mathinner {\textrm{y}_{2k-1}(x)}\mathopen [C_{2k} + I_{2,k}(x)\mathclose ], \end{equation}
where
$I_{1,k}$
and
$I_{2,k}$
denote antiderivatives that propagate leading-order effects through the inhomogeneous forcing functions,
$g_k$
, appearing on the right-hand side of (3.42). These antiderivatives may be written as
\begin{gather} I_{1,k}(x) = \int \frac {\textrm{y}_{2k-1}(x)g_k(x)}{x} {\rm d}{x} \end{gather}
and
\begin{gather} I_{2,k}(x) = \int \frac {\textrm{j}_{2k-1}(x)g_k(x)}{x} {\rm d}{x}, \end{gather}
where the constants of integration are excluded for simplicity. Within (C2),
$k = \mathopen \{ 1,2 \mathclose \}$
, and both
$\textrm{j}_{n}$
and
$\textrm{y}_{n}$
represent spherical Bessel functions of the first and second kinds. The two forcing terms,
$g_k$
, are prescribed by
\begin{align} \mathinner {g_2(x)} &= x \Bigl \{ x^2\mathinner {\textrm{j}_{1}^3(x)} - x \mathinner {\textrm{j}_{1}^2(x)}\mathinner {\textrm{j}_{2}(x)} - 10\mathinner {\textrm{j}_{1}(x)}\mathinner {\textrm{j}_{2}^2(x)} + 4x\mathinner {\textrm{j}_{2}^3(x)} \nonumber\\&\quad - \bigl [ x\mathinner {\textrm{j}_{2}(x)} - 2\mathinner {\textrm{j}_{1}(x)} \bigr ]^2 \mathinner {\textrm{j}_{3}(x)} \Bigr \}, \end{align}
\begin{align} \mathinner {g_1(x)} &= -8F_2(x) - 2\Bigl [ \tilde {C}_B^{\scriptscriptstyle {(1)}} x^3 - x\mathinner {\textrm{j}_{1}^2(x)} - 8\mathinner {\textrm{j}_{1}(x)}\mathinner {\textrm{j}_{2}(x)} + 3x\mathinner {\textrm{j}_{2}^2(x)} \Bigr ]\mathinner {\textrm{j}_{1}(x)}, \end{align}
where
$F_2$
denotes the second radial component of the ansatz in (3.41).
To compute
$g_1$
, it is necessary to first determine
$F_2$
by evaluating (C2) with
$k = 2$
, especially that the functions
$I_{1,2}$
and
$I_{2,1}$
remain central to the definition of
$F_2$
. Solving for
$F_2$
provides the basis for obtaining the solution for
$k = 1$
. For
$k = 2$
, the general solution to (C2) can be represented as a sum of integrals of the form
\begin{equation} I^{(abc\textit{d})}(x) = \int x^a \bigl ( \,\,\textrm{j}_{1}^b \boldsymbol{\cdot }\textrm{j}_{2}^c\boldsymbol{\cdot }\textrm{j}_{3}^d\boldsymbol{\cdot }f_3 \bigr )(x) {\rm d}{x}, \end{equation}
where
$f_3$
corresponds to
$\textrm{j}_{3}$
in (C2a
) and
$\textrm{y}_{3}$
in (C2b
). Thus, for
$k = 2$
, the integrals in (C2) return
\begin{equation} I_{n,2}(x) = I^{(2300)} - I^{(1210)} - 10 I^{(0120)} + 4 I^{(1030)} - I^{(2021)} + 4 I^{(1111)} - 4 I^{(0201)}. \end{equation}
Upon evaluating (C5) and grouping expressions by trigonometric and special functions, explicit relations for the integrals may be conveniently obtained. One gets
\begin{align} I_{1,2}(x) ={}& - \frac {75}{32 x^8} - \frac {47}{8 x^6} + \frac {13}{8 x^4} + \frac {1}{x^2} \end{align}
\begin{align} & + \biggl ( \frac {45}{44 x^{11}} - \frac {277}{44 x^9} + \frac {103}{22 x^7} + \frac {3779}{660 x^5} - \frac {3449}{1980 x^3} - \frac {8}{495 x} \biggr ) \sin (2 x) \end{align}
\begin{align} &- \biggl ( \frac {45}{88 x^{11}} - \frac {547}{88 x^9} + \frac {953}{88 x^7} - \frac {1651}{1320 x^5} + \frac {1}{990 x^3} - \frac {4}{495 x} \biggr ) \sin (4 x) \end{align}
\begin{align} &- \biggl ( \frac {45}{22 x^{10}} - \frac {217}{22 x^8} - \frac {101}{33 x^6} + \frac {5461}{1320 x^4} - \frac {511}{1980 x^2} \biggr ) \cos (2 x) \end{align}
\begin{align} & - \biggl ( \frac {45}{22 x^{10}} + \frac {4087}{352 x^8} - \frac {1361}{264 x^6} + \frac {41}{330 x^4} + \frac {1}{495 x^2} \biggr ) \cos (4 x) \end{align}
\begin{align} & + \frac {16}{495}\mathopen [\textrm{Ci}(2 x) - \textrm{Ci}(4 x)\mathclose ] \end{align}
and
\begin{align} I_{2,2}(x) ={}& - \frac {135}{88 x^{11}}+ \frac {51}{8 x^9} + \frac {27}{4 x^7} - \frac {179}{40 x^5} - \frac {2}{3 x^3} + \frac {1}{2 x} \end{align}
\begin{align} &+ \left (\frac {45}{22 x^{11}}-\frac {277}{22 x^9}-\frac {1}{88 x^7} + \frac {2773}{660 x^5} - \frac {3433}{1980 x^3} - \frac {16}{495 x}\right ) \cos (2 x) \end{align}
\begin{align} &- \left (\frac {45}{88 x^{11}}-\frac {547}{88 x^9}+\frac {953}{88 x^7}-\frac {1651}{1320 x^5}+\frac {1}{990 x^3} - \frac {4}{495 x}\right ) \cos (4 x) \end{align}
\begin{align} &+ \left (\frac {45}{11 x^{10}}-\frac {2647}{176 x^8}-\frac {989}{264 x^6}+\frac {6137}{1320 x^4} - \frac {527}{1980 x^2}\right ) \sin (2 x) \end{align}
\begin{align} &- \left (\frac {45}{22 x^{10}}-\frac {4087}{352 x^8}+\frac {1361}{264 x^6}-\frac {41}{330 x^4} - \frac {1}{495 x^2}\right ) \sin (4 x) \end{align}
\begin{align} &+ \frac {16}{495} \mathopen [\textrm{Si}(4 x) - 2 \textrm{Si}(2 x)\mathclose ], \end{align}
where
$\textrm{Si}(x)$
and
$\textrm{Ci}(x)$
stand for the sine and cosine integral functions, respectively:
\begin{equation} \textrm {Si}(x) = \int \frac {\sin {x}}{x}{\rm d}{x} \quad \text{and} \quad \textrm {Ci}(x) = \int \frac {\cos {x}}{x}{\rm d}{x}. \end{equation}
With
$F_2$
determined, the component integrals
$I_{1,1}$
and
$I_{2,1}$
, which prescribe
$F_1$
, become well-posed. They can be evaluated rather straightforwardly and expressed as
\begin{align} & \mathinner {I_{1,1}(x)} = \frac {7}{11x^6} + \frac {161}{220 x^4} + \frac {10 C_3 - 1}{x^2}\nonumber \\& + {\left ( \frac {2}{11 x^9} - \frac {29}{33 x^7} + \frac {5940 C_3-463}{495 x^5} + \frac {407-5940 C_3}{495 x^3} + \frac {495 \tilde {C}_B^{(1)}-64}{495 x} \right )}\sin (2x)\nonumber \\& + {\left ( \frac {4}{11 x^8} - \frac {16}{11 x^6} + \frac {23760 C_3-3007}{990 x^4} + \frac {64-1980 C_3}{990 x^2} + \frac {\tilde {C}_B^{(1)}}{2} \right )} \cos (2x)\nonumber \\& + {\left ( \frac {1}{11 x^9} - \frac {65}{66 x^7} + \frac {79}{90 x^5} - \frac {32}{495 x^3} - \frac {32}{495 x} \right )}\sin (4 x)\nonumber \\& - {\left ( \frac {4}{11 x^8} - \frac {17}{11 x^6} + \frac {91}{1980 x^4} + \frac {8}{495 x^2} \right )}\cos (4x)\nonumber \\& + {\left [ \frac {32}{495}{\left ( \frac {5}{x^2}+4 \right )} + \frac {32}{495}{\left ( \frac {12}{x^4}-\frac {1}{x^2} \right )}\cos (2x) - \frac {64}{165}{\left ( \frac {1}{x^5}-\frac {1}{x^3} \right )}\sin (2x) \right ]}\nonumber \\& \times {\left [ \textrm {Ci}(4x) - \textrm {Ci}(2x) \right ]}\nonumber \\& - {\left [ \frac {64}{495}{\left ( \frac {3}{x^5}+\frac {3}{x^3}-\frac {1}{x} \right )} + \frac {64}{165}{\left ( \frac {1}{x^5}-\frac {1}{x^3} \right )}\cos (2x) + \frac {32}{495}{\left ( \frac {12}{x^4}-\frac {1}{x^2} \right )}\sin (2x) \right ]}\nonumber \\& \times {\left [ \textrm {Si}(4x) - 2\textrm {Si}(2x) \right ]} \end{align}
and
\begin{align} & \mathinner {I_{2,1}}\!(x) = \frac {3}{11 x^9} - \frac {17}{22 x^7}+\frac {3960 C_3-919}{330 x^5}+\frac {3960 C_3-534}{330 x^3} + \frac {330 \tilde {C}_B^{(1)}\!-1320 C_3}{330 x}-\tilde {C}_B^{(1)} x\nonumber \\& - {\left [ \frac {8}{11 x^8} - \frac {18}{11 x^6} + \frac {23760 C_3-3557}{990 x^4} - \frac {2 (495 C_3+16)}{495 x^2} + \frac {\tilde {C}_B^{(1)}}{2} \right ]}\sin (2x)\nonumber \\& - {\left [ \frac {4}{11 x^9} - \frac {58}{33 x^7} + \frac {4 (1485 C_3-326)}{495 x^5} + \frac {391-5940 C_3}{495 x^3} + \frac {495 \tilde {C}_B^{(1)}-128}{495 x} \right ]}\cos (2x)\nonumber \\& +{\left ( \frac {4}{11 x^8}-\frac {17}{11 x^6}+\frac {91}{1980 x^4}+\frac {8}{495 x^2} \right )}\sin (4x) \nonumber \\& + {\left ( \frac {1}{11 x^9}-\frac {65}{66 x^7}+\frac {79}{90 x^5}-\frac {32}{495 x^3}-\frac {32}{495 x} \right )} \cos (4x) \nonumber \\& + {\left [ \frac {32}{99 x^2}+\frac {128}{495} - {\left ( \frac {128}{165 x^4}-\frac {32}{495 x^2} \right )}\cos (2x) + {\left ( \frac {64}{165 x^5} - \frac {64}{165 x^3} \right )}\sin (2x)\! \right ]}\nonumber \\& \times \left [ 2\mathinner {\textrm {Si}(2x)} - \mathinner {\textrm {Si}(4x)} \right ] + \left [ 2\mathinner {\textrm {Ci}(2x)} - \mathinner {\textrm {Ci}(4x)} \right ] \nonumber \\& \times {\left [ \frac {64}{495}{\left ( \frac {3}{x^5}+\frac {3}{x^3}-\frac {1}{x} \right )} - \frac {64}{165}{\left ( \frac {1}{x^5} - \frac {1}{x^3} \right )}\cos (2x) - \frac {32}{495}{\left ( \frac {12}{x^4} - \frac {1}{x^2} \right )}\sin (2x) \right ]}. \end{align}
This concludes our analytical evaluation of the first-order integrals.
Appendix D. Asymptotic expansion for an implicitly defined surface
In what follows, we describe the asymptotic procedure used to determine the spatial coordinates of a surface that is defined implicitly by the extremum of a scalar function. This approach is undertaken to obtain two-term asymptotic expressions for the mantle radius,
$r_{{m}}(z)$
, as well as the vortex core radius,
$r_{{c}}(z)$
. To begin, let the surface be specified by the condition
\begin{equation} \frac {\partial f(r, z; \epsilon )}{\partial r} = 0, \end{equation}
where
$f$
and
$\epsilon \ll 1$
denote a scalar function and a small perturbation parameter, respectively. In this work, we take
$\epsilon = M_0^2$
and let the surface defining position,
$r(z; \epsilon )$
, be a function of both
$z$
and
$\epsilon$
. To proceed, asymptotic expansions for
$f$
and
$r$
may be sequentially written in powers of
$\epsilon$
using
\begin{align} f(r, z; \epsilon ) & = f_0(r, z) + \epsilon f_1(r, z) + {O}(\epsilon ^2), \end{align}
\begin{align} r(z; \epsilon ) & = r_0(z) + \epsilon r_1(z) + {O}(\epsilon ^2). \end{align}
In order to obtain
$r_0(z)$
and
$r_1(z)$
, (D1) at
$r(z;\epsilon )$
can be decomposed using a Taylor series expansion around the leading-order radius,
$r_0(z)$
. One simply puts
\begin{equation} \!\left . \frac {\partial f}{\partial r} \right |_{r = r_0 +{\kern1pt} \epsilon r_1} = \left . \frac {\partial f}{\partial r} \right |_{r=r_0} + \epsilon r_1(z) \!\left . \frac {\partial ^2 f}{\partial r^2} \right |_{r=r_0} + {O}(\epsilon ^2) = 0. \end{equation}
This condition can be interlaced with the asymptotic expansion of
$f$
in (D2), namely,
\begin{equation} \!\left . \frac {\partial }{\partial r} (f_0 + \epsilon f_1) \right |_{r=r_0} + \epsilon r_1(z) \!\left . \frac {\partial ^2}{\partial r^2} (f_0 + \epsilon f_1) \right |_{r=r_0} + {O}(\epsilon ^2) = 0. \end{equation}
By distributing the derivatives and collecting terms in equal powers of
$\epsilon$
, we obtain
\begin{equation} \!\left ( \!\left . \frac {\partial f_0}{\partial r} \right |_{r=r_0} \right ) + \epsilon \!\left [ \!\left . \frac {\partial f_1}{\partial r} \right |_{r=r_0} + r_1(z) \!\left . \frac {\partial ^2 f_0}{\partial r^2} \right |_{r=r_0} \right ] + {O}(\epsilon ^2) \equiv 0;\quad \forall \epsilon \ll 1. \end{equation}
For this parametric equation to hold, the coefficients of each power of
$\epsilon$
must vanish independently. This enables us to identify the leading and first-order equations, specifically
\begin{align} \!\left . \frac {\partial f_0}{\partial r} \right |_{r=r_0(z)} & = 0, \end{align}
\begin{align} \!\left . \frac {\partial f_1}{\partial r} \right |_{r=r_0(z)} + r_1(z) \!\left . \frac {\partial ^2 f_0}{\partial r^2} \right |_{r=r_0(z)} & = 0. \end{align}
After specifying the leading-order approximation,
$r_0(z)$
, the first-order correction,
$r_1(z)$
, can be readily extracted provided that the extremum remains non-degenerate. We get
\begin{equation} r_1(z) = - \!\left . \frac {\partial f_1 / \partial r }{\partial ^2 f_0 / \partial r^2}\right |_{r = r_0(z)}; \quad \text{} \quad \!\left . \frac {\partial ^2 f_0}{\partial r^2} \right |_{r=r_0(z)} \neq 0. \end{equation}
In view of (D9), and setting
$f=\psi$
, one may simply substitute
$r_0=r_{{m}}^{\scriptscriptstyle {(0)}}$
and
$r_1=r_{{m}}^{\scriptscriptstyle {(1)}}$
in § 4.2.2 or, alternatively,
$r_0=r_{{c}}^{\scriptscriptstyle {(0)}}$
and
$r_1=r_{{c}}^{\scriptscriptstyle {(1)}}$
in § 4.3.4.
Appendix E. Reduction of helicity integrals
Following Moffatt (Reference Moffatt1969), the dimensionless helicity
$\mathcal{H}$
may be written as
\begin{equation} \mathcal{H} \equiv \int _{\mathcal{V}} h \,{\rm d}{V} = \int _{\mathcal{V}} \boldsymbol{\omega }\boldsymbol{\cdot }\boldsymbol{u}\,{\rm d}{V}=-C_B\int _{\mathcal{V}} \rho \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{u}\, {\rm d}{V}, \end{equation}
where the Beltramian relation between
$\boldsymbol{\omega }$
and
$\boldsymbol{u}$
is used. Next, the integrand may be expressed in terms of the streamfunction by setting
\begin{equation} \rho \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u} = \frac {{\boldsymbol{\nabla }}\psi \boldsymbol{\cdot }{\boldsymbol{\nabla }}\psi }{\rho r^2} + C_B^2\frac {\rho \psi ^2}{r^2}, \end{equation}
where
$r = R\sin {\phi }$
. The first term in the above may be further transformed by means of the vector identity
\begin{equation} \frac {{\boldsymbol{\nabla }}\psi \boldsymbol{\cdot }{\boldsymbol{\nabla }}\psi }{\rho r^2} = {\boldsymbol{\nabla }}\boldsymbol{\cdot }{\left ( \psi \frac {{\boldsymbol{\nabla }}\psi }{\rho r^2} \right )} - \psi {\boldsymbol{\nabla }}\boldsymbol{\cdot }{\left ( \frac {{\boldsymbol{\nabla }}\psi }{\rho r^2} \right )}. \end{equation}
Coincidentally, the last term in (E3) may be eliminated using an algebraically rearranged form of the compressible BHE, specifically
\begin{equation} {\boldsymbol{\nabla }}\boldsymbol{\cdot }{\left ( \frac {{\boldsymbol{\nabla }}\psi }{\rho r^2} \right )} =-C_B^2\frac {\rho \psi }{r^2} . \end{equation}
Then, through backward substitution into (E2) and (E1), one arrives at a markedly compact form of the total helicity, namely,
\begin{equation} \mathcal{H} = -C_B\int _{\mathcal{{\partial } V}}\frac {\psi {\boldsymbol{\nabla }}\psi }{\rho r^2} \boldsymbol{\cdot }\hat {\boldsymbol{n}}\,{\rm d}{S} - 2C_B^3 \int _{\mathcal{V}}\frac {\rho \psi ^2}{r^2}{\rm d}{V}, \end{equation}
where the volume integral of the divergence has been converted into a surface integral using Gauss’s theorem. Moreover, it may be shown that the integrand of this first term vanishes identically at all points taken along the boundaries. As for the remaining volume integral, it may be further reduced using the limits of integration to obtain
\begin{equation} \mathcal{H} = - 4\pi C_B^3 \int _0^1\int _0^{\pi /2}\frac {\rho \psi ^2}{\sin {\phi }}{\rm d}{\phi }\,{\rm d}{R}. \end{equation}
Lastly, by taking the RJ expansion of (E6), setting
$\mathcal{H} = \mathcal{H}_0 + M_0^2\mathcal{H}_1 + {O}(M_0^4)$
and evaluating the integrals at each successive order, we recover
\begin{equation} \mathcal{H}_0 = -4\pi {\left ( C_B^{\scriptscriptstyle {(0)}} \right )}^3\int _0^1\int _0^{\pi /2}\frac {\psi _0^2}{\sin \phi }{\rm d}{\phi }\,{\rm d}{R}=-320.4\kappa ^2 \end{equation}
and
\begin{equation} \mathcal{H}_1 = -4\pi {\left ( C_B^{\scriptscriptstyle {(0)}} \right )}^3\int _0^1\int _0^{\pi /2}\frac {\psi _0}{\sin {\phi }}\biggl [ \biggl ( 3\frac {C_B^{\scriptscriptstyle {(1)}}}{C_B^{\scriptscriptstyle {(0)}}} + \rho _1 \biggr )\psi _0 + 2\psi _1 \biggr ]{\rm d}{\phi }\,{\rm d}{R}=-19041\kappa ^4. \end{equation}
We hence arrive at the explicit dependence of the leading- and first-order helicity components on the inflow parameter,
$\kappa$
.
Open access
