On the axisymmetric steady incompressible Beltrami flows

In this paper, Beltrami vector fields in several orthogonal coordinate systems are obtained analytically and numerically. Specifically, axisymmetric incompressible inviscid steady state Beltrami (Trkalian) fluid flows are obtained with the motivation to model flows that have been hypothesized to occur in tornadic flows. The studied coordinate systems include those that appear amenable to modeling such flows: the cylindrical, spherical, paraboloidal, and prolate and oblate spheroidal systems. The usual Euler equations are reformulated using the Bragg--Hawthorne equation for the stream function of the flow, which is solved analytically or numerically in each coordinate system under the assumption of separability of variables. Many of the obtained flows are visualized via contour plots of their stream functions in the $rz$-plane. Finally, the results are combined to provide a qualitative quasi-static model for a progression of flows that one might see in the process of a vortex breakdown. The results in this paper are equally applicable in electromagnetics, where the equivalent concept is that of a force-free magnetic field.


Introduction
In fluid mechanics, Beltrami or helical flows are fluid flows in which the velocity and the vorticity (curl of velocity) of the fluid are parallel to each other at all points and all times. Flows of this nature have been studied since at least the late 1800s and have applications in both fluid dynamics and electromagnetics, where How to cite this paper: Bělík, P., Su, X.Q., Dokken, D.P., Scholz, K. and Shvartsman, M.M. (2020) On the Axisymmetric Steady a force-free magnetic field is one for which the Lorentz (magnetic) force density vanishes, or equivalently, the magnetic field is everywhere parallel to the direction of the current flow [1]. In this paper we will focus mainly on the hydrodynamics case, but many parallels can be drawn between the two, and results from one field can be applied in the other one. An alternative to this approach is to study generalized Beltrami flows in which the curl of the cross product of the velocity and vorticity is zero, rather than the cross product itself [2] [3].
Beltrami fluid flows are of interest for several reasons. They can have complex dynamics [4], and types of Beltrami flows that possess ergodic theoretic properties, such as strong mixing, make for attractive models for turbulent flows [5]. Additionally, "every incompressible fluid flow is a superposition of Beltrami flows" [4]. In nature, classes of rotating thunderstorms (supercells) may exhibit characteristics of Beltrami flows [6]; this is supported by numerical simulations such as in [7]. Flows with low instability (low convective available potential energy or CAPE) and nearly circular hodographs approach Beltrami flows as the hodograph becomes more circular [8]. Highly helical flows are thought to be present in well-developed tornadic flows [9] or resemble the gross supercell structure [10] [11] [12]. Beltrami flows that are solutions of the stationary Euler equations and the decaying nonsteady Beltrami flow solutions of the Navier-Stokes equations (e.g. [11] [12] [13]) could be or have been used to test code in numerical weather models, such as in the Advanced Regional Prediction System (ARPS) [14], the Weather Research and Forecasting Model (WRF) [15], or the Terminal Area Simulation System (TASS) [16]. Additionally, Beltrami flows have been shown to be the equilibrium states of statistical mechanics systems consisting of simplified axisymmetric Euler solutions characterized by only three conserved quantities: microscopic energy, helicity, and angular momentum [17] [18] [19].
Work on Beltrami flows was initiated in the works of Eugenio Beltrami and Ippolit Gromeka in the late 1800s. Interesting explicit incompressible flows in Cartesian coordinates have been found since then including the famous Gromeka-Arnold-Beltrami-Childress (GABC) flow, the cat's eye 3D flow, and others.
These flows and a theorem for construction of other flows can be found in [5].
In this paper we analytically and numerically investigate axisymmetric incompressible inviscid steady state Beltrami (specifically, Trkalian, see below) flows with the goal of potentially developing other test cases for numerical models as well as possible models for tornadic or supercell flows. Steady axisymmetric solutions to Euler's equation have been studied and their general form established in terms of the stream function, angular momentum, and potential vorticity in [17]. Numerical approaches to compute axisymmetric flows with swirl have also been explored (see, e.g. [22]). We explore several geometries characterized by various orthogonal coordinate systems and construct separable Beltrami solutions to the relevant equations in these systems. Some of these solutions are P. Bělík et al. Open Journal of Fluid Dynamics known (cylindrical and spherical coordinate systems), while we are not aware of the existence in the literature of the other solutions developed in this paper (paraboloidal and the two types of spheroidal coordinate systems). We further explore symmetries in the mathematical problem and outline a way to construct infinitely many other Beltrami flows in each coordinate system. We present graphical results of such flows in each coordinate system. Our approach here shares some similarity with that in [23] and its extension in [24], but our emphasis is on developing new purely Beltrami flows.
The paper is organized as follows. In Section 2 we review the governing equations for inviscid fluid flows and discuss Beltrami flows. In Section 3 we reformulate the axisymmetric problem using the Bragg-Hawthorne equation and discuss additional symmetries that can be used to construct new solutions from existing ones. In Section 4, several orthogonal coordinate systems of interest are introduced, and under the assumption of separability of variables the mathematical problem is reformulated as a coupled system of two second-order ordinary differential equations, which are equipped with either zero, one, or two boundary conditions. In several cases, one of the equations results in a boundary value eigenvalue problem. These equations are then analyzed and in some cases analytic solutions are presented. Many of these solutions are given in terms of special functions (classical hypergeometric and Kummer functions) and we refer the reader to [25] for details. In the remaining cases where analytic solutions do not appear available (prolate and oblate spheroidal coordinates in Sections 4.4 and 4.5), a numerical approach is used to approximate the solutions. In Section 5 we combine several solutions from three of the coordinate systems and discuss their similarity to swirling tornadic flows [26] [27].

Beltrami Flows
Inviscid fluid flows are usually modeled by the Euler Equation (plus relevant energy and state equations) is the velocity of the fluid at some point is the mass density, g is the acceleration due to gravity, and z is the upward pointing unit coordinate vector. Gradients are taken with respect to the spatial variable Taking "curl" of the equation for velocity (1), we obtain an equation for the vor-  [5] for some 0, where α is called the abnormality and is, in general, a function of position and time. This condition implies that It is shown in [5] that a steady Beltrami pair u , ω for which solution to the Euler Equation (1) with constant density ρ ; the mass conservation Equation (2) is then trivially satisfied as well. The vorticity Equation (3) then implies that such a flow will be what is known in the atmospheric science community as "barotropic", since the "baroclinic" term p ρ ∇ ×∇ will vanish.
Also, taking divergence of both sides of (4) and using 0 ∇ ⋅ = u results in a necessary condition on the function α 0, α ⋅∇ = u so α is constant on the streamlines of the flow.
A steady Beltrami pair u , ω for which The special case of (4), in which α is independent of the spatial variables, is known as a Trkalian flow [28] [29]. In this case, α is also independent of time [29]. Our work focuses exclusively on Trkalian flows.

The Bragg-Hawthorne Equation
We now formulate the basic equations governing incompressible, steady, axisymmetric Eulerian flows with constant mass density and in the absence of body forces. Additional details can be found, for example, in [30] [31]. . The corresponding velocity of the flow is then Notice the similarity with the laplacian of ψ which in the axisymmetric case is We let 2 1 and 2 p C rw H ρ = = + u (9) be the swirl (sometimes referred to as angular momentum) and energy head, respectively, where ( ) , p r z is the fluid pressure and ρ its mass density. The conditions for a steady flow are [30] ( ) ( )  In what follows, we will focus on the simple case studied, for example, in [31]. In this case, α = u ω , so the proportionality constant is independent of r and z and the flow is Trkalian. The Bragg-Hawthorne Equation (10) then further reduces to where the operator 2 D is defined in (7). Note the similarity to the well-studied

Symmetry Transformations
Note that Equation (12) possesses several properties that can be used to construct additional Beltrami solutions: 1) By linearity and homogeneity of (12), a linear combination of two solutions is another solution.

2) If
( ) ζ ∈  . This follows since only a second derivative in z is present in the operator is a solution and The last three properties can be used to construct flows in the upper half space that satisfy the boundary condition that the flow cannot penetrate the ground.

Equations and Solutions in Various Coordinate Systems
In this section we explore the various appropriate coordinate systems, rewrite the governing Equation (12) in these systems, and, under the assumption of separability of the sought solution, rewrite the problem as a system of two independent ordinary differential equations to be solved.  Figure 1(c)), and eventually bifurcates into several vortices as swirl increases (Figure 1(d)), we will focus on coordinate systems that seem most "friendly" to modeling such flows. Specifically, after discussing cylindrical and spherical coordinate systems, we will focus on paraboloidal and prolate and oblate spheroidal coordinate systems. More about various characteristics of rotating flows in nature as swirl changes can be found in [6] [26], and an equivalent process observed in a tornado vortex simulation chamber is described in [27].

Cylindrical Coordinates
This is a well-studied case and various solutions appear in the literature (see, e.g. [1] [17]). We include it here for completeness and classify all separable solutions.
Assume that the solution to (12) has the form With primes denoting the derivatives with respect to the appropriate variables, Equation (12) can be written as 2 1 , for c ∈  . As we will see below, when equipped with the appropriate boundary conditions described next, this system has solutions for all real values of the constant c.
Recall that the boundary condition on the stream function is and velocity (and vorticity) components can be obtained using (5), (9), and (11).
Note that if we define a function ( )   In all displayed contour plots in this paper, the horizontal axis is the r-axis and the vertical axis is the z-axis. The thicker contours indicate where 0 ψ = .
The color coding (red vs. blue hues) distinguishes between the signs of the stream function ψ , and therefore between the regions in the rz-plane where the flow is clockwise or counterclockwise, and, in view of (9) and (11), it also corresponds to the tangential component of the flow being either into the rz-plane or out of it. Since the stream function can be arbitrarily rescaled due to the linearity of (12), we do not display any color bars indicating magnitudes of ψ . In all plots, we use

Spherical Coordinates
We next consider the spherical coordinates which have also been studied in the literature (some solutions appear, e.g. in [1]).
We again provide this case for completeness and again classify all separable solutions. In this coordinate system, Equation (12) Under the assumption of separability, for c ∈  . As we will see below in Lemma 1, when equipped with the appropri-Open Journal of Fluid Dynamics Note that if we define a function ( ) , we can rewrite the system as to be solved for 0 We now have the following lemma describing the spectrum of the eigenvalue problem for G (or g). Lemma 1. The system (13) with ( ) ( ) ( ) These limiting values will be 0 at the poles of the Γ function, i.e. when the ar-  for some m ∈  . In this case, the limits in (14) can be rewritten as  and the requirement of at least one of the two limits being 0 implies m ∈  .
Since negative values of m produce the same  (15) and it follows from the definition of 2 1  that they are all polynomials (see also [1]).

Paraboloidal Coordinates
The paraboloidal coordinates are The curves along which u and v are constant are shown in Figure 4(a). The curves with constant u are the parabolas opening up and the curves with constant v are the parabolas opening down. It is easy to see that Under the assumption of separability, for c ∈  . As we will see below, when equipped with the appropriate boundary conditions described next The general solution to (17) in the case 0 c ≠ can be written using the hyper- Remark. By using (13.2.4) and (13.4.1) in [25] and trigonometric identities, it can be shown that the solutions in (19) are real valued and, in fact, are given by ( ) Using (19), the stream function in paraboloidal coordinates has the form  (20) when 0 c = , this corresponds to ( ) 2 2 , sin sin , 2 2 which agrees with the special case solution given in (18).
The case with 0 c = is shown in Figure 5(a) and Figure 5

Prolate Spheroidal Coordinates
The prolate spheroidal coordinates are The curves along which u and v are constant are shown in Figure 4(b). The curves with constant u are the hyperbolas and the curves with constant v are the ellipses.
In order to obtain u and v from r and z, we can use the conversion formulas from [32], which, in this case, result in for c ∈  . As we will see below, when equipped with the appropriate boundary conditions described next, this system has solutions only for a discrete spectrum of values of the constant c. This spectrum is bounded below.
Recall that the boundary condition on the stream function is The z-axis consists of three intervals: ( ] the problem for f is underdetermined with a regular singular point at 0, and the Open Journal of Fluid Dynamics problem for g becomes a second-order boundary value eigenvalue problem with regular singular endpoints. The problem for g is a regular Sturm-Louville problem [33], and therefore there exists a countable, real, bounded below spectrum of values for c. Note that if we define functions ( ) , we can rewrite the system as Also note that in this case the two differential equations for F and G are actually the same, though solved on different intervals.
The easy case to solve analytically is when In the general case without the assumption 2 c A = , we have not found explicit solutions analytically and instead approximated them numerically. The idea is to first approximate the eigenvalues c by approximating the solution to the boundary-value problem for G, and once c is approximated, then approximate the solution to the initial-value problem for F. The problem for F can be supplied with a second initial condition ( ) 1 1 F ′ = , because other values simply lead to different scalings of F, and thus of ψ , which does not affect the stream surfaces of ψ . Contour plots for various values of the parameters are shown in Figures 6-9. In all of them 0.01 α = and 2 A ranges from 1 to 100. We note that the white lines visible along the r-axis for the odd eigenmodes are due to the contour plotter in Mathematica struggling with the piecewise function that converts r and z to v. The same is true in the next section.
We again note that for large µ the solution for F in (21) resembles that of 2 0 F A F ′′ + = and therefore exhibits near periodicity in µ . Consequently, the contour plots of the corresponding stream function ψ consist of repeating blocks along some of the curves shown in Figure 4(b) as indicated by the contour plots in Figures 6-9.

Oblate Spheroidal Coordinates
The oblate spheroidal coordinates are The curves along which u and v are constant are shown in Figure 4  In order to obtain u and v from r and z, we can again use the conversion formulas from [32], which now have the form     Under the assumption of separability, for c ∈  . As we will see below, when equipped with the appropriate boundary conditions described next, this system has solutions only for a discrete spectrum of values of the constant c. This spectrum is bounded below.
The boundary condition on the stream function is The positive z-axis corresponds to 2 v = π and the negative z-axis corresponds to 2 v = − π . Therefore, necessary boundary conditions are ( ) ( ) However, several more conditions have to be checked to ensure that ψ is continuous and differentiable in ( )

{ }
, : 0 r z r > . We first observe that the symmetry in the boundary-value eigenvalue problem for g implies that g is either an even or an odd function of v. The r-axis is divided into intervals ( ] ( ) when g is an even function. Hence the problem for f is underdetermined but with no singular points, and the problem for g again becomes a second-order boundary-value eigenvalue problem with regular singular endpoints, whose spectrum has the same properties as in the prolate spheriodal coordinate system case.
Note that if we define functions ( ) to be solved for 0 µ > and 1 1 η − < < , where the boundary conditions become ( ) ( ) when G is even since the symmetry of G is inherited from g. In this coordinate system, again we have not found any explicit solutions analytically and instead approximated them numerically in the way described in the previous section. The "missing" initial condition for F has been supplied by Contour plots for various values of the parameters are shown in Figures 10-13. In all of them 0.01 α = and 2 A ranges from 1 to 100. We again note that for large µ the solution for F in (22) resembles that of

Non-Stability with Respect to Axisymmetric Perturbations
We now address the important question of stability of the flows we have discovered in the previous sections. Since we are studying axisymmetric flows, we will focus on one type of stability, a stability with respect to axisymmetric perturbations. Open Journal of Fluid Dynamics  perturbations [34]. Conversely, there may be an instability if the derivative is negative somewhere in the flow. Here, w is the azimuthal component of the flow as defined in (5).
Recall from (9) and (11)   where C is the swirl, ψ is the stream function, and α is the abnormality con-

Tornado-Like Flows
We noted in Section 4 that the choice of the coordinate systems used in this paper was made with the intention of modeling tornado-like flows. Specifically, focusing on the corner flow near the origin, we can now address the question whether flows similar to those shown in Figure 1 are possible with Beltrami flows.
Note that the flows that correspond to the even eigenmodes for the spherical and both cases of the spheroidal coordinates have the r-axis (or, more accurately the 0 z = plane) as a stream surface, and therefore the part of the flow where 0 z > can be taken to model a flow above the (horizontal) ground. We also note that similar flows can be easily created from the odd eigenmodes by applying transformation 3. on the list of transformations discussed in the Symmetry Transformations subsection of Section 0 with any ζ ∈  , similar to what we showed above in Figure 5(c) in the context of paraboloidal coordinates.
In Figure 14 we show, from left to right, stream functions that correspond to the fourth eigenmodes of the prolate spheroidal coordinates (Figure 14(a) and Figure 14(b)), the spherical coordinates ( Figure 14(c)), and the oblate spheroidal coordinates (Figure 14(d) and Figure 14(e)). These can be thought of as snapshots of a continuous transformation in which a in the prolate spheroidal coordinates decreases to 0, the coordinates becoming spherical coordinates, and then a in the oblate spheroidal coordinates increases from 0. One can visualize this transformation by looking at Figure 4  Finally, the last two images, (d) and (e), correspond to the flow in which the central downdraft reaches the ground and results in a horizontal outflow near the ground. Therefore, the progression shown in Figure 14 can be viewed as a quasi-static model for the transition between the flows shown in Figure 1(b) and Figure 1(c).

Conclusions
In fluid dynamics, Beltrami flows have been, among other purposes, used for software validation and hypothesized to occur in supercell and tornadic flows. Consequently, it would be beneficial to have a rich catalog of such flows. In this paper we have attempted to construct such flows both analytically and numerically by focusing on Trkalian flows in several orthogonal coordinate systems. Motivated by tornado-like flows shown in Figure 1, we focused on incompressible, steady, axisymmetric flows which allowed us to use a stream function formulation. After some simplifying assumptions on the flows, we were able to construct solutions to the linear Bragg-Hawthorne Equation (12) in several suitable coordinate systems by reducing the problem to a system of two ordinary differential equations. In the coordinate systems, for which some solutions exist in the literature, specifically, the cylindrical and spherical coordinate systems, we have provided a complete set of separable solutions. In the paraboloidal coordinate system, we found all possible separable analytic solutions. In the case of the two spheroidal coordinate systems, we could not find analytic solutions, but we established the existence of countably many solutions, and we numerically computed the first few of them. The obtained solutions have been visualized using contour plots of the stream functions in the rz-plane and some of these solutions have been compared to two-cell tornadic flows. Additionally, we proposed ways to generate infinitely many new stream functions that can be constructed from the existing ones by continuously varying a scalar parameter. The richness of the solution set obtained in this paper with a constant abnormality α indicates that many other solutions can be found by exploring other three-dimensional coordinate systems and by allowing α to vary in space.