Paradoxical predictions of swirling jets

Abstract This paper examines the shape of a steady jet with a swirling component, ejected from a circular orifice at an angle to the horizontal. Assuming the Froude number to be large, we derive a set of asymptotic equations for a slender jet. In the inviscid limit, the solutions of the set predict that, if the swirling velocity of the flow exceeds a certain threshold, the jet bends against gravity and rises until the initial supply of the liquid's kinetic energy is used up. This effect is due to the fact that the contributions of the swirl and streamwise velocities to the centrifugal force are of opposite signs, with their sum to be balanced by gravity. As a result, swirl- and streamwise-dominated jets bend in opposite directions. Downward-bending jets also exhibit counter-intuitive behaviour. If the swirling velocity is strong enough (but is still below the above threshold), the streamwise velocity on the jet's axis may decrease with the distance from the orifice, despite the acceleration due to gravity. Eventually, a stagnation point emerges due to this effect, arguably destabilising the jet. Also paradoxically, viscosity-dominated jets can reach higher (if they bend upwards) and farther (in all cases) than their inviscid counterparts, due to the fact that viscosity suppresses formation of stagnation points.


Introduction
Straight and nearly straight jets have been studied for almost 150 years, since the seminal works of Plateau (1873) and Rayleigh (1878). Curved jets, on the other hand, have a much shorter history -probably because they cannot be examined using the usual (cylindrical) coordinates. Alternative curvilinear coordinates associated with the jet centreline were suggested by Entov & Yarin (1984); even though they are non-orthogonal and, thus, make the Navier-Stokes equations several pages long, they have been successfully applied to gravity-affected slender jets (Wallwork 2001;Wallwork et al. 2002;Shikhmurzaev & Sisoev 2017;Decent et al. 2018).
The present paper considers jets where the velocity involves a non-zero swirling component, so that the trajectories of fluid particles are spirals; such jets can be created by † Email address for correspondence: eugene.benilov@ul.ie inserting a propeller into the orifice. We are concerned with oblique swirling jets: those ejected in a direction other than the vertical, so that their trajectories are curved by gravity. A set of asymptotic equations have been derived for slender jets, i.e. those whose radii of curvature exceed their thickness.
Three highly counter-intuitive results will be presented in this paper.
(i) We show that the Navier-Stokes equations admit asymptotic solutions describing swirling jets which bend upwards, i.e. against gravity.
They have yet to be tested for stability, but even if they turn out to be unstable, downward-bending jets are not the alternative, as they do not exist in this range of parameters. The only realistic alternative is unstable, eternally evolving flow.
(Note that solutions describing upward-bending jets have been found before -and not only jets (Wallwork 2001), but also liquid curtains (Keller & Weitz 1957;Benilov 2019Benilov , 2021. These jets and curtains, however, bend upwards due to surface tension, whereas the jets found in this work do so even if the capillary coefficient is zero. It is the swirling velocity that makes them bend upwards, not surface tension.) (ii) We show that the streamwise velocity in some downward-bending jets decreases with the distance from the orifice -to the extent that a stagnation point (SP) emerges in the flow. This behaviour is at odds with the intuitive perception that falling objects accelerate due to gravity. It is further argued that SPs cause instability, so that downward jets developing them terminate due to an instability at a finite distance from the orifice.
(Note that the emergence of SPs in swirling flows (not necessarily jets) has been reported before, and is usually referred to as 'vortex breakdown' (e.g. Benjamin (1962), Escudier (1988), Brown & Lopez (1990), Gelfgat, Bar-Yoseph & Solan (1996), Billant, Chomaz & Huerre (1998), Ruith et al. (2003), Moise & Mathew (2021), Shtern (2012) and references therein). So far, it has been observed in flows not affected (thus, not accelerated) by gravity, so their slow-down is less counter-intuitive -but it is still the same effect as the one described here. The mere fact that our asymptotic model describes a paradoxical, but known, phenomenon, indirectly validates the rest of our results.) (iii) Viscosity-dominated jets may reach higher (if they bend upwards) and propagate farther (in all cases) than their inviscid counterparts. The reason is that the latter are prone to be terminated by SPs, whereas the former are affected by cross-stream homogenisation (either eliminating, or at least postponing, formation of SPs). This paper has the following structure. In § 2, we formulate the problem for an inviscid fluid without surface tension, and in § 3, derive a set of asymptotic equations describing jets with a large Froude number. Examples of solutions of these equations are obtained in § 4, and the results obtained are extended to viscous capillary fluids in § 5. In § 6, we compared our results with those obtained previously for liquid curtains.
2. Formulation of the problem: ideal fluids without surface tension 2.1. Preliminaries Consider a steady flow of fluid of density ρ ejected with a flux F from a circular orifice of radius R 0 . Choosing the characteristic velocity scale to be Figure 1. The setting: a jet ejected from a circular outlet. Here (l, r, φ) are the curvilinear coordinates associated with the jet's centreline, and α(l) is the angle between the centreline and the horizontal.
we introduce the reciprocal of the Froude number, where g is the acceleration due to gravity. It can be shown that the characteristic radius of the jet's curvature due to gravity is Everywhere in this paper, we assume that ε 1 -hence, L R 0 and the slender-jet approximation can be employed.
We shall use the curvilinear coordinates (l, r, φ) related to the Cartesian coordinates (x, y, z) by x =x + r sin α sin φ, y = −r cos φ, z =z − r cos α sin φ, (2.4a-c) wherex(l) andȳ(l) are the coordinates of the jet's centreline, and α(l) is the local angle between the centreline and the horizontal (see figure 1). Observe that the second expression in (2.4a-c) does not include the lateral displacementȳ(l), so that all of the centreline is in the vertical (x, z) plane. It can be readily verified that (l, r, φ) are orthogonal coordinates (which they would not be, shouldȳ(l) be introduced in (2.4a-c) (Entov & Yarin 1984;Shikhmurzaev & Sisoev 2017;Decent et al. 2018)). To identify l with the centreline's arclength, we require Equations (2.5a,b) and boundary conditions (2.6a,b) relate uniquely (x,z) to α, whereas the latter remains mathematically undetermined. It will be fixed later as convenient.
2.2. Governing equations In order to minimise straightforward, yet extremely tedious, algebra associated with the use of curvilinear coordinates, we shall first consider the simplest setting: an ideal fluid without surface tension. Once this limit has been examined, we shall briefly explain how the results can be extended to the general case of viscous capillary jets.
The scaling in the present problem is similar to that used by Benilov (2019) for liquid curtains. In application to jets, it amounts to introducing the following non-dimensional variables: where (u l , u r , u φ ) are the projections of the velocity onto the local axes of the curvilinear coordinates and p is the pressure.
In terms of the non-dimensional variables, the Euler equations have the form (the subscript nd omitted) where the Lamé coefficient for the variable l is The other two Lamé coefficients have been replaced with their values, h r = 1 and h φ = εr. Let the jet's free boundary be determined by the equation r = R(l, φ), where R satisfies the kinematic boundary condition, (2.14) The dynamic boundary condition, in turn, is (2.15)

Asymptotic analysis
3.1. Asymptotic equations Within the framework of the slender-jet approximation, the jet's cross-section is almost circular, and the radial velocity u r is much smaller than u l and u φ . Thus, we seek a solution of the form Substituting these expansions into the exact problem (2.9)-(2.15) and keeping the leading order only, we obtain Evidently, we do not have enough equations to determine the leading-order solutionhence, we have to consider the next order. The next-to-leading order yields the following equations: and the following boundary conditions: Boundary-value problem (3.4)-(3.8a,b) is consistent with the following ansatz: where the functions with the superscripts (1,0) and (1,1) satisfy two non-coupled problems, Expressing u (1,0) r from (3.10) and substituting it into (3.12)-(3.14), we obtain Observe that (3.22)-(3.24) include only the leading-order unknowns and, thus, can be viewed as the solvability conditions for boundary-value problem (3.10)-(3.15).

Discussion: singular points of the equations derived
Observe that (3.23) and boundary condition (3.24) appear to be singular if/where ∂u (0) l /∂r = 0. It turns out, however, that the apparent singularity of the boundary-value problem does not translate into singularity of the solution.
To understand why, integrate (3.23) from r = 0 to r = R (0) , take into account (3.24) and thus obtain d dl Integrating this equation with respect to l and equating (without loss of generality) the constant of integration to unity, we obtain This (non-singular) equality reflects conservation of the mass flux. In what follows, it will replace the singular condition (3.24).
To resolve the singularity in (3.23), we note that the asymptotic boundary-value problem is invariant with respect to simultaneously changing Thus, u (0) l and p (0) can be analytically extended to the region r < 0 as even functions, whereas u (0) φ can be extended as an odd function. Since all three should be analytic at r = 0, this implies Then, as follows from (3.22), Keeping this in mind, we integrate (3.23) from r = 0 to r = r , interchange r ↔ r , and obtain where functions without explicitly stated arguments depend on (l, r). Using (3.34), we reduce (3.22) to In what follows, (3.34)-(3.35) replace (3.22)-(3.23).

Inviscid jets
To simplify the notation, we change Now, (3.2), (3.28)-(3.35) and boundary condition (3.28) become The solution should be analytic at r = 0, which implies It can be shown that, on the (l, r) plane, boundary-value problem (4.2)-(4.7) is hyperbolic and, as such, requires entry conditions at the orifice. We assume that where the jet's initial radius matches that of the orifice, and [u 0 (r), w 0 (r)] and α 0 are the ejection velocity and angle, respectively. Note that u 0 (r) should comply with (4.6), i.e. This requirement is not a restriction, as it can be always enforced via a proper choice of the scale U during non-dimensionalisation. Note also that, near the orifice, the jet may experience certain adjustments -e.g. slight reduction of its radius. These adjustments are not described by the slender-jet approximation, but they can be still accounted for by changing entry conditions (4.9a-c)-(4.10). In this case, the point l = 0 would correspond to a cross-section located at a certain distance from the actual orifice. Finally Exact results Some of the mathematical properties of boundary-value problem (4.2)-(4.13a,b) have far-reaching physical consequences. (4.14) (4.5) implies that dα/dl > 0. Thus, if the swirling velocity w is sufficiently large, the jet bends upwards, against gravity.
To understand the physics behind such a counter-intuitive behaviour, we rewrite (4.5) as a balance of the centrifugal force and gravity, The right-hand side of this equality represents the cross-stream component of the gravitational force acting on a cross-section of radius R -and the first/second term on the left-hand side describes the contribution of the swirl/streamwise velocity to the centrifugal force. To understand the structure of the latter is easy: one should only keep in mind that dα/dl is the curvature of the jet's centreline, and the minus in front shows that the centrifugal force is directed away from the centre of curvature. As for the swirl contribution to the centrifugal force, it is also proportional to dα/dl (because the forces affecting individual particles in a straight swirling jet average out to zero). Most importantly, this term's sign is opposite to that of its streamwise counterpart.
To understand why, consider, say, an upward-bending jet. It is clear from basic geometry that its upper half slightly shrinks, whereas the lower half expands -hence, the trajectories of particles swirling through the former are curved more than those in the latter. Due to this difference and the fact that the centrifugal force is proportional to the curvature of the particles' trajectories, the swirl part of the force affecting the whole cross-section is directed upwards, i.e. opposite to the streamwise part.
Thus, the upward bending of a swirling jet could be understood through an analogy with a system of two interconnected bodies: one with positive mass and another, with negative. If the absolute value of the latter exceeds that of the former, gravity makes the whole system soar.
(ii) Given the counter-intuitive nature of property (i), it should be reassuring to ascertain that our equations are consistent with basic physical principles. Two of these are discussed below.
(a) Boundary-value problem (4.2)-(4.13a,b) is consistent with the following ansatz: where u(l), α(l) and R(l) satisfy (4.17a-c) Using (4.12a,b)-(4.13a,b), one can show that the above equations describe free-falling liquid, such that the particles follow parabolic trajectories. (b) It follows from (4.2)-(4.13a,b) that d dl R 0 u 2 + w 2 2 + p + z ur dr = 0, (4.18) which reflects conservation of the energy flux through the jet's cross-section. This property of the asymptotic model is helpful for controlling the accuracy of numerical simulations. One can also show that d dl These 'local' conservation laws do not seem to have physical meaning, but one of them happens to be useful anyway. (iii) Important information can be obtained by relating ∂u/∂l, ∂w/∂l and dR/dl to u, w and R (in (4.2)-(4.13a,b), the former are related to the latter implicitly). Observe that, if a pair (l, r) exist such that u(l, r) = 0 the coefficients in (4.20) become singular. In what follows, such pairs will be called SPs even if w(l, r) / = 0. (iv) It follows from (4.20) that ∂u/∂l is singular at SPs (see Appendix A.1). The following physical aspects of these singularities are worth noting.
(a) Given the hyperbolic nature of our equations and the fact that the streamlines are their characteristics, it comes as no surprise that the solution becomes singular where one of the characteristics 'stops'. (b) If viscosity is introduced, the solution becomes regular (due to the fact that the governing equations become parabolic, see § 6 below). (c) The presence of an SP invalidates -at least, locally -the assumption that the Froude number be large (which underlies all of our results). Whether this invalidates the global solution is unclear -but is also unimportant in view of item (d), below. (d) Stagnation points in ideal fluids often make the flow unstable (Friedlander & Vishik 1991), and even in viscous fluids they are a destabilising influence (Benilov & Lapin 2014). It is safe to assume that an inviscid or weakly viscous jet exists only between the orifice and the first SP, and it breaks up immediately after it.

Numerical results
Boundary-value problem (4.2)-(4.13a,b) was solved numerically using the method of lines (Schiesser 1978): the r axis was discretised, so that the governing (partial differential) equations turned into a large set of ordinary-differential equations. These equations are implicit with respect to the l-derivatives, which were solved for (at each step in l) using the MATLAB function FSOLVE. The actual integration of the equations was carried out using the function ODE15s (an adaptive-step solver designed for stiff problems). For reasons explained above, the solutions were terminated at the first SP (if any).  (1) is the only one that continues indefinitely, whereas curve (5) continues beyond the boundary of this figure, but eventually stops.
Several particular cases of entry conditions where simulated, but the results will be presented only for the simplest one, where the value of u 0 was chosen to satisfy constraint (4.11). If there were no gravity, conditions (4.25a,b) would describe a non-sheared solid-body-rotating flow (as seen later, this is convenient when comparing inviscid and viscous jets). According to condition (4.14), upward-bending jets exist if and only if ω 0 > 4. (4.26) Figure 2 shows the trajectories of four jets with the same ejection angle, but different ω 0 . The following features can be observed.
(i) As expected, the two jets satisfying condition (4.26) bend upwards.
(ii) Out of the four jets depicted, three are terminated by SPs (a) In the downward-bending jet with ω 0 = 3 and upward-bending jet with ω 0 = 7, the SPs emerge on the jet's axis (see figures 3b and 4b). (b) In the upward-bending jet with ω 0 = 5, the SP emerges somewhere between the axis and free surface (unfortunately, we have been unable to finish this computation due to numerical instability, so figure 3(a) shows only the tendency of developing an SP). (iii) Not all downward-bending jets develop SPs, as illustrated in figure 4(a) for the jet with ω 0 = 1. The fluid particles in this case do accelerate on the way down. Figure 5 illustrates how the behaviour of jets depends on the ejection angle -with the conclusion that α 0 affects the trajectory only near the orifice, but the rest of the behaviour is determined by ω 0 . However, there is a counter-intuitive feature in this case too, as illustrated in figure 6. The jet depicted in this figure is one of those that should bend downwards and develop an SP -but the ejection angle is positive -so, initially, the jet is rising. Given this, its streamwise velocity should be decreasing everywhere in the jet's cross-section, but figure 6 shows that the velocity on the jet's axis initially grows  and reaches its maximum at the top of the trajectory. After that, the jet switches to the behaviour expected for this value of ω 0 and maintains it until it is terminated by an SP. We simulated numerous other entry conditions, and the following tendencies were observed.
(i) If condition (4.14) holds at the orifice, it continues to hold for the rest of the trajectory. Thus, if a jet bends upwards near the orifice, it bends upwards everywhere. (ii) Upward-bending jets typically develop SPs (either on or off the axis). In some cases, however, it is impossible to actually see the SP, as our numerical method becomes unstable due to the emergence of a region with highly sheared w. (Refining the mesh did not help in this case, suggesting that the root cause of the instability is a physical one, e.g. a tangential jump emerging in the cross-stream profile of w. Unfortunately, we have been unable to verify this hypothesis beyond reasonable doubt, as this would require a completely different numerical method: one suitable for weak solutions. Furthermore, even if tangential jumps do occur in reality, they are closely followed by emergence of an SP which breaks the jet up anyway.) (iii) Downward-bending jets with α 0 < 0 develop SPs only at r = 0, and only if ∂u ∂l < 0 at l = 0, r = 0, (4.27) i.e. at the centre of the orifice, particles decelerate in the streamwise direction.
Conclusion (iii) allows one to use boundary-value problem (4.20), (4.22)-(4.23) to derive a sufficient criterion for detecting SPs in downward bending jets. For entry conditions (4.25a,b) and α 0 < 0, criterion (4.27) amounts to (see Appendix A.2) This condition agrees with our numerical results very well; and not only for negative α 0 , but also for most positive ones (in which case criterion (4.28) is formally inapplicable). This is probably due to the fact that, for a downward-bending jet, an initially positive α rapidly becomes negative, so that u and w remain close to their entry profiles.

4.3.
Approximate results: near-critical jets Our asymptotic model (4.2)-(4.13a,b) was derived under the sole assumption that the Froude number is large, but the derived equations are not simple enough to be solved analytically. Thus, it is worthwhile to try to find some extra assumptions allowing one to further simplify (4.2)-(4.13a,b). To avoid confusion, the results of such analyses will be called 'approximate', with the word 'asymptotic' reserved for boundary-value problem (4.2)-(4.13a,b) itself. Denote Admittedly, the extra assumption can make our asymptotic model inconsistent due to the fact that some of the retained terms can now be smaller than the omitted ones. Since the latter are O(ε), one can ensure consistency by requiring |δ| ε. This constraint could be replaced with a weaker one, |δ| ε, if the near-critical regime were examined using the exact Navier-Stokes equations (which is, however, beyond the scope of the present work).
Numerical experiments show that near-critical jets evolve in two stages. First, only α varies with l, whereas u, w, p and R remain virtually unchanged -i.e. the jet's trajectory is curved, but all other characteristics are close to their entry values. Second, α approaches either π/2 or −π/2, after which the rest of the characteristics start to evolve.
The second stage is difficult to examine analytically, but the first stage, is easy (see Appendix B). The approximate results have been compared with the numerical solution of asymptotic problem (4.2)-(4.13a,b) for entry conditions (4.25a,b), see in figure 7. One can see that the approximate solution predicts the shape of the trajectory well even if δ = 1.
Note that the approximate results assume that u is close to its entry value: hence, they are not applicable near the SP. However, if δ is indeed small, the approximate solution does predict the position of the stagnation (termination) point reasonably well (see figure 7).

The asymptotic equations
It is well known that the behaviour and characteristics of an inviscid flow can be qualitatively different from those of a viscous one. Thus, it is essential to verify whether upward-bending jets and jets with SPs can be observed in the presence of viscosity.
Let the fluid's kinematic viscosity and surface tension be ν and σ , respectively. Then, in addition to the Froude number, two extra non-dimensional parameters arise: the so-called 'reduced Reynolds number' and the Weber number. We shall use their reciprocals, namely, where L is defined by (2.3). Omitting the details (which are similar to those of the inviscid non-capillary case (see Dubovskaya 2020)), we only formulate the resulting asymptotic equations.
To incorporate viscosity and surface tension in set (4.2)-(4.13a,b), one needs to replace (4.2)-(4.3), (4.5) and boundary condition (4.7) with respectively. Accordingly, the condition of existence of upward-bending jets becomes This condition does not involve the viscosity coefficient μ, whereas an increase in the capillary coefficient γ makes it easier to hold. However, jets with γ = O(1) break up near the orifice due to the Plateau-Rayleigh instability. Thus, in what follows, we let γ = 0, i.e. assume that surface tension is negligible.
To illustrate this assumption, consider water at 20 • C, in which case ρ = 998.2 kg m −3 , ν = 1.002 × 10 −6 m 2 s −1 , σ = 72.8 mN m −1 . (5.7a-c) Let the jet parameters correspond to those of a garden hose, say, (5.8a,b) in which case μ ≈ 4.1 × 10 −3 , γ ≈ 1.5 × 10 −2 . (5.9a,b) Since μ and γ are small, viscosity and surface tension can be both neglected (i.e. these effects influence the jet at a much larger distance from the orifice than gravity). If, however, water is replaced with glycerol, then and we obtain (for the same jet parameters) μ ≈ 4.6, γ ≈ 1.0 × 10 −2 . (5.11a,b) In this case, viscosity has to be accounted for, whereas surface tension can still be neglected. Generally, for different liquids and jets, μ and γ vary in a wide range and, thus, should be estimated on a case-to-case basis.

Viscosity-dominated jets
We shall examine viscous jets under an extra assumption, μ 1. (5.12) Observe that the leading-order approximation of boundary-value problem comprising (4.4), (4.6), (4.8a,b) and (5.2)-(5.5a-c) is consistent with the following ansatz: which describes an almost shearless/solid-body-rotating flow, with the streamwise velocity u(l) and angular velocity ω(l) varying along the jet. Ansatz (5.13a,b)-(5.14a,b) automatically eliminates SPs, but allows u to vanish in a certain cross-section of a jet, with R simultaneously tending to infinity (the slender-jet approximation fails some distance before that, of course). In order to derive asymptotic equations for u(l) and ω(l), one should generally examine higher orders of the perturbation theory. In the present case, however, there is a simpler option: the governing equations can be rearranged in such a way that the leading-order terms are eliminated. Then, substitution of leading-order solution (5.13a,b)-(5.14a,b) into the rearranged equations will immediately yield the desired equations for u(l) and ω(l).
Following this plan, consider Integrating by parts, taking into account boundary conditions (5.5a-c) and setting γ = 0, we observe that the O(μ) terms disappear, and we obtain Substituting into these equations ansatz (5.13a,b)-(5.14a,b) and keeping the leading-order terms only, we obtain, after straightforward algebra, which allows one to deduce that all viscosity-dominated upward-bending jets reach the same maximum (non-dimensional) height, z max = 2. Equations (5.22)-(5.23) can be readily integrated numerically, and examples of solutions are presented in figure 8 together with the corresponding solutions of the inviscid boundary-value problem (4.2)-(4.13a,b). Interestingly, viscosity-dominated jets may propagate farther and reach higher than their inviscid counterparts: simply because the latter are terminated by SP, whereas the former, cannot.
Note also that the leading-order ansatz (5.13a,b)-(5.14a,b) satisfies the energy conservation law (4.18): hence, energy dissipation due to viscosity is a higher-order effect.

Concluding remarks: jets versus liquid curtains
The main question associated with our results is whether or not upward-bending jets can be observed in an experiment. Since a similar issue has been debated in the context of liquid curtains it is worthwhile to briefly review this work and see if the conclusions drawn for capillary curtains apply to swirling jets.
Liquid curtains are flows originating from a long slot (outlet), the same way jets originate from a circular orifice. Solutions describing upward-bending curtains (which exist when the Weber number We is less than unity) were originally found by Keller & Weitz (1957), and then rediscovered in a more general formulation by Benilov (2019). The latter work, however, was criticised by Weinstein et al. (2019), who pointed out that the hyperbolic equations used by Benilov are such that, for upward-bending curtains, one of the characteristics corresponds to (capillary) waves propagating upstream. As a result, one of the boundary conditions set at the outlet effectively constrains events occurring in the past and, thus, contradicts the causality principle.
This issue was re-examined by Benilov (2021) for the simplest setting comprising an ideal fluid and almost shearless curtains. It was further assumed that We ≈ 1, in which case the phase velocity of the upstream-propagating waves is small, so that their dispersion becomes important. It was claimed that, since the velocity of dispersive waves is not bounded above, events occurring anywhere in the curtain are immediately sensed near the outlet, thus resolving the conflict with the causality principle.
Mathematically, the asymptotic equation derived by Benilov (2021) for curtains with We ≈ 1 is not hyperbolic: as a result, the concept of characteristics becomes irrelevant. Overall, it was argued that upward-bending curtains with We ≈ 1 could be observed experimentally, but curtains with We differing from unity by O(1) are likely to be unstable.
It turns out that none of the results obtained by Benilov (2021) can be immediately extended to swirling jets. To do so, we need an evolutionary model for jets, not just the one describing steady states. If the former were available, we would be able to see if it is hyperbolic in the (t, l, r) space and, if it is, where its characteristics 'go'. We would also be able to examine the parameter region where the phase velocity of waves vanishes (which would be similar to the region We ≈ 1 for curtains).
Unfortunately, there is no simple way to derive an evolutionary model for a swirling jet. The problem is that, unlike curtains, a jet cannot move just in the vertical (x, z) plane: as soon as it starts doing so, the gyroscopic force pushes it laterally. To account for the lateral motion, the curvilinear coordinates (2.4a-c) have to be modified by introducing the displacementȳ(l, t) -which, however, makes the modified coordinates non-orthogonal (Entov & Yarin 1984;Shikhmurzaev & Sisoev 2017). As a result, one has to either cope with the extremely cumbersome algebra due to the use of non-orthogonal coordinates or try to extend the two-dimensional approach of Benilov (2019) to three dimensions. (Benilov (2019) treated the Cartesian coordinates as extra unknowns (depending on the curvilinear coordinates and time), and the orthogonality conditions, as extra governing equations.) The latter appears to be simpler.
Once the technical framework for evolving swirling jets is in place, one can use it to study their stability (so far, this issue has been examined only for straight jets with strong surface tension (e.g. Ponstein 1959;Kubitschek & Weidman 2007;Eggers & Villermaux 2008), which are irrelevant to the problem at hand).
Another potentially important modification of the present setting consists in taking into account of the drag force created by the surrounding air. Its density and viscosity are small, but the stress it generates can be important, especially, for the swirl component of the flow (which is highly sheared).
Funding. We acknowledge the financial support of Science Foundation Ireland under grant number SFI/13/IA/1923.

Declaration of interests. The authors report no conflict of interests.
Author ORCIDs.

Appendix A. Stagnation points and singularities
A.1. The flow near as SP Two kinds of SPs can be distinguished: those located on the jet's axis and those located elsewhere. (In three dimensions, the latter are actually rings, but we will call them points anyway, as that is what they are on the (l, r) plane.) The latter will be examined first.
where u 2 , w 0 and w 1 are constant coefficients. Substituting these expansions into (4.20), one can readily find the asymptotics of the general solution for v, v ∼ C 1 r − r sp cos √ A r − r sp + C 2 r − r sp sin √ A r − r sp as r → r sp , where C 1,2 are constants of integration and which shows that the streamwise derivative of the streamwise velocity is indeed singular.
(ii) Let an SP be located at (l sp , 0). This time, we need the behaviour of the solution near the SP, not just at it; hence, we let u = u 0 + u 2 r 2 + O(r 4 ), w = w 1 r + O(r 3 ) as r → 0, ( where C is a constant of integration (note that C / = 0, otherwise v would be zero for all r and, thus, would not be able to satisfy the boundary condition (4.23) unless α = 0), then (4.24) yields ∂u ∂l ∼ 2C + O(r 2 ) as l → l sp , r → 0.
Thus, the asymptotics of ∂u/∂l at the SP is finite (unlike that in the case r sp / = 0 examined previously). This does not mean, however, that it remains finite after the SP, which can be proved by contradiction.
Assume that ∂u/∂l remains finite in a certain neighbourhood of r = 0, for a certain interval after l = l sp . Since ∂u/∂l was negative before the SP, it will remain negative after it -which means that, for l > l sp , the SP moves away from the axissay, to r = r sp . This situation, however, is different from case (i) , because now the asymptotics of u and w are (compare these equalities with (A1)). Substituting (A11) into (4.20), one can readily find the asymptotics of the general solution for v, v ∼ C 1 r − r sp 1/2+1/2 √ 1−4A + C 2 r − r sp where A = 2w 0 r sp u 2 1 w 1 + w 0 r sp .
Finally, it follows from (4.24) that, regardless of the value and sign of A, ∂u/∂l is singular as r → r sp -which is what we set out to prove.