Effects of nozzle inlet size and curvature on the flow development in a bidirectional vortex chamber

The characterization of wall-bounded cyclonic flows continues to receive attention in the propulsion and fluid dynamics communities, and this is largely due to the favorable attributes associated with the onset of rotating motion in vortex engines, cyclone separators, and other swirl-dominated chambers.¹ By translating the basic characteristics of unconfined tornadoes and hurricanes to propulsive, combusting, and flow separation devices,² one is able to exploit the high energy density, self-sustaining capability, self-cooling properties, prolonged residence time, flow separation efficiency, and effective mixing potential of cyclonic flows to optimize the design space, enhance the controllability, and improve the stability of numerous power production and thrust generation devices.³ 

Depending on the number of end-to-end flow passes that a fluid entering an enclosure will undertake before exiting, confined cyclonic motions can be classified as unidirectional, bidirectional, or multidirectional.⁴ The primary emphasis of the present investigation will be on the bidirectional cyclonic pattern, which contains two distinct swirling regions: an outer annular vortex, where the fluid spirals toward the headwall, and an inner vortex region, where the flow targets the endwall, specifically where nozzles, outlet tubes, and vortex finders can be mounted.⁵ Furthermore, the two characteristic segments of the flow are separated by a virtual interface consisting of a rotating but non-translating layer termed “mantle.” In a cyclonic chamber, swirl is sustained in accordance with the conservation of angular momentum of a tangentially injected fluid, often introduced at one end of the chamber (usually the endwall), and extracted at the same location after completing a bidirectional pass.

The earliest studies of bidirectional swirling flows may be traced back to ter Linden⁶ and his experimental investigations of the geometric dimensions on the efficiency of cyclonic dust collectors. Then after four decades of sporadic studies, Bloor and Ingham⁷ develop a simple incompressible model for wall-bounded conical cyclones from which a useful closed-form solution is retrieved. This is followed by a series of numerical simulations of confined cyclonic motions by Hsieh and Rajamani,⁸ Monredon, Hsieh, and Rajamani⁹ and Hoekstra, Derksen, and van den Akker.¹⁰ In addition to their engagement in numerical investigations, some of these researchers perform experiments to validate their computational results. Later, by redirecting attention back to the development of a viable mathematical framework for wall-bounded cylindrical cyclones, Vyas and Majdalani⁵ and Majdalani and Rienstra¹¹ manage to produce exact inviscid solutions for the bulk cyclonic motion starting from first principles. However, to overcome the singularity that undermines the inviscid motion along its axis of rotation, Majdalani and Chiaverini¹² apply the tools of boundary-layer theory to obtain viscous corrections that capture the forced vortex behavior in the core region of a bidirectional chamber. Along similar lines, Maicke and Majdalani¹³ are able to overcome the inviscid singularity that evolves at the centerline by implementing a constant shear stress model. The necessity of securing the velocity adherence requirement near the sidewall is also investigated using matched-asymptotic expansions in a series of studies that lead to viscous approximations for the core and wall boundary layers.^12,14 Shortly thereafter, several distinct families of helical motions are introduced and shown to exhibit rather compelling cyclonic flow patterns that correspond to either “Beltramian” or “Trkalian” fields.³ Other complementary studies are pursued for the purpose of incorporating the effects of arbitrary headwall injection,^15,16 compressibility,^17,18 uniform sidewall injection,^19,20 and hollow cylindrical air cores.²¹ The development of exact cyclonic solutions for hemispherically bounded motions is also achieved under both irrotational²² and rotational conditions, with the latter leading to essentially one complex-lamellar²³ and two Beltramian mean flow profiles.²⁴ 

Apart from the filtration industry, where dedusters and hydrocyclones are commonplace, several thrust chambers in the propulsion industry have been conceived to operate optimally by leveraging the bidirectionality and swirl advantages of cyclones. Of those, two particular concepts come to mind, namely, the vortex injection hybrid rocket engine (VIHRE), developed by Knuth et al.,^25,26 and the vortex combustion cold-wall chamber (VCCWC), conceived by Chiaverini et al.²⁷ For the reader's convenience, schematics of the VIHRE and VCCWC concepts are furnished in Figs. 1(a) and 1(b). The VCCWC, which is more relevant to this article's underlying theme, consists of a unique design concept wherein the oxidizer is injected tangentially at the endwall of the chamber just above the nozzle. After entering the chamber rather unconventionally, the oxidizer ascends along the walls in a spiraling motion that forms the “outer vortex” region. As it approaches the headwall, the oxidizer then turns inwardly, mixes with the fuel, and combusts vigorously, even as it courses its way down the bore, where it is referred to as the “inner vortex” core. Through this helical trajectory, the outer annular vortex confines the reactants to the core region and prevents the hot products and gaseous mixtures from coming in contact with the walls. In practice, the resulting film cooling mechanism is shown to be very effective at protecting the walls from the high convective and radiative core temperatures.¹ 

It can be inferred from a perfunctory literature review that most of the experimental studies related to confined cyclonic motions have been carried out on cyclone separators or dust collectors. As such, most of the information pertaining to bidirectional vortex engines remains either limited or proprietary; few cases can be found including the cold-flow experimental and computational studies performed by Anderson et al.²⁹ and Rom, Anderson, and Chiaverini.³⁰ Their work considers a laboratory-scale model of a simplified cyclonic chamber and produces tangential velocity data at three specific axial locations using particle image velocimetry (PIV). These PIV data are later used by Maicke and Majdalani^13,31 to guide their viscous and constant shear stress formulations.

The use of numerical models in the treatment of cyclonic flows has continued to receive attention due to the evolution of hardware and affordability of resources relative to experimentation. In this vein, numerical simulations that describe wall-bounded cyclones have continued to appear and these have attempted different turbulence models in the hope of capturing the anisotropic behavior of the flow in the presence of a highly swirling flow environment. Despite these pursuits, investigations of precession, flow instability, and breakdown, which are susceptible to numerical simulations, have remained largely unattempted and unresolved. Moreover, the identification of appropriate turbulence models that are capable of handling highly anisotropic conditions continues to pose open challenges. In what follows, some of the pertinent efforts are briefly reviewed.

We begin with the design optimization study of Boysan, Ayers, and Swithenbank,³² where a two-phase algebraic turbulence model is constructed for gaseous cyclones with entrained solid particles. Their model computes six Reynolds stress components, thus leading to a numerical scheme that is capable of predicting grade-efficiency curves by means of stochastic particle tracking techniques. Nearly a decade later, another detailed study is conducted by Concha et al.,³³ where different sets of equations are applied to each of six individual zones constituting a hydrocyclone. These equations are intended to permit the computation of local velocities with considerable precision. Then using the momentum balance of particle trajectories, the characteristic separation size is accurately deduced.

Along similar lines, Hoekstra, Derksen, and van den Akker¹⁰ carry out a comprehensive numerical investigation of cyclone separators using air as the working fluid. Their study provides a comparison between different turbulence models and experimental results. Accordingly, the Reynolds stress transport model is found capable of producing more realistic predictions relative to other turbulence models. A year later, Derksen and van den Akker³⁴ extend this line of inquiry by using large eddy simulations (LES) to predict the fundamental features of a reverse-flow cyclone separator. Their study also discusses the precession patterns of the inner vortex and its stabilizing role. Advancing further in their investigation, Derksen³⁵ applies this LES code to predict the separation performance of a high efficiency Stairmand-type cyclone using frozen-field, eddy-lifetime, and periodic-flow approaches. Derksen, van den Akker, and Sundaresan³⁶ also examine the effects of incorporating two-phase flow effects on the fidelity of their computations by successfully implementing two-way coupling. In the meantime, Hu et al.³⁷ modify the empirical constants for isotropization of their production and convection model (IPCM) as well as the wall-reflection term of their pressure-strain model to better handle anisotropic turbulence. This enables them to produce an improved Reynolds stress model (RSM) for strongly swirling flows while dodging the computational costs associated with LES. By way of confirmation, their three-dimensional computations are shown to agree with measurements acquired in a cyclone separator using laser Doppler velocimetry (LDV). In fact, this sequence of studies builds up to a comprehensive review by Cortes and Gil,³⁸ where different computational models for reverse-flow cyclone separators are surveyed and cataloged based on their ability to predict the tangential velocity distribution, pressure drop, and collection efficiency.

Despite the existence of several numerical studies of cyclone separators, only a few may be found that devote themselves to the VCCWC flow configuration using, for instance, the laminar, k–ε, and Reynolds stress transport (RSM) models. Among those, one may cite Fang et al.,³⁹ Anderson et al.,²⁹ Murray et al.,⁴⁰ Rom, Anderson, and Chiaverini,³⁰ Hu et al.,³⁷ Maicke and Talamantes,⁴¹ Talamantes and Maicke,⁴² Majdalani and Chiaverini,¹ Sharma and Majdalani,^43,44 and Rajesh, Jothi, and Jayachandran.⁴⁵ In this work, we augment this line of inquiry by taking a close look at the cyclonic motion in an idealized VCCWC chamber using a robust finite-volume solver. Since most previous studies have focused on describing the tangential velocity character, our emphasis will be broadened to perform cold-flow simulations that help to characterize all three components of the velocity and pressure fields while exploring the effects of the chamber outlet radius and nozzle extension size on the positioning and formation of mantles and backflow regions. Moreover, unlike a previous study that has focused on the use of pressure boundary conditions,⁴⁴ the present study will implement velocity boundary conditions that enable us to achieve substantially more realistic flow speeds throughout the chamber.

In previous work,⁴³ the main objective was to compare numerical results to those of an exact complex-lamellar solution derived by Vyas and Majdalani⁵ directly from Euler's equation. To do so, the geometry and boundary conditions in the numerical model were chosen as plainly as possible to mirror those associated with the theoretical framework. In view of the purely tangential injection assumed in the complex-lamellar model, special adjustments had to be made to align the projected injectors in a similar arrangement along the outer circumference with no outward protrusions.⁴³ 

In the present work, we opt for a more realistic model, as depicted in Fig. 2, in order to reproduce a cyclonic motion that mimics the VCCWC bulk flowfield. In this idealization, the focus will be on the fundamental physics of the wall-bounded cyclonic motion under incompressible conditions; it will be carried out in the absence of chemical reactions or gaseous expansions that typically take place through converging-diverging De Laval nozzles. In this vein, the model will not be used to predict thrust performance characteristics. Rather, it will be relied upon to identify the characteristic cyclonic flow features and compare them to those predicted by existing analytical solutions under inviscid conditions.^3–5 To this end, a three-dimensional computational domain is conceived for an idealized depiction of the bidirectional vortex engine with eight protruding injectors. Note that, in Fig. 2(b), the axis of the chamber is made to coincide with the z-axis, which is measured downward starting at the headwall. Our computational domain is protracted using SolidWorks with dimensions that are identical to those used by Anderson et al.²⁹ It is subsequently meshed using a well-established meshing software called Integrated Computer Engineering and Manufacturing (ICEM CFD). The discretized domain is then exported into our finite-volume solver, which is part of a commercially available multiphysics package incorporating an implicit, incompressible, Reynolds-averaged Navier–Stokes solver.⁴⁶ The geometric details of the chamber and its injectors are provided in Table I and illustrated in Fig. 2. In the interest of clarity, Fig. 2(a) provides an isometric view of the idealized vortex chamber, whereas Fig. 2(b) shows a side view that displays the actual dimensions used in the upcoming simulations.

Creating a complete and effective structured mesh remains, perhaps, one of the most labor-intensive components of this study. The difficulty may be ascribed to the complications brought about by the need to carefully resolve and refine the mesh in the region surrounding the tangential injectors as well as the sidewall. For example, the small angle between the injectors and the chamber wall causes the mesh elements in this region to become highly skewed. To overcome this issue, a mesh consisting of tetrahedral elements is generated in a manner to minimize the skewness. The mesher that we use helps to evaluate the normalized mesh quality by reporting zero for a highly skewed element and unity for a perfect, unskewed element. Bearing these factors in mind while creating the mesh, the process is repeated until the normalized quality of the mesh is seen to return a value above the minimum recommended threshold of 0.3.

To make further headway, we recognize that a tetrahedral mesh requires a relatively large number of elements for a given volume, thus increasing both the computation time and power requirements. As is customary, our tetrahedral mesh elements are converted into polyhedral elements using the reverse Cuthill–McKee algorithm.⁴⁷ This step reduces computation time while improving convergence properties; in fact, it leads to simulation outcomes that are characteristic of a fully structured hexahedral mesh. In what follows, Fig. 3(a) provides a magnified view of the meshed geometry near the injector and sidewall regions, whereas Fig. 3(b) displays the geometric mesh over the entire chamber and injector openings.

After mesh generation and the specification of boundary types, the mesh is imported into a robust finite volume solver. Simulations are then carried out using a steady, implicit, 3D, segregated solver. The computations are performed using a strictly inviscid model that admits slip-permitting walls. This is done for the purpose of facilitating the capture of stable and coherent cyclonic patterns that correspond to the available analytical solutions.

Since one of the objectives of this study is to compare our simulation results to existing analytical solutions, we find it helpful to normalize all variables using the same reference values adopted in the corresponding theoretical formulations.^3–5 Accordingly, one may use overbars to denote dimensional variables and write

In the above, $r¯$ and $z¯$ refer to the radial and axial coordinates, $Uw(z)=uθ(1,z)$ represents the swirl velocity sweeping the sidewall, and p_w stands for the pressure at the sidewall.

The conventional equations for conservation of mass and momentum are adapted for the flow in a cyclonic chamber. Using ρ and $u$ to represent the density of the fluid and its velocity, the fundamental equations that we solve consist of

where p stands for the static pressure, $ρg$ symbolizes the gravitational body force, and $F$ accounts for other external forces, such as the chamber's reaction forces exerted on the fluid.⁴⁶ In this work, air is chosen as the working fluid with a density of $ρ≈$ 1.225 kg/m³. The improved mesh and chamber configuration will be shown to facilitate the characterization of all three components of the velocity, including the axial and radial components, which were omitted in two previous studies.^42,43 They will also enable us to explore the effect of changing the nozzle inlet diameter and radius of curvature on the cyclonic flow properties. To this end, the normalized outlet radius, β, will be varied while keeping all other parameters fixed. The normalized outlet radius, β, is simply the ratio $β≡b/a$⁠, where b defines the chamber outlet radius, also denoting the nozzle inlet radius, while a represents the chamber radius.

In our parametric trade analysis study, which is summarized in Table II, the values of β are selected in a manner to coincide with their theoretical values for a cyclonic motion with multiple mantles, as provided for the complex-lamellar family of solutions.⁵ The same simulations can be as easily repeated using the multiple mantle locations associated with the Beltramian family of solutions that is available in the literature.³ The baseline geometry is chosen to match the one used by Anderson et al.²⁹ Since the flow remains inviscid and incompressible, a uniform top-hat velocity profile of 60 m/s is specified as a plug flow at each injector inlet tube, which extends 5 cm tangentially to the chamber. As for the outlet, it is taken to be at atmospheric pressure. With the inlet being set as a velocity boundary condition, the robust SIMPLE algorithm⁴⁷ is then used to couple the pressure with the velocity. In this process, the pressure is discretized according to the so-called “standard” scheme and the momentum equation is discretized using a first-order upwind stencil. Later during execution, as the residuals start to stabilize, and to achieve higher order accuracy, the pressure is rediscretized using a second-order scheme hand-in-hand with the momentum equation, which is recast using a second-order upwind technique. This is performed, in part, to prevent the sudden divergence of residuals due to the complexities that accompany cyclonic flows. Besides constructing scalar values at cell faces, gradients are also implemented to compute secondary diffusion terms and velocity derivatives. Next, the convection and diffusion terms in the flow conservation equations are discretized using the gradients of variables. More specifically, the gradients are computed using the Green–Gauss cell-based scheme.⁴⁷ Finally, to promote a reasonable degree of precision, stringent convergence criteria are set. These ensure that the residuals remain consistently below a $10−6$ tolerance level in all three velocities, as well as the continuity equation. In addition to this convergence criterion, volume integrals of the tangential, axial, and radial velocities are computed and monitored to ensure proper convergence.

The primary aim of performing a numerical investigation of the bulk gaseous motion in a bidirectional vortex engine is to better understand the flow physics and identify its controlling parameters. In this vein, we proceed with a characterization of the axial and radial velocity components to complement previous results that mainly focused on the tangential component.⁴³ In addition, a parametric study is conducted to evaluate the effect of varying the outlet diameter and radius of curvature on the flowfield. As indicated in Eq. (1), all velocities in this study are normalized by the tangential velocity at the wall, $Uw(z)=uθ(1,z)$⁠.

In order to more effectively relate our computational results to existing analytical solutions, effort has been made to ensure that the CFD settings match most of the parameters that are used in the analytical model, including, for example, the modified swirl number

In the above, a denotes the radius of the chamber and $Ainj$ refers to the total injection area. The second characteristic parameter is often dubbed the off-swirl parameter;³ it appears naturally in the analytical solutions and serves to gauge the magnitude of the axial and radial velocities relative to their tangential counterpart. It is given by

where L represents the chamber length. Note that a lower value of κ implies larger tangential-to-radial or tangential-to-axial velocity ratios. In what follows, illustrative simulation results will be presented and these will correspond to our baseline Case 3 in Table II with $β=0.5$ and, in view of Table I, $κ=1/(2π×2.805)×(2.54/6.1)=0.0236$⁠, unless specified otherwise.

We begin with the tangential velocity distribution in Fig. 4(a), where the radial variation of $uθ$ is given at several discrete axial stations, starting at z/L = 0.1 near the headwall and ending with a value of 0.8 near the outlet entrance region. To be consistent with theory and tradition, such as Rankine's free vortex plots, the tangential velocity is normalized using the tangential speed, $Uw(z)=uθ(1,z)$⁠, taken at the respective z/L locations. These results are shown side-by-side with the corresponding swirl velocity contour plots in Fig. 4(b). Despite the elongated scale used in Fig. 4(a), it may be seen that the axial variations of the swirl velocity are minimal. Unsurprisingly, the maximum swirl velocity is slightly higher near the injection site at $z/L=0.8$ and decreases as the headwall is approached. This can be more readily observed in the iso-velocity contours, which appear to be rather uniform along the length of the chamber, with a very minor decrease in the width of the $(uθ)max$ isoband as z is increased. The corresponding behavior is reassuring as it helps to confirm the axial invariance assumption of the swirl velocity, which is often used at the basis of deriving cyclonic flow solutions.³ Furthermore, the existence of a forced vortex motion in the core region, followed by a free vortex tail beyond the point of maximum $uθ$⁠, can be ascertained from the swirl velocity plot given in Fig. 4(a). Note that no values past $z/L=0.8$ are shown to avoid naturally-occurring nozzle entrance effects.

Next, we proceed to examine the pressure distribution in Fig. 5, where it is normalized by its maximum radial value taken at r = 1. Although all cases adopt the same inlet pressure at their injectors entry point, this maximum pressure value at steady state will slightly differ depending on the axial station within the chamber. Specifically, as we travel within the outer annular vortex, the normalizing pressure value taken at the sidewall will slightly diminish as the headwall is approached. Here too, one may infer that the pressure gradient along the radius of the chamber remains nearly independent of the axial distance, particularly in Fig. 5(a), although minor deviations seem to occur near the axis of rotation. Naturally, if we consider the outer vortex region corresponding to $r≥0.5$⁠, the pressure near the headwall must be slightly lower than the pressure near the injection site to force the upward convection of the incoming fluid. In contrast, the pressure within the inner vortex will slightly diminish in the direction of the outflow, namely, as z/L increases from 0.1 to 0.8. Nonetheless, these variations appear to be minimal at steady state if we are to travel axially along a constant r line. A similar trend may be inferred from the pressure contours of Fig. 5(b), where the axial invariance of the pressure is fundamentally established across the entire domain. In fact, a similar pressure distribution is reported by Hoekstra, Derksen, and van den Akker,¹⁰ Hu et al.,³⁷ Murray et al.,⁴⁰ and Majdalani and Chiaverini.¹ In short, the radial values of the pressure are found to be largest at the sidewall and smallest along the chamber axis; axially, however, the sidewall pressure is slightly lower at the headwall than the endwall, which is a requirement to properly siphon the incoming flow upwardly; the converse is true for the core pressure, which slightly decreases in the outflow direction.

As we turn our attention back to the axial velocity, its behavior is illustrated in Fig. 6(a), where the radial distribution of u_z, normalized using the corresponding tangential speed at the wall, $Uw(z)$⁠, is shown at eight different values of z/L. As usual, we start with a chained line at $z/L=0.1$ and stop with a broken blue line at a value of 0.8 near the outlet. The companion Fig. 6(b) displays the axial velocity contours across the chamber. It is clear from the plot that there exist two distinct regions that comprise an inner core and an outer, annular segment. Motion in the outer segment, which is positive upward until roughly $r≈β=0.5$⁠, turns negative in the core region as the flow polarity switches to a downward motion.

The so-called updraft and downdraft regions, with two opposing axial polarities, are the reason a cyclonic motion is often dubbed “bidirectional” or “bipolar.” These flow segments are separated by a rotating but non-translating sheet or mantle, which is depicted using a solid line in Fig. 6(b). This interface corresponds to the locus of u_z = 0, which evolves around a radius of $r≈0.504$⁠, fairly closely to the nozzle entrance radius β for this configuration. It can also be confirmed in Fig. 6 that the magnitude of $|uz|$ increases in the core region as the fluid attempts to leave the chamber. This behavior is further corroborated in the axial velocity contours of Fig. 11(a), where the dark blue color denotes incoming fluid moving upwardly, while the light blue color marks the outgoing core flow.

Finally, the normalized radial velocity distribution in the chamber is explored in Fig. 7(a), where symbols are used to represent different axial stations and thus help to illustrate the statistical scatter of the solution around its average value. Here, it can be seen that the average radial velocity is qualitatively similar to the theoretical form of u_r predicted by the family of analytical solutions provided in the literature.³ This includes the presence of a maximum absolute value for $|ur|$ midway along the radius. Note that the solid line corresponds to a polynomial fit for the scatterplot and, as such, constitutes a rough representation of the average radial velocity obtained for this baseline case. The complementary radial velocity contours in Fig. 7(b) confirm these trends across the whole chamber where the axial invariance of u_r may be observed in the upper three quarters of the domain. In the bottom quarter length, the bulging and contraction of the contours may be attributed to the presence of finite injectors and a nozzle attachment in their immediate vicinity, where a strong flow turning activity is prominent. The axial uniformity of these contours and their confirmatory curves in Fig. 7(a) lend support to the axial invariance of the radial velocity, an assumption that is often made in theoretical investigations of this problem.^3–5 This behavior may be attributed to the imposition of slip walls everywhere in the present simulations and to the sufficient distance from the injectors.

Having showcased the capabilities of the numerically computed flowfield at predicting the axial and radial velocities, we now proceed to vary the outlet diameter to study its effects on the axial velocity, the mantle location, the swirl velocity, and the pressure distribution. Naturally, with a fixed inlet velocity condition, varying the outlet diameter directly affects the exit velocity. The next parametric study is therefore performed to better understand the influence of the normalized outlet radius, β, on the solution variables. Having already explored the flow characteristics using a benchmark $β=0.5$⁠, this value will serve as the control case for the remaining tests. As we proceed, the normalized outlet radius is varied between 0.354 and 0.935, as shown in Figs. 8–10, where the tangential velocity, pressure distribution, and axial velocity are examined individually. We recall that the specific values of β are chosen because they correspond to the first theoretical investigation of multiple mantles carried out by Vyas et al.,⁴⁸ namely, where the locus of multiple mantles is identified and subsequently confirmed to a certain degree by the experimental measurements and CFD simulations of Anderson et al.²⁹ and Rom, Anderson, and Chiaverini.³⁰ 

Forthwith, Fig. 8 illustrates the behavior of the inviscid form of the tangential velocity, normalized by $Uw(z)$⁠, for various outlet diameters and a fixed station of $z/L=0.5$⁠. As expected, successive increments in the normalized outlet radius serve to decrease the peak swirl velocity in the chamber such that $(uθ)max$ is largest for $β=0.354$⁠. Rather interestingly, for $β>0.707$⁠, the peak swirl velocity moves consistently outwardly closer to the sidewall in an attempt to compensate for an expanding outlet area. A similar pattern is reported by Talamantes and Maicke⁴² in their confirmatory computational results. The magnitude of the normalized $uθ$ near the sidewall also approaches a value of unity with successive increases in β, another trait that may be graphically inferred from Fig. 8.

Next, a pressure comparison is featured in Fig. 9 for the different cases of β and $z/L=0.5$⁠. As before, the pressure is normalized by its value taken at the wall. The first observation that we make in Fig. 9(a) is that the net pressure drop at the centerline increases as the outlet area is enlarged. Although the wall pressure is not shown, we find that, in all cases considered, the lowest pressure occurs at the centerline and the largest sidewall pressure is highest for the case with the smallest outlet area, or lowest β. Moreover, for $β>0.707$ in Fig. 9(b), the pressure near the centerline seems to flatten: the pressure gradient decreases, the peak pressure drops, and the pressure in the inner core plateaus to an essentially constant value. This behavior can be attributed to the presence of a central recirculation zone (CRZ) that becomes more appreciable as the outlet area is enlarged.^49–54 The reason for recirculation to occur will become apparent in the discussion below and will be associated with the emergence of a backflow region for $β≥0.707$⁠. Except for the backflow effects, these trends seem to be comparable to those predicted theoretically by the Trkalian profile with no slip.³ 

To further verify the trends exhibited by the numerical simulations, $uθ$ is compared to the tangential velocity profiles obtained using the constant core shear model.¹³ As the name implies, the latter assumes a constant shear stress to be dominant in the forced vortex region, which typically extends between the chamber axis and $r=rmax$⁠, where $uθ=(uθ)max$⁠. It then transitions over the $rmax≤r≤X$ interval into a free vortex motion that becomes prevalent for r > X. This comparison is showcased in Figs. 8(c) and 8(d) for $β=0.628$ and 0.707, respectively. As one can gather from these plots, the profiles agree quite well with the constant core shear model for $uθ$⁠, namely,

where $X=rmaxe, Ue=12e(uθ)maxrmax≈1.359(uθ)maxrmax$⁠, and $rmax$ denotes the locus of $(uθ)max$⁠. In this comparison, we find the dimensionless values of $X≈0.75$ and 0.88 with $rmax=0.455$ and 0.535 for $β=0.628$ and 0.707, respectively. As for the normalizing speeds associated with our numerical simulations, we compute $Ue/Uw=0.728$ and 0.820 for the two cases in question.

Next, by switching to the axial velocity in Fig. 10, and its respective contours in Fig. 11, one can infer that the peak $|uz(0,z)|$ at the centerline decreases monotonically as the outlet area is increased. This behavior is consistent with continuity and concurrent with the theoretical finding that a Beltramian profile with $β=0.628$ will be accompanied by a higher axial velocity magnitude at the chamber outlet than the complex-lamellar profile with a (larger) $β=0.707$⁠. Such behavior is also corroborated in the axial velocity comparison presented by Majdalani⁵⁵ for the two types of profiles. Naturally, the magnitude of the axial velocity at the sidewall follows suit, and increases with β, as shown in Fig. 10. Note that for $β≥0.707$⁠, flow recirculation can be surmised from the axial velocity distribution in Fig. 10(b), and this is caused by the development of a backflow region that protrudes deeper into the chamber with successive expansions of the outlet diameter. More detail on the backflow region, and how to suppress it, is provided further below.

It can be gleaned from the polarized axial velocity contours in Fig. 11 that the mantle tends to adapt to the outlet size, particularly, by slightly repositioning itself to merge smoothly with the available outlet area at the bottom. In this case, the polarized values of u_z are tagged using a binary color scheme of either dark blue for upward or light blue for downward motions. For outlet radii greater than $β=0.707$⁠, the mantle turns outwardly, as the flow approaches the outlet, to accommodate the larger outlet size. In contrast, for $β<0.628$⁠, it gradually turns inwardly while merging with the nozzle inlet. For in-between values of $0.628≤β≤0.707$⁠, the diameter of the mantle remains nearly uniform as the flow exits the chamber.

It can also be remarked that as the exit radius is expanded in Fig. 11, a backflow region starts to form and move upstream into the chamber with successive increases in the outlet diameter. Specifically, the backflow remains secondary and confined within the nozzle section for $0.628≤β<0.707$⁠, i.e., with limited impact on the computed flowfield. However, for $β≥0.707$⁠, the backflow region becomes gradually more prominent in size as it protrudes deeper into the chamber with successive increases in β. The development of a backflow region may be ascribed to the negative pressures that can occur in the core region, as predicted by Fig. 9. In fact, similar behavior has been reported in other cyclonic flow studies including those by Smith,⁴⁹ Hall,⁵⁰ Faler and Leibovich,⁵¹ Zhang, Guo, and Yang,⁵² Talamantes and Maicke⁴² (for large values of β), Wang, Wang, and Yang,⁵³ and Khan.⁵⁴ Such backflow, however, is highly unlikely to evolve in the case of a choked compressible motion through a converging-diverging outlet nozzle, as in the actual VCCWC chamber.

Overall, the case with the largest $β=0.935$ is found to be the most unstable as observed in Fig. 11(l). Therein, the mantle tries to form but the outlet radius appears to be excessively large for the mantle to realign itself with the available outlet area at the bottom wall. Suppressing the backflow in this case, as in other such cases of internal backflow developments, requires geometric adjustments that must be judiciously investigated.

Before leaving this section, it may be instructive to note that the present simulations, which are based on fixed velocity inlet conditions, do not give rise to multiple mantles. Based on a previous investigation, we also find that the use of fixed pressure boundary conditions do not promote the development of multiple mantles either.⁴⁴ This leads us to speculate that the development of multiple mantles is not solely dependent on the outlet diameter, but rather influenced by other factors that remain unexplored. These may include the chamber aspect ratio, the nozzle diameter, the nozzle length, and the nozzle radius of curvature, as well as other operating conditions. In what follows, the effect of the nozzle length and radius of curvature are examined in more detail.

As we explore ways to prevent the backflow from forming, our initial step is to change the length of the nozzle extension. However, we find that increasing the nozzle length further can be problematic, as it can lead to unstable simulations. Conversely, reducing the nozzle length causes the backflow to penetrate deeper into the chamber, thus leading to larger CRZ segments. It thus appears that a finite range of nozzle lengths will produce a stable cyclonic flowfield for any $β≥0.628$⁠.

In lieu of the nozzle length, we find that increasing R_N, the radius of curvature of the converging section of the nozzle, can be quite effective at suppressing the undesirable backflow. After several iterations, a value of 0.84 cm is found to be suitable for $β=0.628$ and, similarly, a value of 1.04 cm is found to be adequate for $β=0.707$⁠. Recalling that the radius of the chamber is 2.54 cm, these curvatures translate into relative proportions of 0.33 and 0.41, respectively. Then noting that the normalized distances measured from the sidewall correspond to $(1−0.628)=0.372$ and $(1−0.707)=0.293$⁠, we realize that the effective radii of curvature for both Beltramian and complex-lamellar configurations correspond to the same constant value of $0.33×0.372=0.41×0.293=0.12$⁠. We can thus proceed to eliminate the backflow from all other cases using an effective radius of curvature that is deducible from a simple empirical correlation, namely,

The immediate effect of suppressing the backflow is illustrated in Fig. 12; therein, the simulations are repeated with a more optimal R_N, and used to generate the axial velocity contours along with their three-dimensional streamlines for the two baseline cases of $β=0.628$ and 0.707. On comparing Figs. 12(a) and 12(c) with Figs. 11(d) and 11(i), the elimination of recirculatory zones is evident. These contours also display quite distinctly the computationally-resolved mantle lines separating the outer (dark blue) and inner (light blue) regions. Although these mantles remain fairly uniform, which is consistent with theory, their numerically computed averages of 0.51 and 0.53 are slightly reduced due to the inward shifting of the mantle that is required to accommodate a rounded nozzle entrance.

These observations are reinforced in Fig. 13, where the radial distributions of the normalized u_z midway through the chamber are computed and displayed using our finite volume solver for the two nozzle configurations with and without backflow. The CFD results are also compared to available theoretical predictions that correspond in Fig. 13(a) to the Beltramian profile with $β=0.628$ and, in Fig. 13(b), to the complex-lamellar profile with $β=0.707$⁠. In the interest of clarity, the analytical velocity profiles are given by

where $J0(x)$ denotes the Bessel function and $λ1≈3.831 706$ represents the first zero of $J1(x)$⁠.

As before, both CFD and theoretical predictions are normalized using the tangential velocity at the sidewall, $Uw(z)$⁠. Based on the resulting graphs, one may note an overall agreement between theory and computation in the flow patterns, especially in comparing the theoretical profiles to their CFD predictions with no backflow. Almost paradoxically at first, the average estimation of the mantle location using CFD simulations with the original radius of curvature, i.e., with backflow development, leads to closer agreement with the theoretically predicted mantle locations. Upon closer look, this behavior may be explained by the larger deviation between the theoretical outflow configuration and the modified outlet shape resulting from further decreases in the nozzle's final diameter. More specifically, in the process of suppressing the backflow, an increased nozzle radius of curvature is imposed; this requires a smaller effective outlet diameter including a smaller final diameter. Comparatively, the effective diameter of the nozzle configuration that entails a backflow remains larger and, therefore, closer to the theoretical value of the targeted β. It can thus be realized that by increasing the radius of curvature of the outlet nozzle, the constricted outlet area causes the mantle location to shift radially inwardly relative to its theoretically predicted value. Moreover, although the R_N correlation is obtained empirically, its validity may be readily verified for various cases of β as shown in the contours depicted in Figs. 14(a)–14(d); therein, the presence of backflow regions is wholly eliminated. It should be noted, however, that for the case of $β≥0.760$⁠, the radius of curvature prescribed by R_N becomes so large that the converging nozzle walls end up crossing over the vertical z-axis, thus leading to an infeasible configuration. So long as this range is avoided, and by taking into account the complementary study that uses pressure inlet conditions,⁴⁴ the present investigation confirms that the R_N correlation is not only effective over a practical range of $β<0.866$ but also for both velocity and pressure inlet boundary conditions.

Finally, before concluding our flowfield characterization, we turn our attention to the behavior of the normalized radial velocity. As such, we compare the radial distribution of u_r in Figs. 15(a) and 15(b) for the two characteristic values of β, midway across the chamber, using the same CFD solutions that either allow or suppress the re-entrant backflow. Also included are the strictly inviscid analytical predictions based on the Beltramian and complex-lamellar models. For the reader's convenience, these are given by

Despite the similarities in the spatial patterns displayed by these solutions, the numerically computed radial velocities appear to be smaller in magnitude than their theoretical counterparts. Such differences may be attributed to the full inclusion of swirling effects in the CFD equations, unlike the analytical models, where the effects of swirl are essentially decoupled. In particular, the computational framework retains $uθ$ in the radial momentum equation; as such, the radial velocity remains computationally coupled to both the axial and tangential velocities in the attendant simulations, which have a strong bearing on the computation of u_r. On the contrary, $uθ$ is taken to be independent of u_r in the analytical models as a prerequisite for obtaining closed-form solutions. Consequently, the velocity profiles generated by the numerical model are seen to slightly differ from their theoretical counterparts.

Of the various computational studies that have undertaken the task of describing wall-bounded cyclonic flows, most have concentrated on industrial cyclone separators. A few have been specific to the bidirectional vortex motion associated with a vortex combustion cold-wall chamber (VCCWC) configuration. Furthermore, parametric studies related to outlet size variations remain either scarce or inaccessible. By imposing slip-wall conditions in concert with an inviscid flow solver, we show in this work that a judiciously constructed numerical model is capable of predicting the forced vortex behavior in the inner core region, as well as the free vortex motion in the outer region. These particular regions are well approximated analytically by a constant core shear stress model that entertains slip velocities along the walls and a constant shear stress around the core region. In a sense, it emulates to a certain degree the quasi-inviscid assumptions used by Maicke and Majdalani.¹³ The work also confirms that a numerical model is somewhat effective at characterizing the pressure distribution as well as other flow variables of interest. For example, the parametric trade analysis conducted so far suggests that increasing the normalized outlet radius, β, can have an appreciable impact on the swirl velocity, axial velocity, mantle location, and pressure distribution throughout the chamber. In fact, as the outlet area is increased, the peak tangential velocity is reduced and drawn radially outwardly. The pressure drop along the chamber axis is also increased until $β=0.707$⁠, beyond which the pressure curves start to asymptote in the inner core region. The same can be said of the mantle location, which seems to adapt itself to the expanding outlet diameter, until it can no longer attach itself to the nozzle inlet circumference. Incidentally, when the outlet diameter becomes sufficiently large, a backflow region seems to evolve in the nozzle area and then creep upstream into the chamber in the form of a central recirculation zone (CRZ), namely, with successive increases in β. Pursuant to a cursory investigation, the backflow formation is found to be quite sensitive to geometric changes, especially in the nozzle region. To suppress the emergence of recirculation zones, we find it useful to increase the radius of curvature at the nozzle entrance region to the extent of reducing the final outlet diameter. In this process, we identify a nozzle inlet curvature, $RN=0.12a/(1−β)$⁠, that will prevent the backflow from occurring. Unsurprisingly, increasing this curvature eliminates the backflow that has plagued numerous cyclonic flow investigations; however, it also causes the mantles to shift radially inwardly while compelling the corresponding axial velocities to increase near the centerline, as in the classical case of vena contracta. Naturally, any area contraction must be accompanied by an increase in the axial velocity magnitude, which is necessary to evacuate the same amount of fluid across a narrower passage. This trend is further reflected in the higher axial velocities observed in both CFD and analytical models of the Beltramian motion with $β=0.628$⁠, i.e., when gauged against the computational and complex-lamellar solutions with a larger $β=0.707$⁠. Finally, as far as the development of multiple mantles is concerned, we initially speculate that multiple flow passes may be triggered by the narrowing of the outlet diameter. Nonetheless, despite the pre-selection of outlet diameters that are guided by the theoretical solutions of Vyas et al.,⁴⁸ we are unable to induce multiple mantle combinations. Upon careful analysis and review of the available data, we now speculate that the inconclusive indicators of perceived multiple mantles that have been reported in former experimental and numerical studies, such as those by Vyas et al.,⁴⁸ Anderson et al.,²⁹ Rom, Anderson, and Chiaverini,³⁰ Rom,⁵⁶ and others⁴³ are likely caused by the presence of backflow and CRZ development. In actual numerical simulations, we find that the imposition of either velocity or pressure inlet conditions with a properly rounded outlet nozzle will effectively mitigate the evolution of multiple mantles and backflow regions. Now that the backflow problem has been essentially resolved, it is possible to conduct, in future work, a wide range of parametric trade permutations that include the dynamic flow control parameters to the extent of making it feasible to determine the threshold values of the non-dimensional parameters that affect the inception, precession, and breakdown of wall-bounded cyclonic flows under both laminar and turbulent conditions.

This work was supported partly by the National Science Foundation, through Grant No. CMMI-1761675, and partly by the Hugh and Loeda Francis Chair of Excellence, Department of Aerospace Engineering, Auburn University. The authors are deeply indebted to M. J. Chiaverini, Director of Propulsion Systems at Sierra Space Corporation, and to D. Benner, B. Pomeroy, and A. Sauer, for numerous technical exchanges and for their unwavering support of their cyclonic flow investigations.

The authors report no conflict of interest.

The data that support the findings of this study are available from the corresponding author upon reasonable request.