2.1 Normalization

Before deriving the compressible relations, the governing equations may be conveniently converted to a dimensionless form. Using overbars to designate dimensional quantities, we take

(2.1a-d ) $$\begin{eqnarray}r=\frac{\bar{r}}{a};\quad z=\frac{\bar{z}}{a};\quad \unicode[STIX]{x1D735}=a\bar{\unicode[STIX]{x1D735}};\quad \unicode[STIX]{x1D6FD}=\frac{b}{a};\end{eqnarray}$$

(2.2a-f ) $$\begin{eqnarray}u=\frac{\bar{u}}{U};\quad v=\frac{\bar{v}}{U};\quad w=\frac{\bar{w}}{U};\quad \unicode[STIX]{x1D734}=\frac{a\bar{\unicode[STIX]{x1D734}}}{U};\quad \unicode[STIX]{x1D713}=\frac{\bar{\unicode[STIX]{x1D713}}}{\unicode[STIX]{x1D70C}_{0}Ua^{2}};\quad H=\frac{\bar{H}}{U^{2}};\end{eqnarray}$$

(2.3a-f ) $$\begin{eqnarray}p=\frac{\bar{p}}{p_{0}};\quad \unicode[STIX]{x1D70C}=\frac{\bar{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x1D70C}_{0}};\quad Q_{i}=\frac{\bar{Q}_{i}}{Ua^{2}}=\frac{A_{i}}{a^{2}};\quad Q_{o}=\frac{\bar{Q}_{o}}{Ua^{2}};\quad {\dot{m}}_{i}=\frac{\dot{\bar{m}}_{i}}{\unicode[STIX]{x1D70C}_{0}Ua^{2}};\quad {\dot{m}}_{o}=\frac{\dot{\bar{m}}_{o}}{\unicode[STIX]{x1D70C}_{0}Ua^{2}},\end{eqnarray}$$

where $(r,z)$ represent the two primary spatial coordinates and $b$ , the chamber exit radius. In (2.2), $(u,v,w)$ denote the radial, tangential and axial components of the velocity, whereas $\unicode[STIX]{x1D734}$ , $\unicode[STIX]{x1D713}$ and $H$ stand for the vorticity, streamfunction and stagnation enthalpy, respectively. As usual, $p$ and $\unicode[STIX]{x1D70C}$ refer to the pressure and density, while $Q$ and ${\dot{m}}$ refer to the volumetric and mass flow rates, respectively. In (2.1), all spatial parameters are normalized by the chamber radius, $a$ . Similarly, the wall-tangential injection velocity, $U$ , is used to non-dimensionalize the aforementioned velocities and other related variables as needed. Based on this choice of normalized quantities, the equations of motion become

(2.4) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{u})=0\quad \text{(conservation of mass);} & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}H\boldsymbol{u})=0\quad \text{(conservation of energy);} & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-{\displaystyle \frac{\unicode[STIX]{x1D735}p}{\unicode[STIX]{x1D6FE}M_{0}^{2}\unicode[STIX]{x1D70C}}}\quad \text{(conservation of momentum);} & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & H={\displaystyle \frac{1}{2}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}+{\displaystyle \frac{1}{M_{0}^{2}(\unicode[STIX]{x1D6FE}-1)}}{\displaystyle \frac{p}{\unicode[STIX]{x1D70C}}}\quad \text{(stagnation enthalpy)}, & \displaystyle\end{eqnarray}$$

where $M_{0}=U/c_{0}$ , with $c_{0}$ representing the reference speed of sound in the chamber. For convenience, we define the axisymmetric Stokes operator $\text{D}^{2}$ and the compressible streamfunction–velocity relations as

(2.8a-c ) $$\begin{eqnarray}\text{D}^{2}\equiv \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}-\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}};\quad u=-\frac{1}{\unicode[STIX]{x1D70C}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z};\quad w=\frac{1}{\unicode[STIX]{x1D70C}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r}.\end{eqnarray}$$

2.2 Density–streamfunction formulation

For the compressible Bragg–Hawthorne framework, it is useful to recast the stagnation enthalpy and tangential angular momentum in terms of $\unicode[STIX]{x1D713}$ . By substituting the conservation of mass equation into the energy conservation equation, one may write

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D735}H\boldsymbol{\cdot }\boldsymbol{u}=0.\end{eqnarray}$$

Equation (2.9) suggests that $H$ can only vary along directions that are orthogonal to the velocity field. Hence, $H=H(\unicode[STIX]{x1D713})$ remains constant along streamlines.

With the stagnation enthalpy in hand, the expanded tangential component of the momentum conservation expression in (2.6) may be revisited. Given the underlying assumption of axisymmetry, the $\unicode[STIX]{x1D703}$ -momentum equation may be reduced to

(2.10) $$\begin{eqnarray}u\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}r}+w\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}z}+\frac{uv}{r}=0.\end{eqnarray}$$

Equation (2.10) can be further multiplied by $r$ to the extent of reproducing the material derivative of $B\equiv rv$ , where B represents the tangential angular momentum

(2.11) $$\begin{eqnarray}ru\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}r}+rw\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}z}+uv=\frac{\text{d}(rv)}{\text{d}t}=\frac{\text{d}B}{\text{d}t}=0.\end{eqnarray}$$

A vanishing material derivative confirms that $B=B(\unicode[STIX]{x1D713})$ must remain invariant along streamlines. Having $H$ and $B$ in streamfunction form proves helpful in simplifying the momentum equation. In fact, this step is facilitated when paired with the isentropic flow relation, namely, $p=K\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}}$ , where $K$ denotes a general constant.

At this juncture, one may recognize that the momentum relation given by (2.6) incorporates a pressure gradient divided by the density, and so an equivalent form may be realized by manipulating the isentropic relation to recover

(2.12) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D6FE}-1}\unicode[STIX]{x1D735}\left(\frac{p}{\unicode[STIX]{x1D70C}}\right)=\frac{\unicode[STIX]{x1D735}p}{\unicode[STIX]{x1D70C}}.\end{eqnarray}$$

This expression may be inserted on the right-hand side of (2.6). Then using the vector identity, $\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=(1/2)\unicode[STIX]{x1D735}\boldsymbol{u}^{2}-\boldsymbol{u}\times \unicode[STIX]{x1D735}\times \boldsymbol{u}$ , we extract

(2.13) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D735}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u})}{2}+\frac{1}{M_{0}^{2}(\unicode[STIX]{x1D6FE}-1)}\unicode[STIX]{x1D735}\left(\frac{p}{\unicode[STIX]{x1D70C}}\right)-\boldsymbol{u}\times \unicode[STIX]{x1D734}=0.\end{eqnarray}$$

By inspection of (2.7), it may be recognized that the first two terms in the above correspond to the spatial gradient of the total enthalpy. We hence arrive at

(2.14) $$\begin{eqnarray}\unicode[STIX]{x1D735}H=\boldsymbol{u}\times \unicode[STIX]{x1D734}.\end{eqnarray}$$

As with the vorticity–streamfunction, the axial component of (2.14) may be segregated after imposing the axisymmetry condition. To eliminate the vorticity, the right-hand side of (2.14) may be expanded in terms of the velocity, viz.

(2.15) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}z}=-u\left(\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}r}-\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}\right)+v\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}z}.\end{eqnarray}$$

Next, the velocities may be eliminated in favour of the streamfunction and the angular momentum is factored out to obtain

(2.16) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}z}=\frac{1}{\unicode[STIX]{x1D70C}^{2}r^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z^{2}}-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r^{2}}-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}r}-\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r}\right)+\frac{B}{r^{2}}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}z}.\end{eqnarray}$$

Recognizing that $\unicode[STIX]{x2202}H/\unicode[STIX]{x2202}z=(\text{d}H/\text{d}\unicode[STIX]{x1D713})(\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}z)$ and $\unicode[STIX]{x2202}B/\unicode[STIX]{x2202}z=(\text{d}B/\text{d}\unicode[STIX]{x1D713})(\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}z)$ and simplifying leads to

(2.17) $$\begin{eqnarray}r^{2}\frac{\text{d}H}{\text{d}\unicode[STIX]{x1D713}}-B\frac{\text{d}B}{\text{d}\unicode[STIX]{x1D713}}=\frac{1}{\unicode[STIX]{x1D70C}^{2}}\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z^{2}}-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r^{2}}-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}r}-\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r}\right).\end{eqnarray}$$

At length, taking advantage of vector notations and using the $\text{D}^{2}$ operator defined in (2.8), we arrive at the compact and convenient form

(2.18) $$\begin{eqnarray}\text{D}^{2}\unicode[STIX]{x1D713}+\unicode[STIX]{x1D70C}^{2}\left(B\frac{\text{d}B}{\text{d}\unicode[STIX]{x1D713}}-r^{2}\frac{\text{d}H}{\text{d}\unicode[STIX]{x1D713}}\right)=\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}.\end{eqnarray}$$

In the above, conservation of mass (2.4) and energy (2.5) have been combined into (2.9), namely, to show that the stagnation enthalpy, $H$ , must be necessarily a function of $\unicode[STIX]{x1D713}$ . In a similar manner, the tangential component of (2.6) may be considered in the context of axisymmetric motion to ascertain the exclusive dependence of the tangential angular momentum, $B$ , on the streamfunction. In fact, the use of axisymmetry leads to the decoupling of the axial component of the momentum equation from its radial and tangential counterparts in (2.14). Furthermore, the isentropic relation (2.12) enables us to express the outcome in terms of the stagnation enthalpy, vorticity and velocity. These, in turn, can be written as a function of $\unicode[STIX]{x1D713}$ , namely, after introducing the latter and taking into account the properties of $B$ and $H$ , thus leading to the compressible Bragg–Hawthorne equation (2.18). Therefore, given the general dependence of these quantities on $\unicode[STIX]{x1D713}$ , physically suitable forms of $H$ and $B$ may be specified analogously to the manner by which velocity components and streamfunctions are often posited while seeking similarity solutions to the Navier–Stokes equations. It is this flexibility that sets the Bragg–Hawthorne technique apart, particularly as a versatile and promising framework that can help to unravel multiple solutions for the same geometry and physical model.

2.3 Compressible energy relation

To achieve closure in the compressible Bragg–Hawthorne equation, a density relation is required. Hence, in the spirit of complementing (2.18), the expression for total enthalpy may be used. By rewriting (2.7) in terms of $\unicode[STIX]{x1D713}$ and using the isentropic assumption to eliminate pressure in favour of density, we find

(2.19) $$\begin{eqnarray}H-\frac{B^{2}}{2r^{2}}=\frac{1}{2\unicode[STIX]{x1D70C}^{2}r^{2}}\left[\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r}\right)^{2}\right]+\frac{1}{M_{0}^{2}(\unicode[STIX]{x1D6FE}-1)}\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}-1}.\end{eqnarray}$$

Equations (2.18) and (2.19) provide the basis for the compressible framework that we plan to establish for the steady-state analysis of axisymmetric gaseous motions in a friction-free environment.

3.1 Rayleigh–Janzen expanded equations

Applying the Rayleigh–Janzen series expansions in (3.1) to the compressible Bragg–Hawthorne equation and segregating leading and first-order quantities, one gets

(3.2) $$\begin{eqnarray}O(1):\text{D}^{2}\unicode[STIX]{x1D713}_{0}+B_{0}{\displaystyle \frac{\text{d}B_{0}}{\text{d}\unicode[STIX]{x1D713}}}-r^{2}{\displaystyle \frac{\text{d}H_{0}}{\text{d}\unicode[STIX]{x1D713}}}=0;\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle O(M_{0}^{2}):\text{D}^{2}\unicode[STIX]{x1D713}_{1}+B_{1}{\displaystyle \frac{\text{d}B_{1}}{\text{d}\unicode[STIX]{x1D713}}}-r^{2}{\displaystyle \frac{\text{d}H_{1}}{\text{d}\unicode[STIX]{x1D713}}} & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x2202}z}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}}+{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x2202}r}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D70C}_{1}\left[\text{D}^{2}\unicode[STIX]{x1D713}_{0}+3\left(B_{0}{\displaystyle \frac{\text{d}B_{0}}{\text{d}\unicode[STIX]{x1D713}}}-r^{2}{\displaystyle \frac{\text{d}H_{0}}{\text{d}\unicode[STIX]{x1D713}}}\right)\right].\end{eqnarray}$$

Consistent with conventional perturbation theory, the leading-order equation reduces to the traditional incompressible Bragg–Hawthorne equation. The first-order correction, however, encapsulates the $O(M_{0}^{2})$ compressible contribution. At first order, its left-hand side mirrors the leading-order operator while the terms on the right-hand side give rise to a non-homogeneous partial differential equation (PDE).

The same procedure may be straightforwardly applied to the stagnation enthalpy in (2.19). As usual, by segregating terms of the same order, we recover

(3.4) $$\begin{eqnarray}O(1):H_{0}-\frac{B_{0}^{2}}{2r^{2}}=\frac{1}{2r^{2}}\left[\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}\right)^{2}\right]+\frac{\unicode[STIX]{x1D6FE}+1}{\unicode[STIX]{x1D6FE}-1}\unicode[STIX]{x1D70C}_{1};\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle O(M_{0}^{2}):2\unicode[STIX]{x1D70C}_{1}\left(H_{0}-\frac{B_{0}}{2r^{2}}\right)+H_{1}-\frac{B_{0}B_{1}}{2r^{2}} & = & \displaystyle \frac{1}{2r^{2}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}r}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x1D6FE}+1}{\unicode[STIX]{x1D6FE}-1}(\unicode[STIX]{x1D70C}_{2}+\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D70C}_{1}^{2}).\end{eqnarray}$$

When (3.2) is used to solve for $\unicode[STIX]{x1D713}_{0}$ , substitution into (3.4) directly unravels the density correction, $\unicode[STIX]{x1D70C}_{1}$ . With the density in hand, the right-hand side of (3.3) becomes fully determined and the resulting non-homogeneous PDE may be solved for the first compressible streamfunction correction. In principle, this sequence may be repeated until a satisfactory truncation error is reached. Presently, the procedure will enable us to extract closed-form expressions for the leading- and first-order corrections. Generally speaking, the complexity of the particular solutions can grow rapidly to the extent that a compressible approximation at the second order or beyond may require considerable effort. In most problems, however, the first-order compressible correction will be sufficiently accurate to convey the bulk compressibility effects, and this may be largely attributed to the typical values of $M_{0}^{2}$ . This provision is especially true for swirl-dominated flows such as those arising in the context of a bidirectional vortex engine in which the reference Mach number remains much smaller than unity.

3.2 Specification of $B$ and $H$

Modelling the bidirectional vortex, or any other motion for that matter, begins with the selection of suitable forms for $B$ and $H$ in the compressible Bragg–Hawthorne equations. To facilitate analytical closure, several test functions may be considered, specifically

(3.6a,b ) $$\begin{eqnarray}\displaystyle B\frac{\text{d}B}{\text{d}\unicode[STIX]{x1D713}}=\text{const.}\quad B=\sqrt{B_{0}\unicode[STIX]{x1D713}+B_{1}};\quad B\frac{\text{d}B}{\text{d}\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D713}\quad B=\sqrt{B_{0}\unicode[STIX]{x1D713}^{2}+B_{1}}; & & \displaystyle\end{eqnarray}$$

(3.7a,b ) $$\begin{eqnarray}\displaystyle \frac{\text{d}H}{\text{d}\unicode[STIX]{x1D713}}=\text{const.}\quad H=H_{0}\unicode[STIX]{x1D713}+H_{1};\quad \frac{\text{d}H}{\text{d}\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D713}\quad H=H_{0}\unicode[STIX]{x1D713}^{2}+H_{1}. & & \displaystyle\end{eqnarray}$$

Although the number of candidate functions may be limitless, the selections above lead to linear relations that increase the likelihood of producing analytical formulations. Higher-order polynomial relations may require a numerical treatment of the density–streamfunction equations. For example, the (original) incompressible model of the bidirectional vortex by Vyas & Majdalani (Reference Vyas and Majdalani2006) may be recovered by setting $B=1$ and $\text{d}H/\text{d}\unicode[STIX]{x1D713}=-C_{n}^{2}\unicode[STIX]{x1D713}$ , where $C_{n}$ is a constant. To make further headway in illustrating this procedure, one may attempt to follow Bloor & Ingham (Reference Bloor and Ingham1987) or Majdalani (Reference Majdalani2012) by specifying $B$ and $H$ such that

(3.8a-c ) $$\begin{eqnarray}\frac{\text{d}H}{\text{d}\unicode[STIX]{x1D713}}=0;\quad B=\sqrt{C_{0}^{2}\unicode[STIX]{x1D713}^{2}+C_{1}^{2}};\quad \frac{\text{d}B}{\text{d}\unicode[STIX]{x1D713}}=\frac{C_{0}^{2}\unicode[STIX]{x1D713}}{\sqrt{C_{0}^{2}\unicode[STIX]{x1D713}^{2}+C_{1}^{2}}}.\end{eqnarray}$$

In addition to being consistent with the axisymmetric system of equations described in the previous section, the foregoing expressions for $H$ and $B$ may be physically associated with a homenergic flow field that is capable of satisfying the problem’s boundary conditions while still leading to a linear Bragg–Hawthorne equation. By considering a fluid that originates from a pressurized feed system, which may be viewed as a uniform source of energy, a homenergic expression of $\text{d}H/\text{d}\unicode[STIX]{x1D713}=0$ may be used. The choice of $B$ may be similarly ascribed to the need to consider the simplest expression that remains compliant with the angular momentum constraints and corresponding boundary conditions, while still giving rise to a linear problem.

Interestingly, it turns out that, in the compressible case, these declarations prove insufficient in reproducing a congruent first-order system of equations. The source of this disparity may be traced back to the right-hand side of (3.3), where third-order multiples of the streamfunction emerge. As per (3.4), the density correction contains $\unicode[STIX]{x1D713}_{0}^{2}$ terms and these are multiplied by another $\unicode[STIX]{x1D713}_{0}$ during final book keeping. To compensate for these additional powers of $\unicode[STIX]{x1D713}$ , a modification of (3.8) is warranted. This may be accomplished by taking $\text{d}H/\text{d}\unicode[STIX]{x1D713}=0$ and writing

(3.9) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}B=\sqrt{C_{0}^{2}\unicode[STIX]{x1D713}^{2}+C_{1}^{2}+M_{0}^{2}\unicode[STIX]{x1D713}^{2}\left(C_{2}^{2}+{\textstyle \frac{1}{2}}C_{3}^{2}\unicode[STIX]{x1D713}^{2}\right)};\\ {\displaystyle \frac{\text{d}B}{\text{d}\unicode[STIX]{x1D713}}}={\displaystyle \frac{C_{0}^{2}\unicode[STIX]{x1D713}+M_{0}^{2}\unicode[STIX]{x1D713}(C_{2}^{2}+C_{3}^{2}\unicode[STIX]{x1D713}^{2})}{\sqrt{C_{0}^{2}\unicode[STIX]{x1D713}^{2}+C_{1}^{2}+M_{0}^{2}\unicode[STIX]{x1D713}^{2}\left(C_{2}^{2}+{\textstyle \frac{1}{2}}C_{3}^{2}\unicode[STIX]{x1D713}^{2}\right)}}}.\end{array}\right\}\end{eqnarray}$$

It may be instructive to note that the reference Mach number, $M_{0}$ , remains invariant under steady-state flow conditions. As such, its inclusion in the fundamental definition of $B$ does not violate in any way the streamfunction constraint. Furthermore, realizing that $B$ and $\text{d}B/\text{d}\unicode[STIX]{x1D713}$ appear only as a product in the Bragg–Hawthorne equation, their combination may be expanded as:

(3.10) $$\begin{eqnarray}B\frac{\text{d}B}{\text{d}\unicode[STIX]{x1D713}}=C_{0}^{2}\unicode[STIX]{x1D713}+M_{0}^{2}\unicode[STIX]{x1D713}(C_{2}^{2}+C_{3}^{2}\unicode[STIX]{x1D713}^{2}).\end{eqnarray}$$

From an asymptotic standpoint, equation (3.10) does not entail a loss of generality. It is obtained by expanding the angular momentum and its derivative to the appropriate truncation order before substituting the outcome into the streamfunction relation. The next step is to insert the perturbed form of $\unicode[STIX]{x1D713}$ and write:

(3.11) $$\begin{eqnarray}\displaystyle B\frac{\text{d}B}{\text{d}\unicode[STIX]{x1D713}} & = & \displaystyle C_{0}^{2}(\unicode[STIX]{x1D713}_{0}+M_{0}^{2}\unicode[STIX]{x1D713}_{1}+M_{0}^{4}\unicode[STIX]{x1D713}^{2})\nonumber\\ \displaystyle & & \displaystyle +\,M_{0}^{2}[C_{2}^{2}(\unicode[STIX]{x1D713}_{0}+M_{0}^{2}\unicode[STIX]{x1D713}_{1}+M_{0}^{4}\unicode[STIX]{x1D713}^{2})+C_{3}^{2}(\unicode[STIX]{x1D713}_{0}+M_{0}^{2}\unicode[STIX]{x1D713}_{1}+M_{0}^{4}\unicode[STIX]{x1D713}^{2})^{3}],\end{eqnarray}$$

where the constants $C_{0}$ , $C_{2}$ and $C_{3}$ may be determined from the fundamental boundary conditions, as shown in the following section. Additionally, $H$ may be specified in a manner to achieve a congruent expansion of the density relation. If we take $p_{0}$ and $\unicode[STIX]{x1D70C}_{0}$ to be stagnation properties, we can write $H=1/[(\unicode[STIX]{x1D6FE}-1)M_{0}^{2}]$ .

3.3 Compressible Bragg–Hawthorne solution strategy

By gathering $O(M_{0}^{2})$ in $B\,\text{d}B/\text{d}\unicode[STIX]{x1D713}$ , we retrieve

(3.12) $$\begin{eqnarray}B\frac{\text{d}B}{\text{d}\unicode[STIX]{x1D713}}=C_{0}^{2}\unicode[STIX]{x1D713}_{0}+M_{0}^{2}(C_{0}^{2}\unicode[STIX]{x1D713}_{1}+C_{2}^{2}\unicode[STIX]{x1D713}_{0}+C_{3}^{2}\unicode[STIX]{x1D713}_{0}^{3})+O(M_{0}^{4}).\end{eqnarray}$$

Inserting these contributions back into (3.2) and (3.3) gives rise to a congruent set of linearized Bragg–Hawthorne equations at the first two successive perturbation orders, namely,

(3.13) $$\begin{eqnarray}\displaystyle O(1):\text{D}^{2}\unicode[STIX]{x1D713}_{0}+C_{0}^{2}\unicode[STIX]{x1D713}_{0}=0; & & \displaystyle\end{eqnarray}$$

(3.14) $$\begin{eqnarray}\displaystyle O(M_{0}^{2}):\text{D}^{2}\unicode[STIX]{x1D713}_{1}+C_{0}\unicode[STIX]{x1D713}_{1} & = & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D70C}_{1}(\text{D}^{2}\unicode[STIX]{x1D713}_{0}+3C_{0}^{2}\unicode[STIX]{x1D713}_{0})-C_{2}^{2}\unicode[STIX]{x1D713}_{0}-C_{3}^{2}\unicode[STIX]{x1D713}_{0}^{3}.\end{eqnarray}$$

Equation (3.14) can be further simplified by realizing that the left-hand side of (3.13) partially appears on its right-hand side. This permits reducing (3.14) into

(3.15) $$\begin{eqnarray}\text{D}^{2}\unicode[STIX]{x1D713}_{1}+C_{0}\unicode[STIX]{x1D713}_{1}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}-2\unicode[STIX]{x1D70C}_{1}C_{0}^{2}\unicode[STIX]{x1D713}_{0}-C_{2}^{2}\unicode[STIX]{x1D713}_{0}-C_{3}^{2}\unicode[STIX]{x1D713}_{0}^{3}.\end{eqnarray}$$

Similar substitutions may be implemented in the density relation to unravel

(3.16) $$\begin{eqnarray}O(M_{0}^{2}):\unicode[STIX]{x1D70C}_{1}=-\frac{1}{2r^{2}}\left[\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}\right)^{2}+B_{0}^{2}\unicode[STIX]{x1D713}_{0}^{2}\right];\end{eqnarray}$$

(3.17) $$\begin{eqnarray}\displaystyle O(M_{0}^{4}):\unicode[STIX]{x1D70C}_{2} & = & \displaystyle -\frac{2+\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D70C}_{1}^{2}-\frac{\unicode[STIX]{x1D70C}_{1}}{r^{2}}[C_{0}^{2}\unicode[STIX]{x1D713}_{0}^{2}+C_{1}^{2}]\nonumber\\ \displaystyle & & \displaystyle -\,\frac{\unicode[STIX]{x1D713}_{0}}{2r^{2}}[2C_{0}^{2}\unicode[STIX]{x1D713}_{1}+C_{2}^{2}\unicode[STIX]{x1D713}_{0}^{2}+C_{3}^{2}\unicode[STIX]{x1D713}_{0}^{3}]-\frac{1}{r^{2}}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}r}\right].\qquad\end{eqnarray}$$

In seeking a compressible mean flow approximation, our procedure consists of solving (3.13), (3.14), (3.16) and (3.17) in this staggered sequence. A flowchart describing this process is posted as figure 2.

Figure 2. Flowchart for the density–streamfunction formulation needed to obtain a compressible Bragg–Hawthorne solution.

4.1 Geometric idealization

As alluded to earlier, the motivation for the present study centres on the bidirectional vortex flow field that is relevant to a number of cyclonic-flow applications, including the VCCWC engine prototype. In its basic form, the VCCWC may be modelled as a closed–open cylinder of radius $a$ and height $L$ . When all spatial coordinates are normalized by the radius, an idealized chamber of unit radius emerges, as depicted in figure 3. The origin of the coordinate system may be conveniently placed at the centre of the headwall with the dimensionless radial and axial coordinates being $r$ and $z$ , respectively. The chamber contains a confined vortex that is dominated by its swirl velocity, $U$ ; the latter may be used to normalize all velocities to the extent of producing a unit tangential velocity at the inlet $(v=1)$ . In practice, this average inflow velocity is induced by tangential injectors that are evenly distributed around the engine base, just upstream of the nozzle. Due to unavoidable collisions at the boundaries, the fluid develops both axial and radial components of velocity as it swirls towards the headwall while occupying the outer, annular region. When the flow reaches the headwall, it reverses axial direction and returns along the chamber bore until it discharges through the partially open base of radial fraction  $\unicode[STIX]{x1D6FD}$ . As usual, a spinning, non-translating mantle separates the two vortex regions. In its simplest mathematical idealization, two geometric parameters emerge in the analysis, namely, the aspect ratio $l=L/a$ and the open radial fraction at the base, $\unicode[STIX]{x1D6FD}=b/a$ . The latter is often taken to coincide with the mantle location in such a manner to allow the outflow diameter to match that of the outlet cross-section. The resulting arrangement is captured in figure 3 where the broken lines of radius $\unicode[STIX]{x1D6FD}$ delineate the interfacial boundary between the inner, core and outer, annular vortex regions. Such an optimal configuration mitigates the onset of flow collisions and the establishment of undesirable recirculation zones.

Figure 3. Dimensionless coordinates and key parameters associated with the mathematical idealization of the bidirectional vortex chamber.

4.2 Analytical model

In modelling the bidirectional vortex, the following boundary conditions are customarily adopted: (i) no axial flow at the headwall, (ii) no radial flow across the centreline, (iii) no radial flow at the sidewall and (iv) a mass balance between the inflow and the outflow. As depicted on the graph in dimensionless form, these conditions translate into

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle w(r,0)=0;\quad \frac{1}{\unicode[STIX]{x1D70C}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r}=0; & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle u(0,z)=0;\quad -\frac{1}{\unicode[STIX]{x1D70C}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}=0; & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle u(1,z)=0;\quad -\frac{1}{\unicode[STIX]{x1D70C}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}=0 & \displaystyle\end{eqnarray}$$

and

(4.4) $$\begin{eqnarray}\displaystyle {\dot{m}}_{i}=2\unicode[STIX]{x03C0}\int _{0}^{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D70C}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}r\,\text{d}r. & & \displaystyle\end{eqnarray}$$

To remain consistent with the Rayleigh–Janzen perturbation technique, we expand the boundary conditions in terms of the Mach number squared. At leading order, we collect

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle w^{(0)}(r,0)=0;\quad \frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}=0; & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle u^{(0)}(0,z)=0;\quad -\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}=0; & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle u^{(0)}(1,z)=0;\quad -\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}=0 & \displaystyle\end{eqnarray}$$

and

(4.8) $$\begin{eqnarray}Q_{i}=2\unicode[STIX]{x03C0}\int _{0}^{\unicode[STIX]{x1D6FD}}w^{(0)}r\,\text{d}r=2\unicode[STIX]{x03C0}\int _{0}^{\unicode[STIX]{x1D6FD}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}\,\text{d}r.\end{eqnarray}$$

In the above, the boundary conditions used by Majdalani (Reference Majdalani2012) to derive the incompressible Trkalian vortex are restored identically.

In conformance with perturbation theory, the constraints at the first compressible order must not unduly influence the incompressible motion. To ensure that the solution remains valid $\forall \;M_{0}$ , a set of homogeneous boundary conditions may hence be imposed on all upper-level approximations. Furthermore, implementation of these conditions will have to be carried out while taking into account the additional terms that are produced from their higher-order expansions. For the bidirectional vortex, the first-order counterparts to (4.1)–(4.4) may be expressed as:

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle w^{(1)}(r,0)=0:\quad \frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}r}-\frac{\unicode[STIX]{x1D70C}_{1}}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}=0\quad (\text{impermeable headwall}); & \displaystyle\end{eqnarray}$$

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle u^{(1)}(0,z)=0:\quad \frac{\unicode[STIX]{x1D70C}_{1}}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}-\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z}=0\quad (\text{centreline symmetry}); & \displaystyle\end{eqnarray}$$

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle u^{(1)}(1,z)=0:\quad \frac{\unicode[STIX]{x1D70C}_{1}}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}-\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z}=0\quad (\text{impermeable sidewall}) & \displaystyle\end{eqnarray}$$

and

(4.12) $$\begin{eqnarray}2\unicode[STIX]{x03C0}\int _{0}^{\unicode[STIX]{x1D6FD}}[\unicode[STIX]{x1D70C}_{1}w^{(0)}+w^{(1)}]r\,\text{d}r=0\quad (\text{mass conservation}).\end{eqnarray}$$

These relations represent the first-order constraints for the compressible, bidirectional flow field. Should higher-order corrections be required, homogeneous constraints that resemble those at the first order could be generated at $O(M_{0}^{4})$ , $O(M_{0}^{6})$ , and so on.

5.1 Leading-order solution

The leading-order streamfunction must be consistent with the incompressible solution for the same problem. In this spirit, equation (3.13) may be treated with separation of variables. Assuming $\unicode[STIX]{x1D713}_{0}=f(r)g(z)$ , equation (3.13) becomes

(5.1) $$\begin{eqnarray}-\frac{g^{\prime \prime }}{g}=\frac{1}{f}\left(f^{\prime \prime }-\frac{1}{r}f^{\prime }+C_{0}^{2}f\right)=\unicode[STIX]{x1D708}^{2},\end{eqnarray}$$

where $\unicode[STIX]{x1D708}^{2}$ can be positive, negative or zero. Depending on the value chosen for $\unicode[STIX]{x1D708}^{2}$ , three solutions may be conceived, namely,

(5.2) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{0}=\left\{\begin{array}{@{}ll@{}}r(k_{1}z+k_{2})[k_{3}J_{1}(C_{0}r)+k_{4}Y_{1}(C_{0}r)];\quad & \unicode[STIX]{x1D708}^{2}=0\\ r[k_{1}\sin (\unicode[STIX]{x1D708}z)+k_{2}\cos (\unicode[STIX]{x1D708}z)]\left[k_{3}J_{1}(r\sqrt{C_{0}^{2}-\unicode[STIX]{x1D708}^{2}})+k_{4}Y_{1}(r\sqrt{C_{0}^{2}-\unicode[STIX]{x1D708}^{2}})\right];\quad & C_{0}^{2}>\unicode[STIX]{x1D708}^{2}\\ r[k_{1}\sinh (\unicode[STIX]{x1D708}z)+k_{2}\cosh (\unicode[STIX]{x1D708}z)]\left[k_{3}J_{1}(r\sqrt{C_{0}^{2}+\unicode[STIX]{x1D708}^{2}})+k_{4}Y_{1}(r\sqrt{C_{0}^{2}+\unicode[STIX]{x1D708}^{2}})\right];\quad & C_{0}^{2}<\unicode[STIX]{x1D708}^{2}.\end{array}\right.\end{eqnarray}$$

In reality, the last two variations prove to be equivalent as one can be reproduced from the other by simply replacing $\unicode[STIX]{x1D708}$ with $\pm \text{i}\unicode[STIX]{x1D708}$ . When accounting for the imaginary part, the hyperbolic functions reduce to their regular trigonometric counterparts as arguments of the Bessel functions become identical when the imaginary $\unicode[STIX]{x1D708}^{2}$ switches its sign. For brevity, the remainder of this study will focus on the axially linear case as a vehicle for developing a compressible approximation.

5.2 Leading-order boundary conditions

To satisfy the centreline boundary conditions for all values of $z$ , we set $k_{4}=0$ everywhere. Furthermore, applying (4.1) leads to

(5.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}(r,0)}{\unicode[STIX]{x2202}r}=k_{2}k_{3}C_{0}J_{0}(C_{0}r)=0.\end{eqnarray}$$

Since equating either $k_{3}$ or $C_{0}$ to zero results in a trivial outcome, we take $k_{2}=0$ . Substituting the resultant streamfunction back into the sidewall boundary condition produces

(5.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}(1,z)}{\unicode[STIX]{x2202}z}=k_{1}k_{3}J_{1}(C_{0}).\end{eqnarray}$$

For (5.4) to be true $\forall z$ , $C_{0}$ must be a root of the Bessel function of the first kind, or

(5.5) $$\begin{eqnarray}C_{0}=\unicode[STIX]{x1D706}_{n};\quad n=0,1,2,\ldots .\end{eqnarray}$$

Increments in $n$ will effectively trigger an increasing number of axial reversals in the flow, specifically $(n+1)$ reversals. In practice, only an odd number of reversals will be applicable to the problem at hand and so, to recover the standard bidirectional vortex model of Majdalani (Reference Majdalani2012), we restrict our analysis to the $n=0$ case.

At this juncture, we are left with the lumped constant $k_{1}k_{3}$ that must be determined by matching the inflow and outflow mass fluxes. At the leading order, this may be written as

(5.6) $$\begin{eqnarray}2\unicode[STIX]{x03C0}\int _{0}^{\unicode[STIX]{x1D6FD}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}r\,\text{d}r=2\unicode[STIX]{x03C0}\int _{0}^{\unicode[STIX]{x1D6FD}}w^{(0)}(r,L)r\,\text{d}r=\frac{1}{\unicode[STIX]{x1D70E}}.\end{eqnarray}$$

From the mass balance in (5.6) we deduce

(5.7) $$\begin{eqnarray}k_{1}k_{3}=\frac{1}{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}LJ_{1}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}=\frac{\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D6FD}J_{1}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})};\quad \unicode[STIX]{x1D705}\equiv \frac{1}{2\unicode[STIX]{x03C0}L\unicode[STIX]{x1D70E}}.\end{eqnarray}$$

This converts the streamfunction form into

(5.8) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{0}=z\frac{r\unicode[STIX]{x1D705}J_{1}(\unicode[STIX]{x1D706}_{0}r)}{\unicode[STIX]{x1D6FD}J_{1}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}.\end{eqnarray}$$

As to be expected from a leading-order approximation, equation (5.8) reproduces the incompressible Trkalian profiles obtained by Majdalani (Reference Majdalani2012) in a right-cylindrical cyclone. With the streamfunction solution being fully determined, the density may be retrieved from (3.16) following a straightforward substitution. After some manipulations, the density may be presented as:

(5.9) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{1}=-\frac{\unicode[STIX]{x1D705}^{2}}{2\unicode[STIX]{x1D6FD}^{2}J_{1}^{2}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}\{J_{1}^{2}(\unicode[STIX]{x1D706}_{0}r)+z^{2}\unicode[STIX]{x1D706}_{0}^{2}[J_{0}^{2}(\unicode[STIX]{x1D706}_{0}r)+J_{1}^{2}(\unicode[STIX]{x1D706}_{0}r)]\}.\end{eqnarray}$$

Equation (5.9) is quite illuminating. In fact, it confirms the need for higher powers of $\unicode[STIX]{x1D713}$ within the fundamental definition of $B$ in (3.9). Clearly, Bessel functions that are elevated to the second power appear thrice in the density. These, in turn, multiply a single Bessel function in (2.13), the first-order streamfunction relation. Therefore, in seeking appropriate candidate functions for the particular solution, terms that may be expressed in multiples of three Bessel functions must be attempted. This step is prompted by the requirement to write $B$ in terms of $\unicode[STIX]{x1D713}$ at the basis of the Bragg–Hawthorne procedure.

5.3 First-order streamfunction solution

The first-order correction follows a similar roadmap, albeit with increased complexity. Instead of a homogeneous equation, a particular solution must be determined in such a way to accommodate the terms appearing on the right-hand side of (2.13). For the axially linear case, suitable substitutions of $\unicode[STIX]{x1D713}_{0}$ and $\unicode[STIX]{x1D70C}_{1}$ lead to the first-order, compressible Bragg–Hawthorne equation, namely,

(5.10) $$\begin{eqnarray}\displaystyle \text{D}^{2}\unicode[STIX]{x1D713}_{1}+\unicode[STIX]{x1D706}_{0}^{2}\unicode[STIX]{x1D713}_{1} & = & \displaystyle \frac{z\unicode[STIX]{x1D705}^{3}\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}J_{1}(\unicode[STIX]{x1D706}_{0}r)(\!J_{0}(\unicode[STIX]{x1D706}_{0}r)J_{1}(\unicode[STIX]{x1D706}_{0}r)-2\unicode[STIX]{x1D706}_{0}rJ_{0}^{2}(\unicode[STIX]{x1D706}_{0}r)\nonumber\\ \displaystyle & & \displaystyle +\,z^{2}\unicode[STIX]{x1D706}_{0}^{2}\{J_{0}(\unicode[STIX]{x1D706}_{0}r)J_{1}(\unicode[STIX]{x1D706}_{0}r)+\unicode[STIX]{x1D706}_{0}r[J_{0}^{2}(\unicode[STIX]{x1D706}_{0}r)+J_{1}^{2}(\unicode[STIX]{x1D706}_{0}r)]\}\!)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{rz\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D6FD}J_{1}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}C_{2}^{2}J_{1}(\unicode[STIX]{x1D706}_{0}r)-\frac{r^{3}z^{3}\unicode[STIX]{x1D705}^{3}}{\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}C_{3}^{2}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}r).\end{eqnarray}$$

Rather than a standard separation of variables approach, we employ an ansatz that is guided by the content of the non-homogeneous terms. Recognizing that both $z$ and $z^{3}$ appear in (5.10), we let

(5.11) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{1}=zR_{a}+z^{3}R_{b}.\end{eqnarray}$$

By virtue of (5.11), our single PDE gives rise to two ODEs with one-way coupling through

(5.12) $$\begin{eqnarray}\displaystyle z^{3}:\;R_{b}^{\prime \prime }-\frac{1}{r}R_{b}^{\prime }+\unicode[STIX]{x1D706}_{0}^{2}R_{b} & = & \displaystyle -C_{3}^{2}\frac{\unicode[STIX]{x1D705}^{3}}{\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}r^{3}J_{1}^{3}(r\unicode[STIX]{x1D706}_{0})+\frac{\unicode[STIX]{x1D705}^{3}\unicode[STIX]{x1D706}_{0}^{3}}{\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}\nonumber\\ \displaystyle & & \displaystyle \times \,\{J_{0}(r\unicode[STIX]{x1D706}_{0})J_{1}(r\unicode[STIX]{x1D706}_{0})+r\unicode[STIX]{x1D706}_{0}[J_{0}^{2}(r\unicode[STIX]{x1D706}_{0})+J_{1}^{2}(r\unicode[STIX]{x1D706}_{0})]\};\end{eqnarray}$$

(5.13) $$\begin{eqnarray}\displaystyle z:\;R_{a}^{\prime \prime }-\frac{1}{r}R_{a}^{\prime }+\unicode[STIX]{x1D706}_{0}^{2}R_{a}+6R_{b} & = & \displaystyle \frac{\unicode[STIX]{x1D705}^{3}\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}J_{1}(r\unicode[STIX]{x1D706}_{0})[J_{0}(r\unicode[STIX]{x1D706}_{0})J_{1}(r\unicode[STIX]{x1D706}_{0})-2\unicode[STIX]{x1D706}_{0}rJ_{1}(r\unicode[STIX]{x1D706}_{0})]\nonumber\\ \displaystyle & & \displaystyle -\,C_{2}^{2}\frac{\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D6FD}J_{1}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}rJ_{1}(r\unicode[STIX]{x1D706}_{0}).\end{eqnarray}$$

Our next step is to first solve (5.12), being a sole function of $R_{b}$ . The ensuing solution may be then substituted back into (5.13) to extract $R_{a}$ and, with it, a complete compressible correction.

In practice, the solution to (5.12) is exacerbated by its dependence on $J_{0}^{3}(x)$ and $J_{1}^{3}(x)$ terms. While Bessel function integrals remain straightforward to evaluate in closed form, integrals for multiplicative Bessel functions can be elusive. In lieu of a completely analytical closure, our correction becomes limited to a semi-analytical formulation that requires the numerical evaluation of a handful of integrals. To overcome this difficulty, the integrals themselves will be isolated and specified as functions that may be differentiated or integrated at will, so that the boundary conditions can still be determined analytically. In essence, these new integrals may be viewed as special functions that enable us to retain the analytical character of our formulation. After some effort, the $z^{3}$ multiplier is found to be

(5.14) $$\begin{eqnarray}\displaystyle R_{b} & = & \displaystyle rJ_{1}(\unicode[STIX]{x1D706}_{0}r)\left[\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}^{3}}{2\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}(C_{3}^{2}I_{2}-\unicode[STIX]{x1D706}_{0}^{3}I_{1})+k_{7}\right]\nonumber\\ \displaystyle & & \displaystyle +\,rY_{1}(\unicode[STIX]{x1D706}_{0}r)\left[\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}^{3}}{2\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}(\unicode[STIX]{x1D706}_{0}^{3}I_{3}-C_{3}^{2}I_{4})+k_{8}\right].\end{eqnarray}$$

Here $k_{7}$ and $k_{8}$ represent integration constants while $r_{1}$ and $r_{2}$ denote variable substitutions in the radial integrals; furthermore, $I_{n}$ designates the $n$ th integral in the first-order solution. For the reader’s convenience, these are defined in appendix A.

The $z$ multiplier may be obtained along similar lines. Inserting (5.14) into (5.13) yields

(5.15) $$\begin{eqnarray}\displaystyle R_{a} & = & \displaystyle rJ_{1}(\unicode[STIX]{x1D706}_{0}r)\left[\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}}{2\unicode[STIX]{x1D6FD}J_{1}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}\left(\frac{\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x1D6FD}^{2}J_{1}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}I_{5}+C_{2}^{2}I_{6}\right)+3\unicode[STIX]{x03C0}I_{7}+k_{5}\right]\nonumber\\ \displaystyle & & \displaystyle +\,rY_{1}(\unicode[STIX]{x1D706}_{0}r)\left[\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}}{2\unicode[STIX]{x1D6FD}J_{1}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}\left(\frac{\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x1D6FD}^{2}J_{1}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}\unicode[STIX]{x1D706}_{0}I_{8}-C_{2}^{2}I_{9}\right)-3\unicode[STIX]{x03C0}I_{10}+k_{6}\right].\end{eqnarray}$$

By substituting (5.14) and (5.15) back into (5.11), one arrives at the general compressible correction. To make further headway, an appropriate set of boundary conditions will be considered and discussed.

5.4 First-order boundary conditions

Compared to the leading order, the boundary conditions at the first order change slightly. In fact, ensuring that the compressible correction does not unduly influence the solution warrants the use of homogeneous constraints. Because our boundary conditions are written in terms of the velocity, it is useful to revisit the expanded velocity–streamfunction relationship. At the first order, we recover

(5.16a,b ) $$\begin{eqnarray}\displaystyle u^{(1)}=\frac{\unicode[STIX]{x1D70C}_{1}}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}z}-\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z};\quad w^{(1)}=\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}r}-\frac{\unicode[STIX]{x1D70C}_{1}}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}r}. & & \displaystyle\end{eqnarray}$$

To avoid lengthy streamfunction expressions, we omit the general expansion of (5.16). Instead, we examine each boundary condition individually. Examining the sidewall condition on the radial velocity we find

(5.17) $$\begin{eqnarray}\displaystyle & & \displaystyle z^{2}[6k_{8}\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})Y_{1}(\unicode[STIX]{x1D706}_{0})+3\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}^{3}\unicode[STIX]{x1D706}_{0}^{3}Y_{1}(\unicode[STIX]{x1D706}_{0})I_{3}(1)-3C_{3}^{2}\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}^{3}Y_{1}(\unicode[STIX]{x1D706}_{0})I_{4}(1)]\nonumber\\ \displaystyle & & \displaystyle \quad +\, [\!2k_{6}\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})Y_{1}(\unicode[STIX]{x1D706}_{0})-6\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})Y_{1}(\unicode[STIX]{x1D706}_{0})I_{10}(1)\nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}^{3}\unicode[STIX]{x1D706}_{0}Y_{1}(\unicode[STIX]{x1D706}_{0})I_{8}(1)-C_{2}^{2}\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FD}^{2}\unicode[STIX]{x1D705}J_{1}^{2}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})Y_{1}(\unicode[STIX]{x1D706}_{0})I_{9}(1)\!]=0,\end{eqnarray}$$

which must hold for all $z$ . Additionally, since all $I_{n}(1)$ must vanish identically, we deduce that $k_{6}=k_{8}=0$ . To find $k_{7}$ , we expand $I_{7}$ and $I_{10}$ ; we then apply the centreline condition to obtain

(5.18) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}6\unicode[STIX]{x03C0}[k_{7}I_{10a}(0)+I_{10b}(0)]+{\displaystyle \frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D6FD}^{3}J_{1}^{3}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})}}[C_{2}^{2}\unicode[STIX]{x1D6FD}^{2}\unicode[STIX]{x1D705}J_{1}^{2}(\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D6FD})I_{9}(0)-\unicode[STIX]{x1D705}^{3}\unicode[STIX]{x1D706}_{0}I_{8}(0)]=0;\\ C_{3}^{2}=\unicode[STIX]{x1D706}_{0}^{3}{\displaystyle \frac{I_{3}(0)}{I_{4}(0)}}.\end{array}\right\}\end{eqnarray}$$

Along similar lines, we may expand the continuity equation at first order and collect

(5.19) $$\begin{eqnarray}\int _{0}^{\unicode[STIX]{x1D6FD}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}r}\,\text{d}r=0.\end{eqnarray}$$

The detailed form of the above expression is long and, as such, of minimal interest to the reader. However, the remaining integral may be straightforwardly handled using symbolic programming. The evaluation of (5.19) completes our analysis of the first-order compressible streamfunction from which all other flow parameters may be retrieved.

6.1 Comparison to a numerical simulation of a similar cyclonic chamber

Before presenting and discussing our findings, a simple comparison is carried out using a fully nonlinear Navier–Stokes simulation. To this end, the main chamber in a right circular cylinder is used to define the computational domain with a length of $L=508$  mm and a diameter $2a$ of 457.2 mm, thus resulting in a chamber with a practical aspect ratio $L/a$ of 2.22. In order to realistically approximate the injection conditions, the purely tangential injection of the analytical model is replaced with a circular array of eight 38.1 mm injectors, tangentially distributed around the base of the chamber. The chamber’s outlet diameter is set to 320 mm before transitioning into a conical nozzle. As shown in figure 4, the mesh that is used consists of 1 078 632 cells, which comprise 2 383 058 faces, 324 966 nodes and 31 partitions.

Figure 4. The mesh structure and computational domain using (a) isometric, (b) top, (c) side and (d) bottom views.

A robust, finite-volume, computational fluid dynamics (CFD) solver is chosen to perform the attendant simulation. The inlet boundary condition that specifies the velocity normal to the injector face is set to a value of $44.8~\text{m}~\text{s}^{-1}$ , thus corresponding to an injection Mach number of $M_{0}=0.1$ . At the outlet, a pressure boundary condition is set to 10 kPa. The SIMPLE scheme is employed with a Green–Gauss cell based discretization for gradients and first-order upwind discretizations for the remaining terms. Furthermore, convergence is established through an examination of residual checks, with the lowest being a residual of $10^{-6}$ or lower in the energy equation, as well as a visual inspection of basic flow parameters, such as the pressure and swirl velocities, to ensure their proper development. Fundamentally, the model assumes laminar, compressible and ideal gas conditions with standard air as the working fluid. A strictly inviscid simulation could not be achieved in this swirl-driven configuration, mainly because of the classic singularity along the centreline, which could not be handled by the numerical solver without including the effects of viscosity, no matter how small.

A comparison of the axial, radial, swirl and total velocity profiles is presented in figure 5. In all four cases, a qualitative agreement may be seen between the numerically and analytically obtained velocity shapes. In fact, the agreement between the present formulation and computations appears to improve in the core region of the cyclonic vortex, while deteriorating in the wall region. The discrepancy near the sidewall may be attributed to the inability of an inviscid solution to satisfy the no-slip requirement at the solid boundary.

Figure 5. Comparison between analytical predictions and numerical simulations for the (a) radial velocity $u$ , (b) tangential velocity $v$ , (c) axial velocity $w$ and (d) total velocity for a Mach number of 0.1.

By moving our attention to the thermodynamic variables, both pressure and density distributions shown in figure 6 seem to exhibit similar qualitative agreement with the simulated data. In both cases, the agreement is seen to improve near the walls, an artefact that may be attributed to the two solutions being normalized by their peak pressure and density values. In comparing the density to the pressure curves, the latter become more tightly grouped near the walls. Despite their inevitable differences, the overall shapes of the numerically computed pressure and density profiles remain similar to their analytical counterparts, as both exhibit a sharp drop as the centreline of the chamber is approached.

Several reasons may be causing the observed disparities. Firstly, the present study remains limited to an inviscid model, while the numerical simulation incorporates viscous effects. As alluded to earlier, excluding viscosity in the computational model leads to an unbounded swirl velocity at the centreline, thus making an accurately converged solution unattainable. Secondly, the purely circumferential line injection assumed in the analytical framework can only be approximated in the computational model using eight finite injectors that are evenly arrayed around the base of the chamber. While this arrangement leads to a more practical injection method, the dissimilarities with the analytical model become more appreciable near the wall. In fact, the effect of using a discrete number of injectors becomes also noticeable near the aft end of the chamber, where the dimensional axial coordinate approaches $L$ . Thirdly, the analytical approximation incorporates a single, asymptotically linearized, compressible correction to serve as a bulk estimate for dilatational effects, unlike the Navier–Stokes solver, which incorporates all nonlinear interactions. It is therefore plausible that higher-order corrections may be required to improve the model’s capabilities. Finally, the configuration used in this work consists of a straight cylinder, whereas the chamber geometry used by the solver entails a nozzle attachment. Evidently, differences can be expected in the vicinity of the nozzle, which is needed in the computational domain to promote convergence.

Although the qualitative agreement in figures 5 and 6 is encouraging, it is only meant to be serve as a simple check that the analytical model exhibits similar features to those obtained from a simulated chamber, albeit slightly different from the idealized geometric configuration adopted within the mathematical framework. Clearly, a more detailed endeavour will be required to fully characterize the bidirectional vortex flow motion computationally, namely, an effort that falls beyond the scope of this work.

Figure 6. A comparison between analytical and numerical simulations for the (a) pressure and (b) density at three axial locations and a Mach number of 0.1.

6.2 Compressible velocity steepening

Pursuant to the streamfunction determination, the compressible motion may be characterized in all three spatial directions. To avoid unnecessary collisions and potential recirculation sites, the open fraction at the base, $\unicode[STIX]{x1D6FD}$ , may be conveniently equated to the dimensionless mantle radius, with the effect of allowing the outgoing stream to exit the chamber unobstructed. To illustrate the effects of compressibility, we examine a range of $\unicode[STIX]{x1D705}$ extending from 0.1 to 0.4. To facilitate comparisons relative to previous studies, the aspect ratio of the chamber is taken to be $l=4/3$ , while typical values of $\unicode[STIX]{x1D6FE}=\{1.2,1.6\}$ are assigned to the ratio of specific heats. As for the injection Mach number, we anchor our analysis around $M_{0}=0.1$ and 0.2, being two commonly used values in propulsive applications.

Figure 7. Axial velocity variation in the exit plane showcasing (a) the actual profile for $\unicode[STIX]{x1D705}=0.3$ and (b) the spatial distribution of its compressible correction.

Thus motivated by the need to characterize the VCCWC bulk flow field, we begin with the axial velocity, $w$ , which drives engine performance after expansion. In figure 7(a), we consider the axial profile at the chamber exit, $z=l$ , for reference Mach numbers of 0.1 and 0.2 with $\unicode[STIX]{x1D705}=0.3$ . At $M_{0}=0.1$ , the compressible contribution seems to induce a minor although still visible variation in the velocity profile; however, by increasing the injection speed to $M_{0}=0.2$ , a more appreciable effect is realized. In comparison to the compressible Taylor–Culick solution by Majdalani (Reference Majdalani2007a ), which is derived in the context of solid rocket motor internal motion, the first-order velocity corrections also display a weak dependence on $\unicode[STIX]{x1D6FE}$ . Along similar lines, the axial velocity exhibits a steepening effect that intensifies with the growth of the injection Mach number. When this occurs, the polarity transition that accompanies mantle formation acquires a blunter slope even as $w$ crosses the radial axis vertically. This finding is consistent with the axial steepening observed in compressible models of solid rocket motors (e.g. Traineau, Hervat & Kuentamann Reference Traineau, Hervat and Kuentamann1986; Balakrishnan, Liñán & Williams Reference Balakrishnan, Liñán and Williams1992; Wasistho, Balachandar & Moser Reference Wasistho, Balachandar and Moser2004; Majdalani Reference Majdalani2007a ; Maicke & Majdalani Reference Maicke and Majdalani2008; Akiki & Majdalani Reference Akiki and Majdalani2012, Reference Akiki and Majdalani2016). Contrary to the aforementioned studies, no net amplification of the axial velocity may be noted here, aside from a reshaping of the profile itself. As the conservation constraint at the exit must be identically satisfied at the leading order, the resulting mass exiting the chamber at the first order must be self-cancelling when integrated over the flow cross-section. This requirement compels the morphing of the velocity contour without affecting the overall mass flux. As for the compressible contribution $w(1)$ depicted in figure 7(b), it is featured for $\unicode[STIX]{x1D705}=0.3$ , a value that has been used in previous studies in the context of the VCCWC engine simulation. Based on the Bragg–Hawthorne equation framework, the first-order correction becomes more pronounced at higher values of $\unicode[STIX]{x1D705}$ , $\forall M_{0}$ . Additionally, it can be seen that the compressible contribution vanishes at two distinct points, namely, at approximately $r=0.2$ and 0.75. These sites derive their locations from the mass balance relation which, when applied to the compressible correction, gives rise to two polarity nodes in the axial velocity where a zero net flux may be achieved across the chamber cross-section.

Figure 8. Radial velocity variation in the exit plane showcasing (a) the actual profile at $\unicode[STIX]{x1D705}=0.3$ and (b) the spatial distribution of its compressible correction.

By virtue of continuity and momentum balances, elements of the steepening mechanism observed in $w$ are transferred to the radial velocity, $u$ , which is illustrated in figure 8(a). For small deviations from the incompressible case of $M_{0}=0$ , the solution seems to be fairly well guided by the shape of the unperturbed solution. In the present case, the effect of compression causes a spatial shifting of the peak magnitude in $u$ towards $r=1$ . As the Mach number is further increased to $M_{0}=0.2$ , the outward shift in peak amplitudes is accompanied by a more visible increase $\Vert u_{max}\Vert$ beyond its incompressible value. This particular amplification of $\Vert u\Vert$ in the vicinity of the sidewall can be so pronounced that it must be offset by an appropriate attenuation in the radial velocity within the core region. The corresponding shift in $u$ that is experienced near the centreline causes the radial velocity profile to switch polarity while returning to $r=0$ . Mathematically, because the radial velocity is written as a $z$ -derivative of the streamfunction, it will be strongly influenced by the reversing nature imposed by the conservation principle in the exit plane. Here too, the compressible radial contribution vanishes at $r=0.55$ as clearly depicted in figure 8(b). Even for $\unicode[STIX]{x1D705}=0.3$ , the total radial velocity remains significantly smaller than the axial or tangential velocities, and this behaviour may be attributed to the sidewalls being non-injecting. Nonetheless, the compressible correction itself becomes of the same order in both axial and radial directions, hence leading to a proportionately larger effect on the radial velocity. The ensuing trend is reversed, however, when $\unicode[STIX]{x1D705}$ is reduced in a manner to mitigate the actual compressible contribution.

The third and most prominent component of the compressible velocity is showcased in figure 9(a) for the two representative Mach numbers and $\unicode[STIX]{x1D705}$ . In comparison to $u$ and $w$ , the behaviour of the swirl velocity $v$ seems to mimic that of its radial counterpart; its profile draws closer to the sidewall with each additional growth in $M_{0}$ . Here too, the maximum swirling speed increases at higher values of the Mach number. Note that in figure 9(b), only the compressible correction is featured at $\unicode[STIX]{x1D705}=0.3$ . In both cases, the compressible correction vanishes at $r=0.5$ . A closer look at $v^{(1)}$ reveals that its shape resembles that of the radial profile, except for being of opposite sign so long as $0<r<1$ . Mathematically, differences in magnitudes between $u$ and $v$ may be attributed to the reduced $z$ dependence of the radial velocity; unlike $u,$ the swirl velocity retains a $z^{3}$ dependence through the streamfunction expression.

Figure 9. Tangential velocity variation in the exit plane showcasing (a) the actual profile at $\unicode[STIX]{x1D705}=0.3$ and (b) the spatial distribution of its compressible correction.

6.3 Compressible sliding of the mantle interface

Another interesting feature of the compressible solution stems from its mantle gaining an axial dependence that cannot be accounted for by the incompressible model. At the leading order, the mantle maintains a constant radial position for all values of $z$ at approximately $r=0.627$ . In the compressible model, the mantle location gains a $z$ dependence that is clearly illustrated in figure 10 at several representative values of the Mach number and $\unicode[STIX]{x1D705}$ . Accordingly, the mantle is seen to shift outwardly in the proximity of the headwall $(z=0)$ , with the shift widening at larger values of $\unicode[STIX]{x1D705}$ as well as the reference Mach number. Then as the fluid travels towards the exit plane, the mantle continues to slide outwardly, almost linearly in $z$ . Such behaviour appears to be consistent with the theoretical findings of Maicke & Majdalani (Reference Maicke and Majdalani2008), although the two studies are based on dissimilar compressible flow models. At first glance, the linear character of the mantle translation away from the headwall may be viewed as somewhat perplexing because of the solution’s explicit dependence on $z^{3}$ . In the present situation, however, the low aspect ratio of $l=4/3$ may be responsible for the linear behaviour up to the exit plane. In longer chambers, it is quite likely for the linear behaviour to become superseded by a cubic dependence, especially in the presence of sufficiently large reference Mach numbers.

In what concerns experimental evidence, Smith (Reference Smith1962a ,Reference Smith b ) reports two antithetical cases: one in which the mantle slides inwardly, towards the centreline, as the distance from the headwall is increased, and one expanding outwardly, towards the sidewall. Without further scrutiny, it may only be possible to speculate over the factors leading to mantle variability. For example, it may be conceivable for viscous effects to compete with compressibility to the extent of one overpowering the other in a given configuration. It is more likely, however, for the geometric design, along with its inlet and outlet arrangements, including the presence of a vortex finder, to influence the stable position of the mantle interface between the inner and outer vortex regions (Akiki & Majdalani Reference Akiki and Majdalani2010). The presence of a protrusion into the flow, such as the submerged vortex finder appearing in experiments by Smith (Reference Smith1962a ,Reference Smith b ), may have an appreciable bearing on the final mantle location.

Figure 10. Compressible mantle location for (a) $M_{0}=0.1$ and (b) 0.2 using $\unicode[STIX]{x1D705}=0.1$ , 0.2, 0.3, and 0.4. The vertical line corresponds to the fixed, incompressible mantle site.

6.4 Density and pressure variations

In the compressible Bragg–Hawthorne framework, all thermodynamic quantities may be restored from the density. In view of the isentropic relation used at the basis of the density–streamfunction formulation, the pressure and temperature may be straightforwardly deduced from the density. With this in mind, retrieving and characterizing the compressible density correction is paramount to the determination and interpretation of the corresponding pressure and temperature fields. Both $p_{1}$ and $T_{1}$ differ from $\unicode[STIX]{x1D70C}_{1}$ by a constant multiplier, namely,

(6.1a,b ) $$\begin{eqnarray}p_{1}=\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D70C}_{1}\quad \text{and}\quad T_{1}=(\unicode[STIX]{x1D6FE}-1)\unicode[STIX]{x1D70C}_{1}.\end{eqnarray}$$

In figure 11, density variations are shown in the exit plane for the two representative injection Mach numbers of $M_{0}=0.1$ and 0.2, using $\unicode[STIX]{x1D705}=0.1$ , 0.2, 0.3, and 0.4. As it may be surmised from the graphs, the density appears to be sensitive to both variations in the Mach number and $\unicode[STIX]{x1D705}$ . However, the sensitivity to $\unicode[STIX]{x1D705}$ is amplified substantially when the Mach number is incremented from 0.1 to 0.2. This may be attributed to the former being closer to the leading-order benchmark than the latter. In both cases, however, as the Mach number and $\unicode[STIX]{x1D705}$ are augmented, the normalized density undergoes successive decreases throughout the chamber, with the most significant depreciation occurring along the centreline. It is this particular drop in density that drives, in part, the variation in the axial velocity correction at the first order through its contribution to the mass conservation requirement at $z=l$ .

Figure 11. Density distribution for injection Mach numbers of (a) $M_{0}=0.1$ and (b) 0.2 using $z=l$ and $\unicode[STIX]{x1D705}=0.1$ , 0.2, 0.3 and 0.4.

As for the pressure companion, which regains a dependence on $\unicode[STIX]{x1D6FE}$ , similar trends are depicted in figure 12, where the dimensionless pressure distribution is displayed for $\unicode[STIX]{x1D705}=0.1$ and 0.3 as well as $\unicode[STIX]{x1D6FE}$ values of 1.2 and 1.6. Here too, the largest depreciation in the pressure is realized near the centreline, and this effect is accentuated at higher values of $M_{0}$ or $\unicode[STIX]{x1D6FE}$ . In the $M_{0}=0.2$ case, the compressible correction causes the pressure near the centreline to drop precipitously, thus leading to low suction conditions that become even more pronounced with successive increases in $M_{0}$ or $\unicode[STIX]{x1D6FE}$ . At this point, it may be useful to recall that, for cyclonic motions to be stable, the upward streaming of the incoming fluid through a siphoning process is necessary to avoid premature short-circuiting or early spillage, as it were, out of the open base. It can hence be seen that suction conditions near the centreline can be beneficial to the proper and stable formation of a bidirectional vortex. In consequence, it may be ascertained that increasing the injection Mach number or the ratio of specific heats will enhance the suction level in the core region, a condition that can lead to a more stable cyclonic-flow field.

Figure 12. Pressure distribution for injection Mach numbers of (a) $M_{0}=0.1$ and (b) 0.2 using $\unicode[STIX]{x1D705}=0.1$ and 0.3 at $z=l$ .

6.5 Sensitivity to the inflow swirl parameter $\unicode[STIX]{x1D705}$

Up to this point, most solutions have been evaluated at values of $\unicode[STIX]{x1D705}$ that range from 0.1 to 0.4. This convention has enabled us to amplify the effects of compressibility to the extent of better isolating and capturing the specific features associated with each variable of interest. Realistically speaking, it is possible for $\unicode[STIX]{x1D705}$ to take on smaller values, and these will lead to a reduction in the compressible axial and radial speeds along with their compressible counterparts relative to the tangential velocity. From this perspective, the sensitivity of the compressible approximation to variations in $\unicode[STIX]{x1D705}$ may be helpful to explore.

Figure 13. Sensitivity of the axial velocity in the exit plane for injection Mach numbers of (a) $M_{0}=0.1$ and (b) 0.2 using $\unicode[STIX]{x1D705}=0.1$ , 0.2, 0.3 and 0.4 at $z=l$ .

To study the $\unicode[STIX]{x1D705}$ sensitivity, the axial velocity profile is re-examined at $z=l$ using both $M_{0}=0.1$ and 0.2. This is accomplished over a range of $\unicode[STIX]{x1D705}$ extending from 0.1 to 0.4 as depicted in figure 13. It may be safely argued that the remaining dynamic and thermodynamic quantities will exhibit similar trends by virtue of their sensitivity to the swirl parameter $\unicode[STIX]{x1D705}$ being analogous to that of the axial velocity. As clearly illustrated in these graphs, decreasing $\unicode[STIX]{x1D705}$ leads to a corresponding drop in both compressible and incompressible parts of the axial velocity. The compressible contributions diminish even more rapidly, owing to their cubic dependence on $\unicode[STIX]{x1D705}$ , to the extent of approaching the incompressible approximation. Conversely, increasing the injection Mach number to 0.2 or higher stands to offset the effect of decreasing $\unicode[STIX]{x1D705}$ .

At low values of $\unicode[STIX]{x1D705}$ , the axial and radial velocities, which can directly absorb the effects of compression in the absence of restrictions in the $r$ and $z$ directions, become overwhelmingly dominated by the tangential motion. Their overall magnitudes become small relative to $v$ . The latter cannot experience compression in the tangential direction without bypassing the condition of axisymmetry. Its sensitivity to density variations can only be realized through its spatial dependence on the first-order streamfunction, and this association proves to be commensurate with the size of both $\unicode[STIX]{x1D705}$ and $M_{0}$ . Naturally, this coupling weakens at decreasing values of $\unicode[STIX]{x1D705}$ which, physically, imply the existence of higher levels of swirl and, therefore, stronger tendency to promote an axisymmetric distribution of flow field properties. So while higher levels of swirl increase the resistance to compression in the tangential direction, higher injection Mach numbers serve to counterbalance this effect, with the overall motion being controlled almost exclusively by these two contending factors. This behaviour seems to support the tradition of relying on incompressible models for mean flow Mach numbers below 0.3 irrespective of the flow detail.

6.6 Vorticity

Given the helical nature of this motion, it may be instructive to examine the influence of compressibility on the flow vorticity. In this vein, the radial, tangential and axial vorticity components, along with their total vorticity magnitude are plotted in figures 14(a)–14(d), respectively. Of the component vorticities, the axial and tangential contributions appear to be comparable in strength, while the radial component persistently lags behind by one order of magnitude. The swirl velocity, being the dominant component, seems to be the primary contributor to the axial vorticity, by virtue of its strong radial dependence. Similarly, the axial velocity, being the next largest contributor, forms the bulk of the tangential vorticity component. Since the swirl velocity exhibits only a weak dependence on the axial coordinate, its contribution to the radial vorticity is seen to be minimal compared to the other components.

Figure 14. The effect of compressibility on the vorticity in the (a) radial, (b) tangential, (c) axial direction and (d) vorticity magnitude for $\unicode[STIX]{x1D705}=0.1$ , 0.2 and 0.3.

Upon closer examination of figure 14, it may be inferred that re-introducing compressibility tends to shift the vorticity in the flow outwardly, towards the walls. This behaviour seems to be characteristic of all spatial components whose peak values are shifted towards the wall with successive increases in $\unicode[STIX]{x1D705}$ . Interestingly both radial and tangential vorticity components may be viewed as being anchored at both the centreline and the wall. The axial vorticity displays no such anchor, but rather assumes a positive, finite value near the centreline and a negative, finite value near the wall. Furthermore, when accounting for dilatational effects, the peak vorticity at the centreline may be seen to diminish as the vorticity redistributes itself near the sidewall.

To amplify the compressible correction, one may increase either $\unicode[STIX]{x1D705}$ or the Mach number. At $\unicode[STIX]{x1D705}=0.1$ , the Mach number does not seem to have an appreciable bearing on the tangential, axial or total vorticity profiles. Although it displays a visible effect on the radial vorticity, the magnitude of this contribution remains much smaller than the other components to the extent of producing a negligible impact on the total vorticity.

However, as the value of $\unicode[STIX]{x1D705}$ is increased, the influence of the Mach number becomes gradually more noticeable. For the $\unicode[STIX]{x1D705}=0.2$ case, a distinct difference arises in the $M_{0}=0.2$ plot compared to both the incompressible and $M_{0}=0.1$ findings. This effect becomes even more pronounced at $\unicode[STIX]{x1D705}=0.3$ , where deviations become magnified in comparison to the smaller $\unicode[STIX]{x1D705}$ values.

At his juncture, it may be useful to digress and discuss the practical range of $\unicode[STIX]{x1D705}$ based on laboratory experiments. In actuality, most cold-flow studies, such as those by Rom (Reference Rom2006), seem to rely on modular vortex chambers with values of $\unicode[STIX]{x1D705}$ that range from approximately 0.007–0.04. While it is possible for a design to exhibit values of $\unicode[STIX]{x1D705}$ that are as large as 0.3–0.4, laboratory experiments that employ water or air as the working fluid will likely operate at smaller values of $\unicode[STIX]{x1D705}$ . It may be hence speculated that compressibility will not likely have an appreciable effect on the vorticity field in non-reactive flow configurations.

6.7 Compressibility criterion

Before attempting to model a given flow field, it is often desirable to determine whether or not the extra effort entailed in pursuing a compressible solution will be warranted. This argument may be made both for analytical models, such as the one developed here, and for numerical simulations, where the inclusion of dilatational effects can significantly prolong the ensuing computation. Bearing this concern in mind, one may conceive of a compressibility criterion that captures the effective change in velocity due to density variations. Algebraically, this criterion may be taken as

(6.2) $$\begin{eqnarray}\unicode[STIX]{x1D712}=\frac{u^{2}+v^{2}+w^{2}}{u_{0}^{2}+v_{0}^{2}+w_{0}^{2}},\end{eqnarray}$$

where $\unicode[STIX]{x1D712}$ represents the change in the compressible velocity relative to its incompressible form. The resulting comparison of velocity magnitudes provides a single value that can be used to quantify the effects of fluid compression on the flow field. For example, one may consider that a ten per cent (or larger) deviation between compressible and incompressible velocities will justify the use of a fully compressible model. To illustrate this variability in the case of a bidirectional vortex, the values of $\unicode[STIX]{x1D712}$ in the $r$ – $z$ plane are computed and presented in figure 15. For the case of $M_{0}=0.1$ and $\unicode[STIX]{x1D705}$ ranging from 0.1–0.4, the compressible ratio may be seen to always fall within approximately five per cent of the incompressible approximation. Although the shape of the contours visibly change when $\unicode[STIX]{x1D705}$ reaches 0.4, the values of $\unicode[STIX]{x1D712}$ still fluctuate between five per cent below to five per cent above the incompressible prediction. We conclude that, at this level of the injection Mach number, numerical simulations could be safely carried out assuming incompressible conditions.

Figure 15. Contours of the compressible criterion for $M_{0}=0.1$ with (a) $\unicode[STIX]{x1D705}=0.1$ and (b) 0.4, as well as $M_{0}=0.2$ with (c) $\unicode[STIX]{x1D705}=0.1$ and (d) 0.4.

Interestingly, in reference to figures 15(a) and 15(c), it appears that when $\unicode[STIX]{x1D705}=0.1$ , which corresponds to a tangential velocity that is one order of magnitude larger than its axial and radial counterparts, a unidirectional velocity modification is realized across the chamber width, with the compressible velocity mirroring the incompressible profile near the exit plane, albeit with increasing divergence in the vicinity of the headwall. When $\unicode[STIX]{x1D705}$ is incremented to 0.4, the iso-contour of unity (along which the compressible and incompressible models coincide), is no longer located near the exit plane. Instead, it realigns itself almost vertically, thus bisecting the chamber into two regions: a core region where the compressible solution exhibits a reduction in kinetic energy, and a larger, outer region extending from approximately one third of the radius to the wall, where the compressible velocity exceeds its incompressible counterpart. The existence of two regions with dissimilar kinetic energies may be traced back to the mass conservation requirement and its enforcement within the asymptotic analysis. Since mass must be conserved at the leading order, the mass balance prescribed by (4.4) must also be secured at the first order. This implies that any velocity augmentation near the open end must be offset by a velocity deficit in order to maintain a balance between the inflow and the outflow. The clearest example of this behaviour may be inferred from figure 15(b). In this case, the compressible motion near the centreline undergoes a reduction in magnitude, while the flow adjacent to the walls is slightly accelerated. Inside the outflow region, which is delineated by the mantle interface, i.e. $0<r<0.628$ , the lower velocities that occur near the centreline must be offset by appropriately higher velocities within the outer portion of the outflow region. The outcome may be seen in a steepened velocity profile, where the core velocities are reduced as $r\rightarrow 0$ , while the outlying velocities are increased to preserve mass conservation throughout the exit port.

Naturally, the variability from the incompressible baseline is accentuated with successive increases in the Mach number. In figure 15(c,d) the magnification effect that accompanies a larger value of the Mach number may be clearly seen. Overall, in comparing the iso-lines in figure 15, one can conclude that while the shape of the contours is prescribed by the strength of $\unicode[STIX]{x1D705}$ , the size of the deviation depends on the value of the injection Mach number. By raising the injection Mach number to 0.2, the contour shapes remain similar, but the deviations at the extremes are amplified to approximately twenty five per cent, thus marking an order of magnitude increase relative to the $M_{0}=0.1$ case. The inclusion of compressibility effects will hence be essential in bidirectional vortex applications, such as the VCCWC studies, which entail injection Mach numbers that exceed 0.2.