On the Axisymmetric Steady Incompressible Beltrami Flows

Pavel Bělík¹, Xueqing Su², Douglas P. Dokken³, Kurt Scholz³, Mikhail M. Shvartsman³

¹Department of Mathematics, Statistics, and Computer Science, Augsburg University, Minneapolis, MN, USA.
²Department of Statistics, University of Michigan, Ann Arbor, MI, USA.
³Department of Mathematics, University of St. Thomas, Saint Paul, MN, USA.

DOI: 10.4236/ojfd.2020.103014                

Abstract

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 tornado-like flows that develop as swirl increases. The results in this paper are equally applicable in electromagnetics, where the equivalent concept is that of a force-free magnetic field.

Keywords

Axisymmetric Beltrami Flow, Trkalian Flow, Bragg-Hawthorne Equation, Cylindrical Coordinates, Spherical Coordinates, Paraboloidal Coordinates, Prolate Spheroidal Coordinates, Oblate Spheroidal Coordinates, Vorticity

1. 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 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]. Some Beltrami flows in cylindrical and spherical coordinates are described in, for example [1] [20] [21].

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 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].

2. Beltrami Flows

Inviscid fluid flows are usually modeled by the Euler Equation (plus relevant energy and state equations)

$\frac{\partial u}{\partial t}+\left(u\cdot \nabla \right)u=-\frac{\nabla p}{\rho }-gz\mathrm{,}$ (1)

where $u\left(x\mathrm{,}t\right)$ is the velocity of the fluid at some point $x\in {ℝ}^3$ and time $t\in ℝ$, $p\left(x\mathrm{,}t\right)$ is the pressure, $\rho \left(x\mathrm{,}t\right)$ 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 $x\in {ℝ}^3$. The mass conservation equation is

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho u\right)=\frac{\partial \rho }{\partial t}+u\cdot \nabla \rho +\rho \left(\nabla \cdot u\right)=0.$ (2)

Taking “curl” of the equation for velocity (1), we obtain an equation for the vorticity, $\omega =\nabla \times u$,

$\frac{\partial \omega }{\partial t}=\nabla \times \left(u\times \omega \right)+\frac{1}{{\rho }^2}\left(\nabla \rho \times \nabla p\right),$ (3)

also known as the “vorticity” equation. A Beltrami flow is one in which the vorticity and velocity satisfy the Beltrami condition [5]

$\omega =\alpha u\text{for}\text{some}\alpha \ne \mathrm{0,}$ (4)

where $\alpha$ is called the abnormality and is, in general, a function of position and time. This condition implies that $u\times \omega =0$ everywhere.

It is shown in [5] that a steady Beltrami pair $u$, $\omega$ for which $\nabla \cdot u=0$ is a solution to the Euler Equation (1) with constant density $\rho$ ; 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 $\nabla \rho \times \nabla p$ will vanish. Also, taking divergence of both sides of (4) and using $\nabla \cdot u=0$ results in a necessary condition on the function $\alpha$

$u\cdot \nabla \alpha =\mathrm{0,}$

so $\alpha$ is constant on the streamlines of the flow.

A steady Beltrami pair $u$, $\omega$ for which $\nabla \cdot u=0$ is a solution to the governing equations with nonconstant (steady) mass density as well. From (1) we can solve for $\frac{\nabla p}{\rho }$, and (2) implies $u\cdot \nabla \rho =0$. Thus any smooth enough mass density field that is constant on the streamlines will satisfy the governing equations as well.

Finally, even if $u$ is not divergence free, (1) still provides a solution for $\frac{\nabla p}{\rho }$, and if $\rho$ can be found that satisfies the mass Equation (2), then a complete solution to the governing equations is obtained. The mass density would need to satisfy the ordinary differential equation

$\frac{\text{d}\rho }{\text{d}s}=-\rho \left(\nabla \cdot u\right)$

along each streamline, parameterized with the parameter s.

The special case of (4), in which $\alpha$ is independent of the spatial variables, is known as a Trkalian flow [28] [29]. In this case, $\alpha$ is also independent of time [29]. Our work focuses exclusively on Trkalian flows.

3. 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 assumptions of incompressibility, axisymmetry, and steady state allow us to use the stream function formulation in cylindrical coordinates $\left(r\mathrm{,}\theta \mathrm{,}z\right)$ with a stream function $\psi \left(r\mathrm{,}z\right)$. The corresponding velocity of the flow is then

$u=-\frac{1}{r}\frac{\partial \psi }{\partial z}r+w\theta +\frac{1}{r}\frac{\partial \psi }{\partial r}z,$ (5)

where $r$, $\theta$, and $z$ are the cylindrical coordinates basis vectors, and $w\left(r\mathrm{,}z\right)$ is the tangential component of the velocity. In this stream function formulation the incompressibility condition, $\nabla \cdot u=0$, is automatically satisfied provided $\psi$ is smooth enough. The vorticity, $\omega =\nabla \times u$, then satisfies

$\omega =-\frac{1}{r}\frac{\partial \left(rw\right)}{\partial z}r-\frac{1}{r}{D}^2\psi \theta +\frac{1}{r}\frac{\partial \left(rw\right)}{\partial r}z\mathrm{,}$ (6)

where

$D^2\psi =r\frac{\partial }{\partial r}\left(\frac{1}{r}\frac{\partial \psi }{\partial r}\right)+\frac{{\partial }^2\psi }{\partial z^2}=\frac{{\partial }^2\psi }{\partial r^2}-\frac{1}{r}\frac{\partial \psi }{\partial r}+\frac{{\partial }^2\psi }{\partial z^2}.$ (7)

Notice the similarity with the laplacian of $\psi$ which in the axisymmetric case is

$\Delta \psi =\frac{{\partial }^2\psi }{\partial r^2}+\frac{1}{r}\frac{\partial \psi }{\partial r}+\frac{{\partial }^2\psi }{\partial z^2}.$ (8)

We let

$C=rw\text{and}H=\frac{p}{\rho }+\frac{1}{2}{\left|u\right|}^2$ (9)

be the swirl (sometimes referred to as angular momentum) and energy head, respectively, where $p\left(r\mathrm{,}z\right)$ is the fluid pressure and $\rho$ its mass density. The conditions for a steady flow are [30]

$C=C\left(\psi \right),H=H\left(\psi \right),$

i.e. both the swirl and the head are constant on the stream surfaces $\psi =\text{const}\text{.}$, and hence along streamlines. The momentum Equation (1) can be written as

$D^2\psi =r^2\frac{\text{d}H}{\text{d}\psi }-C\frac{\text{d}C}{\text{d}\psi }\mathrm{,}$ (10)

which is known as the Bragg-Hawthorne equation.

When the energy head is constant in space so that $\text{d}H/\text{d}\psi =0$, Equation (10) reduces to $D^2\psi =-C\frac{\text{d}C}{\text{d}\psi }$ and $-\frac{1}{r}{D}^2\psi =w\frac{\text{d}C}{\text{d}\psi }$. The vorticity (6) can then be rewritten as

$\omega =\frac{\text{d}C}{\text{d}\psi }\left(-\frac{1}{r}\frac{\partial \psi }{\partial z}r+w\theta +\frac{1}{r}\frac{\partial \psi }{\partial r}z\right)=\frac{\text{d}C}{\text{d}\psi }u\mathrm{,}$

so the corresponding flow is a Beltrami flow, in which velocity and vorticity are parallel everywhere in space.

In what follows, we will focus on the simple case

$C\left(\psi \right)=\alpha \psi$ (11)

studied, for example, in [31]. In this case, $\omega =\alpha 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

$D^2\psi =-{\alpha }^2\psi ,$ (12)

where the operator $D^2$ is defined in (7). Note the similarity to the well-studied and well-understood Helmholtz equation, in which the differential operator on the left is the laplacian (8). Due to the assumption of axisymmetry, the stream function $\psi \left(r\mathrm{,}z\right)$ solving (12) is sought in the half plane $\left\{\left(r,z\right):r>0,z\in ℝ\right\}$ with a constant boundary condition on the z-axis in order for the axis to be a streamline. Additionally, we set the constant to be zero, since otherwise $C\left(\psi \right)=\alpha \psi$ would imply nonzero swirl C on the z-axis, and the expression $C=rw$ would lead to an unbounded azimuthal component w of the velocity on the z-axis. The condition that $\psi =0$ for $r=0$ will be used in the next sections to specify boundary conditions on the functions arising by separating variables.

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 $\psi \left(r\mathrm{,}z\right)$ is a solution, then so are $\psi \left(r\mathrm{,}-z\right)$ and $\psi \left(r\mathrm{,}z+\zeta \right)$ for any $\zeta \in ℝ$. This follows since only a second derivative in z is present in the operator $D^2$.

3) If $\psi \left(r\mathrm{,}z\right)$ is a solution and $\psi \left(r,-z\right)=\psi \left(r,z\right)$, then $\Psi \left(r,z\right)=\psi \left(r,z-\zeta \right)-\psi \left(r,z+\zeta \right)$ is a solution that satisfies $\Psi \left(r\mathrm{,0}\right)=0$ and $\Psi \left(r,-z\right)=-\Psi \left(r,z\right)$ for any $\zeta \in ℝ$.

4) If $\psi \left(r\mathrm{,}z\right)$ is a solution and $\psi \left(r,-z\right)=-\psi \left(r,z\right)$, then $\Psi \left(r\mathrm{,}z\right)=\psi \left(r\mathrm{,}z-\zeta \right)+\psi \left(r\mathrm{,}z+\zeta \right)$ is a solution that satisfies $\Psi \left(r\mathrm{,0}\right)=0$ and $\Psi \left(r,-z\right)=-\Psi \left(r,z\right)$ for any $\zeta \in ℝ$.

5) If $\psi \left(r\mathrm{,}z\right)$ is a solution, then $\Psi \left(r\mathrm{,}z\right)=\psi \left(r\mathrm{,}\zeta +z\right)-\psi \left(r\mathrm{,}\zeta -z\right)$ is a solution that satisfies $\Psi \left(r,0\right)=0$ and $\Psi \left(r,-z\right)=-\Psi \left(r,z\right)$ for any $\zeta \in ℝ$.

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.

4. 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.

Motivated by the schematics of a progression of tornadic flows shown in Figure 1 and discussed in more detail, for example, in [6], in which a single-cell updraft flow (Figure 1(a)) becomes a two-cell flow with a central downdraft (Figure 1(b) and 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

Figure 1. Schematic examples of tornadic flows illustrating a typical progression with increasing swirl. A single-cell updraft flow in (a) bifurcates into two-cell flows in (b) and (c), and further bifurcates into several vortices in (d) as swirl increases.

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].

4.1. 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 $\psi \left(r,z\right)=f\left(r\right)g\left(z\right)$ with $r\ge 0$ and $-\infty <z<\infty$. With primes denoting the derivatives with respect to the appropriate variables, Equation (12) can be written as

$f^{″}g-\frac{1}{r}{f}^{\prime }g+f{g}^{″}=-{\alpha }^{2}fg,$

which can be separated into the equations

$f^{″}-\frac{1}{r}{f}^{\prime }+cf=0,$

$g^{″}+\left({\alpha }^2-c\right)g=0$

for $c\in ℝ$. 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 $\psi \left(\mathrm{0,}z\right)=f\left(0\right)g\left(z\right)=0$ which implies that $f\left(0\right)=0$. With this boundary condition $\psi$ is a ${\mathcal{C}}^{\infty }$ function on $\left\{\left(r,z\right):r>0\right\}$ and velocity (and vorticity) components can be obtained using (5), (9), and (11).

Note that if we define a function $F\left(\mu \right)$ with $\mu =r^2$ via $f\left(r\right)=F\left(\mu \right)$, we can rewrite the system as

$F^{″}+\frac{c}{4\mu }F=0,$

$g^{″}+\left({\alpha }^2-c\right)g=0,$

to be solved for $\mu >0$ and $-\infty <z<\infty$, where the boundary condition becomes $F\left(0\right)=0$.

4.1.1 Case ${\alpha }^2-c=0$

In this case we immediately have $g\left(z\right)=a+bz$ as a solution for $z\in ℝ$. The equation for f now has the general solution $f\left(r\right)=c_{1}r{J}_1\left(\left|\alpha \right|r\right)+c_{2}r{Y}_1\left(\left|\alpha \right|r\right)$, where $J_1$ an $Y_1$ are the Bessel functions of the first and second kind, respectively. The initial condition $f\left(0\right)=0$ forces $c_2=0$, and we have the solution for the stream function

$\psi \left(r\mathrm{,}z\right)=\left(a+bz\right)r{J}_1\left(\left|\alpha \right|r\right)\mathrm{.}$

4.1.2 Case ${\alpha }^2-c<0$

Note that in this case $c>0$ and we have the solution $f\left(r\right)=c_{1}r{J}_1\left(r\sqrt{c}\right)$ and $g\left(z\right)=a{\text{e}}^{z\sqrt{c-{\alpha }^2}}+b{\text{e}}^{-z\sqrt{c-{\alpha }^2}}$. The stream function is then given by

$\psi \left(r,z\right)=r{J}_1\left(r\sqrt{c}\right)\left(a{\text{e}}^{z\sqrt{c-{\alpha }^2}}+b{\text{e}}^{-z\sqrt{c-{\alpha }^2}}\right).$

4.1.3 Case ${\alpha }^2-c>0$

The general solution for g is $g\left(z\right)=a\sin \left(z\sqrt{{\alpha }^2-c}\right)+b\cos \left(z\sqrt{{\alpha }^2-c}\right)$ and the solution for f depends on the sign of c. If $c=0$, then $f\left(r\right)=c_1{r}^2$. If $c>0$, then, as above, $f\left(r\right)=c_{1}r{J}_1\left(r\sqrt{c}\right)$. Finally, if $c<0$, then $f\left(r\right)=c_{1}ir{J}_1\left(ir\sqrt{-c}\right)$, which is a real-valued function. The stream function can then take one of the three forms,

$\psi \left(r\mathrm{,}z\right)=\left(asin\left(z\sqrt{{\alpha }^2-c}\right)+bcos\left(z\sqrt{{\alpha }^2-c}\right)\right)\cdot (\begin{array}{l}r{J}_1\left(r\sqrt{c}\right)\text{for}c>0,\hfill \\ r^2\text{for}c=0,\hfill \\ ir{J}_1\left(ir\sqrt{-c}\right)\text{for}c<0.\hfill \end{array}$

It follows from the above results that there are three types of structures that these solutions have as shown in Figure 2. When $c\ge {\alpha }^2$, the contour plots of $\psi$ consist of vertical infinite or semi-infinite strips separated by vertical lines where $\psi =0$ (Figure 2(a)). When $0<c<{\alpha }^2$, the contour plots consist of rectangular blocks (Figure 2(b)). When $c\le 0$, the contour plots consist of horizontal semi-infinite strips (Figure 2(c)).

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 $\psi =0$. The color coding (red vs. blue hues) distinguishes between the signs of the stream function $\psi$, 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 $\psi$. In all plots, we use $\alpha =0.01$. This choice is affected by the observation that in some tornadic flows with velocity on the order of tens of meters per second the vorticity is on the order of tenths per second. The displayed window is motivated by a tornado scale and the length units can be thought of as meters.

Figure 2. Contour plots of the stream functions $\psi$ obtained using cylindrical coordinates. (a) Case 4.1.1 with $c={\alpha }^2$, qualitatively similar to Case 4.1.2 with $c>{\alpha }^2$ ; (b) Case 4.1.3 with $0<c<{\alpha }^2$ ; and (c) Case 4.1.3 with $c<0$, qualitatively similar to Case 4.1.3 with $c=0$.

4.2. Spherical Coordinates

We next consider the spherical coordinates

$r=R\sin \left(\varphi \right),$

$z=R\cos \left(\varphi \right),$

$R\ge 0,0\le \varphi \le \pi ,$

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) becomes

$\frac{{\partial }^2\psi }{\partial R^2}+\frac{1}{R^2}\frac{{\partial }^2\psi }{\partial {\varphi }^2}-\frac{\cot \left(\varphi \right)}{R^2}\frac{\partial \psi }{\partial \varphi }=-{\alpha }^2\psi .$

Under the assumption of separability, $\psi \left(R\mathrm{,}\varphi \right)=f\left(R\right)g\left(\varphi \right)$, this equation becomes

$f^{″}g+\frac{1}{R^2}f{g}^{″}-\frac{cot\left(\varphi \right)}{R^2}f{g}^{\prime }=-{\alpha }^{2}fg\mathrm{,}$

which can be separated into the equations

$R^2{f}^{″}+\left({\alpha }^2{R}^2-c\right)f=\mathrm{0,}$

$g^{″}-cot\left(\varphi \right)g^{\prime }+cg=0$

for $c\in ℝ$. As we will see below in Lemma 1, when equipped with the appropriate boundary conditions described next, this system has solutions only for a discrete spectrum of positive integer values of the constant c.

The boundary condition on the stream function, $\psi \left(\mathrm{0,}z\right)=0$, immediately implies that $g\left(0\right)=g\left(\pi \right)=0$. Additionally, if $f\left(0\right)\ne 0$ and $g\cancel{\equiv }0$, then $\psi$ is not continuous at the origin as can be seen by computing limits of $f\left(R\right)g\left(\varphi \right)$ as $R\to 0$ with various values of $\varphi$. Therefore, the required boundary conditions are

$f\left(0\right)=0\text{and}g\left(0\right)=g\left(\pi \right)=0.$

With these boundary conditions $\psi$ is a ${\mathcal{C}}^{\infty }$ function on $\left\{\left(r,z\right):r>0\right\}$. The problem for f is underdetermined with a regular singular point $R=0$, and for g we have a second-order boundary value eigenvalue problem with regular singular endpoints.

Note that if we define a function $G\left(\eta \right)$ with $\eta =\cos \left(\varphi \right)$ for $0\le \varphi \le \pi$ via $g\left(\varphi \right)=G\left(\eta \right)$, we can rewrite the system as

$f^{″}+\left({\alpha }^2-\frac{c}{R^2}\right)f=0,G^{″}+\frac{c}{1-{\eta }^2}G=0,$ (13)

to be solved for $R>0$ and $-1<\eta <1$, where the boundary conditions become

$f\left(0\right)=0\text{and}G\left(-1\right)=G\left(1\right)=0.$

We now have the following lemma describing the spectrum of the eigenvalue problem for G (or g).

Lemma 1. The system (13) with $\left(R\mathrm{,}\eta \right)\in \left(\mathrm{0,}\infty \right)\times \left(-\mathrm{1,1}\right)$ and the boundary conditions $f\left(0\right)=G\left(\pm 1\right)=0$ has nontrivial solutions if and only if $c=m\left(m+1\right)$ for $m\in \mathbb{N}$.

Proof. When $c=0$, the second equation in (13) has only the trivial solution, so let’s assume $c\ne 0$. The general solution can then be written in terms of the hypergeometric function ${}_2\mathcal{F}{}_1$

$G\left(\eta \right)=c_1\cdot {}_2\mathcal{F}{}_1\left(-\frac{d_1}{4}\mathrm{,}\frac{d_2}{4}\mathrm{,}\frac{1}{2}\mathrm{,}{\eta }^2\right)+c_2\eta \cdot {}_2\mathcal{F}{}_1\left(-\frac{d_2}{4}\mathrm{,}\frac{d_1}{4}\mathrm{,}\frac{3}{2}\mathrm{,}{\eta }^2\right)\mathrm{,}$

where $c_1\mathrm{,}{c}_2\in ℝ$ and $d_1=\sqrt{1+4c}+1$ and $d_2=\sqrt{1+4c}-1$. We now have

$\begin{array}{l}\underset{\eta \to 1}{lim}{}_2\mathcal{F}{}_1\left(-\frac{d_1}{4}\mathrm{,}\frac{d_2}{4}\mathrm{,}\frac{1}{2}\mathrm{,}{\eta }^2\right)=\frac{\sqrt{\pi }}{\Gamma \left(\frac{3-\sqrt{1+4c}}{4}\right)\Gamma \left(\frac{3+\sqrt{1+4c}}{4}\right)}\mathrm{,}\\ \underset{\eta \to 1}{lim}{}_2\mathcal{F}{}_1\left(-\frac{d_2}{4}\mathrm{,}\frac{d_1}{4}\mathrm{,}\frac{3}{2}\mathrm{,}{\eta }^2\right)=-\frac{2\sqrt{\pi }}{c\Gamma \left(\frac{1-\sqrt{1+4c}}{4}\right)\Gamma \left(\frac{1+\sqrt{1+4c}}{4}\right)}\mathrm{.}\end{array}$ (14)

These limiting values will be 0 at the poles of the $\Gamma$ function, i.e. when the arguments of $\Gamma$ are non-positive integers. In all other cases, the limiting values of the first fundamental solution at ±1 are some nonzero value a, while the limits of the second fundamental solution will be ±b for some nonzero b. No linear combination of these solutions will satisfy both boundary conditions. When $c<-1/4$, all of the arguments in the $\Gamma$ functions are non-real, so we can assume $c\ge -1/4$. Any such c can be written as a product $c=m\left(m+1\right)$ for some $m\in ℝ$. In this case, the limits in (14) can be rewritten as

$\underset{\eta \to 1}{lim}{}_2\mathcal{F}{}_1\left(-\frac{m+1}{2}\mathrm{,}\frac{m}{2}\mathrm{,}\frac{1}{2}\mathrm{,}{\eta }^2\right)=\frac{\sqrt{\pi }}{\Gamma \left(\frac{1-m}{2}\right)\Gamma \left(\frac{m+2}{2}\right)}\mathrm{,}$

$\underset{\eta \to 1}{lim}{}_2\mathcal{F}{}_1\left(-\frac{m}{2}\mathrm{,}\frac{m+1}{2}\mathrm{,}\frac{3}{2}\mathrm{,}{\eta }^2\right)=-\frac{2\sqrt{\pi }}{m\left(m+1\right)\Gamma \left(-\frac{m}{2}\right)\Gamma \left(\frac{m+1}{2}\right)}\mathrm{,}$

and the requirement of at least one of the two limits being 0 implies $m\in \mathbb{Z}$. Since negative values of m produce the same $c=m\left(m+1\right)$ as nonnegative ones, and since $m=0$ corresponds to $c=0$, we have that only $m\in \mathbb{N}$ need to be considered to produce nontrivial solutions for G. These solutions are, up to a multiplicative constant,

$G\left(\eta \right)=(\begin{array}{ll}{}_2\mathcal{F}{}_1\left(-\frac{m+1}{2}\mathrm{,}\frac{m}{2}\mathrm{,}\frac{1}{2}\mathrm{,}{\eta }^2\right)\hfill & \text{for}m\text{odd,}\hfill \\ {}_2\mathcal{F}{}_1\left(-\frac{m}{2}\mathrm{,}\frac{m+1}{2}\mathrm{,}\frac{3}{2}\mathrm{,}{\eta }^2\right)\cdot \eta \hfill & \text{for}m\text{even}\mathrm{,}\hfill \end{array}$ (15)

and it follows from the definition of ${}_2\mathcal{F}{}_1$ that they are all polynomials (see also [1]).

For any $m\in \mathbb{N}$ and $c=m\left(m+1\right)$, the equation for f in (13) together with its boundary condition has a nontrivial solution in terms of the Bessel function of the first kind, again up to a multiplicative constant,

$f\left(R\right)=\sqrt{R}{J}_{\beta }\left(\alpha R\right)\text{with}\beta =\frac{2m+1}{2}\mathrm{.}$ (16)

□

Using (15) and (16), the stream function in spherical coordinates has the form

$\psi \left(R\mathrm{,}\varphi \right)=\sqrt{R}{J}_{\beta }\left(\alpha R\right)\cdot (\begin{array}{ll}{}_2\mathcal{F}{}_1\left(-\frac{m+1}{2}\mathrm{,}\frac{m}{2}\mathrm{,}\frac{1}{2}\mathrm{,}co{s}^2\left(\varphi \right)\right)\hfill & \text{for}m\text{odd,}\hfill \\ {}_2\mathcal{F}{}_1\left(-\frac{m}{2}\mathrm{,}\frac{m+1}{2}\mathrm{,}\frac{3}{2}\mathrm{,}co{s}^2\left(\varphi \right)\right)cos\left(\varphi \right)\hfill & \text{for}m\text{even}\mathrm{,}\hfill \end{array}$

where $m\in \mathbb{N}$ and $\beta =\left(2m+1\right)/2$. In the original cylindrical coordinates $\left(r\mathrm{,}z\right)$, we have

$\psi \left(r\mathrm{,}z\right)=(\begin{array}{l}{\left(r^2+z^2\right)}^{1/4}{J}_{\beta }\left(\alpha \sqrt{r^2+z^2}\right)\cdot {}_2\mathcal{F}{}_1\left(-\frac{m+1}{2}\mathrm{,}\frac{m}{2}\mathrm{,}\frac{1}{2}\mathrm{,}\frac{z^2}{r^2+z^2}\right)\text{for}m\text{odd}\mathrm{,}\hfill \\ z{\left(r^2+z^2\right)}^{-1/4}{J}_{\beta }\left(\alpha \sqrt{r^2+z^2}\right)\cdot {}_2\mathcal{F}{}_1\left(-\frac{m}{2}\mathrm{,}\frac{m+1}{2}\mathrm{,}\frac{3}{2}\mathrm{,}\frac{z^2}{r^2+z^2}\right)\text{for}m\text{even}\text{.}\hfill \end{array}$

Contour plots of some of these stream functions have been shown in literature (see, e.g. [1]), but we include them in Figure 3 for completeness and comparison to solutions in other coordinate systems.

Figure 3. Contour plots of the stream functions $\psi$ obtained using spherical coordinates. Results for the first six eigenmodes ( $m=1,\cdots ,6$ ) are shown left to right and top to bottom.

4.3. Paraboloidal Coordinates

The paraboloidal coordinates are

$r=uv,$

$z=\frac{1}{2}\left(v^2-u^2\right),$

$u\ge \mathrm{0,}v\ge 0.$

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

$u^2=\sqrt{r^2+z^2}-z,v^2=\sqrt{r^2+z^2}+z.$

In this coordinate system, Equation (12) becomes

$\frac{1}{u^2+v^2}\left(\frac{{\partial }^2\psi }{\partial u^2}+\frac{{\partial }^2\psi }{\partial v^2}-\frac{1}{u}\frac{\partial \psi }{\partial u}-\frac{1}{v}\frac{\partial \psi }{\partial v}\right)=-{\alpha }^2\psi .$

Under the assumption of separability, $\psi \left(u\mathrm{,}v\right)=f\left(u\right)g\left(v\right)$, becomes

$f^{″}g+f{g}^{″}-\frac{1}{u}{f}^{\prime }g-\frac{1}{v}f{g}^{\prime }=-{\alpha }^2\left(u^2+v^2\right)fg,$

which can be separated into the equations

$f^{″}-\frac{1}{u}{f}^{\prime }+\left({\alpha }^2{u}^2-c\right)f=\mathrm{0,}$

$g^{″}-\frac{1}{v}{g}^{\prime }+\left({\alpha }^2{v}^2+c\right)g=0$

for $c\in ℝ$. 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.

The boundary condition on the stream function, $\psi \left(\mathrm{0,}z\right)=0$, immediately implies

$f\left(0\right)=0\text{and}g\left(0\right)=0.$

with these boundary conditions $\psi$ is a ${\mathcal{C}}^{\infty }$ function on $\left\{\left(r,z\right):r>0\right\}$. Hence

Figure 4. Visualization of the (a) paraboloidal, (b) prolate spheroidal, and (c) oblate spheroidal coordinate systems in the rz-plane, the last two both with $a=1$. In the case of the paraboloidal coordinates, parabolas opening up correspond to u constant and parabolas opening down to v constant. For both spheroidal systems, ellipses correspond to u constant and hyperbolas to v constant.

both problems for f and g are underdetermined, and both problems have a regular singular point at 0.

Note that if we define functions $F\left(\mu \right)$ and $G\left(\eta \right)$ with $\mu =u^2$ and $\eta =v^2$ via $f\left(u\right)=F\left(\mu \right)$ and $g\left(v\right)=G\left(\eta \right)$, we can rewrite the system as

$F^{″}+\frac{{\alpha }^2\mu -c}{4\mu }F=0,G^{″}+\frac{{\alpha }^2\eta +c}{4\eta }G=0,$ (17)

to be solved for $\mu >0$ and $\eta >0$, where the boundary conditions become

$F\left(0\right)=0\text{and}G\left(0\right)=0.$

The special case when $c=0$ immediately results in

$F\left(\mu \right)=c_{1}sin\left(\frac{\alpha \mu }{2}\right)\mathrm{,}G\left(\eta \right)=c_{2}sin\left(\frac{\alpha \eta }{2}\right)\mathrm{.}$ (18)

The general solution to (17) in the case $c\ne 0$ can be written using the hypergeometric functions ${}_1\mathcal{F}{}_1$ and U

$F\left(\mu \right)=\mu {\text{e}}^{-\frac{i\left|\alpha \right|\mu }{2}}\left[{}_1\mathcal{F}{}_1\left(\frac{4i\left|\alpha \right|+c}{4i\left|\alpha \right|}\mathrm{,2,}i\left|\alpha \right|\mu \right)c_1+U\left(\frac{4i\left|\alpha \right|+c}{4i\left|\alpha \right|}\mathrm{,2,}i\left|\alpha \right|\mu \right)c_3\right]$

$G\left(\eta \right)=\eta {\text{e}}^{-\frac{i\left|\alpha \right|\eta }{2}}\left[{}_1\mathcal{F}{}_1\left(\frac{4i\left|\alpha \right|-c}{4i\left|\alpha \right|}\mathrm{,2,}i\left|\alpha \right|\eta \right)c_2+U\left(\frac{4i\left|\alpha \right|-c}{4i\left|\alpha \right|}\mathrm{,2,}i\left|\alpha \right|\eta \right)c_4\right]$

where $c_i\in ℝ$ for $i=1,2,3,4$. We now have

$\underset{\mu \to 0}{lim}F\left(\mu \right)=\frac{4c_3}{c\Gamma \left(\frac{c}{4i\left|\alpha \right|}\right)}\mathrm{,}\underset{\eta \to 0}{lim}G\left(\eta \right)=-\frac{4c_4}{c\Gamma \left(-\frac{c}{4i\left|\alpha \right|}\right)}\mathrm{,}$

and therefore to satisfy the boundary conditions we need to have $c_3=c_4=0$. The solution for $c\ne 0$ is then

$\begin{array}{l}F\left(\mu \right)=c_1\mu {\text{e}}^{-\frac{i\left|\alpha \right|\mu }{2}}{}_1\mathcal{F}{}_1\left(\frac{4i\left|\alpha \right|+c}{4i\left|\alpha \right|}\mathrm{,2,}i\left|\alpha \right|\mu \right)\mathrm{,}{c}_1\in ℝ\mathrm{,}\\ G\left(\eta \right)=c_2\eta {\text{e}}^{-\frac{i\left|\alpha \right|\eta }{2}}{}_1\mathcal{F}{}_1\left(\frac{4i\left|\alpha \right|-c}{4i\left|\alpha \right|}\mathrm{,2,}i\left|\alpha \right|\eta \right)\mathrm{,}{c}_2\in ℝ\mathrm{.}\end{array}$ (19)

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

$F\left(\mu \right)=C_1\mu {\displaystyle {\int }_0^1}cos\left[\frac{c}{4\left|\alpha \right|}ln\left(\frac{1-u}{u}\right)+\left|\alpha \right|\mu \left(u-\frac{1}{2}\right)\right]\text{d}u\mathrm{,}$

$G\left(\eta \right)=C_2\eta {\displaystyle {\int }_0^1}cos\left[-\frac{c}{4\left|\alpha \right|}ln\left(\frac{1-u}{u}\right)+\left|\alpha \right|\eta \left(u-\frac{1}{2}\right)\right]\text{d}u\mathrm{,}$

where $C_i=c_i\frac{4\left|\alpha \right|}{c\pi }sinh\left(\frac{c\pi }{4\left|\alpha \right|}\right)$ for $i=1,2$.

Using (19), the stream function in paraboloidal coordinates has the form

$\begin{array}{c}\psi \left(u\mathrm{,}v\right)=u^2{\text{e}}^{-\frac{i\left|\alpha \right|u^2}{2}}\cdot {}_1\mathcal{F}{}_1\left(1+\frac{c}{4i\left|\alpha \right|}\mathrm{,2,}i\left|\alpha \right|u^2\right)\\ \times v^2{\text{e}}^{-\frac{i\left|\alpha \right|v^2}{2}}\cdot {}_1\mathcal{F}{}_1\left(1-\frac{c}{4i\left|\alpha \right|}\mathrm{,2,}i\left|\alpha \right|v^2\right)\mathrm{.}\end{array}$ (20)

when $c=0$, this corresponds to

$\psi \left(u\mathrm{,}v\right)=sin\left(\frac{\alpha u^2}{2}\right)sin\left(\frac{\alpha v^2}{2}\right)\mathrm{,}$

which agrees with the special case solution given in (18).

The case with $c=0$ is shown in Figure 5(a) and Figure 5(b) corresponds to $c=-0.05$. It is clear from (20) and the definition of the paraboloidal coordinates that changing the sign of c corresponds to interchanging the roles of u and v and therefore changing the sign of z. Consequently, a contour plot for $c=0.05$ can be obtained from that with $c=-0.05$. Since the flows with $c\ne 0$ do not appear to have any kind of symmetry with respect to the $z=0$ plane, they can be used, together with transformation 5. From the Symmetry Transformations subsection of Section 0 to generate new flows for which $z=0$ is a stream surface. This is illustrated in Figure 5(c) in which (20) with $c=-0.05$ and transformation 5. with $\zeta =0$ were used.

Finally, we note that for large $\mu$ the solution for F in (17) resembles that of $F^{″}+\frac{{\alpha }^2}{4}F=0$ and therefore exhibits near periodicity in $\mu$. Consequently, the

Figure 5. Contour plots of stream functions using paraboloidal coordinates. The stream function (20) with (a) $c=0$ and (b) $c=-0.05$ ; (c) the stream function $\Psi$ obtained from $\psi$ in (20) with $c=-0.05$ by applying transformation 5. from the Symmetry Transformations subsection of Section 0 with $\zeta =0$.

contour plots of the corresponding stream function $\psi$ consist of repeating blocks along some of the curves shown in Figure 4(a) as indicated by the contour plots in Figure 5.

4.4. Prolate Spheroidal Coordinates

The prolate spheroidal coordinates are

$r=asinh\left(u\right)sin\left(v\right)\mathrm{,}$

$z=acosh\left(u\right)cos\left(v\right)\mathrm{,}$

$u\ge \mathrm{0,}0\le v\le \pi \mathrm{.}$

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

$u=\frac{1}{2}ln\left(1-2q+2\sqrt{q^2-q}\right)\mathrm{,}v=(\begin{array}{ll}arcsin\left(\sqrt{p}\right)\hfill & \text{if}z\ge \mathrm{0,}\hfill \\ \pi -arcsin\left(\sqrt{p}\right)\hfill & \text{if}z<\mathrm{0,}\hfill \end{array}$

where

$p=\frac{-B+\sqrt{B^2+4a^2{r}^2}}{2a^2},q=\frac{-B-\sqrt{B^2+4a^2{r}^2}}{2a^2},B=r^2+z^2-a^2.$

In this coordinate system, Equation (12) becomes

$\frac{1}{a^2\left(cos{h}^2\left(u\right)-co{s}^2\left(v\right)\right)}\left(\frac{{\partial }^2\psi }{\partial u^2}+\frac{{\partial }^2\psi }{\partial v^2}-coth\left(u\right)\frac{\partial \psi }{\partial u}-cot\left(v\right)\frac{\partial \psi }{\partial v}\right)=-{\alpha }^2\psi \mathrm{.}$

Under the assumption of separability, $\psi \left(u\mathrm{,}v\right)=f\left(u\right)g\left(v\right)$, becomes

$f^{″}g+f{g}^{″}-coth\left(u\right)f^{\prime }g-cot\left(v\right)f{g}^{\prime }=-{\alpha }^2{a}^2\left(cos{h}^2\left(u\right)-co{s}^2\left(v\right)\right)fg\mathrm{,}$

which, with $A=\alpha a$, can be separated into the equations

$f^{″}-coth\left(u\right)f^{\prime }+\left(A^{2}cos{h}^2\left(u\right)-c\right)f=\mathrm{0,}$

$g^{″}-cot\left(v\right)g^{\prime }-\left(A^{2}co{s}^2\left(v\right)-c\right)g=0$

for $c\in ℝ$. 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 $\psi \left(\mathrm{0,}z\right)=0$. The z-axis consists of three intervals: $\left(-\infty \mathrm{,}-a\right]$ corresponding to $v=\pi$, $\left[-a\mathrm{,}a\right]$ corresponding to $u=0$, and $\left[a\mathrm{,}\infty \right)$ corresponding to $v=0$. Therefore, necessary boundary conditions for f and g are

$f\left(0\right)=0\text{and}g\left(0\right)=g\left(\pi \right)=0.$

with these boundary conditions $\psi$ is a ${\mathcal{C}}^{\infty }$ function on $\left\{\left(r,z\right):r>0\right\}$. Hence the problem for f is underdetermined with a regular singular point at 0, and the 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 $F\left(\mu \right)$ and $G\left(\eta \right)$ with $\mu =cosh\left(u\right)$ and $\eta =cos\left(v\right)$ via $f\left(u\right)=F\left(\mu \right)$ and $g\left(v\right)=G\left(\eta \right)$, we can rewrite the system as

$F^{″}+\frac{A^2{\mu }^2-c}{{\mu }^2-1}F=0,G^{″}+\frac{c-A^2{\eta }^2}{1-{\eta }^2}G=0,$ (21)

to be solved for $\mu >1$ and $-1<\eta <1$, where the boundary conditions become

$F\left(1\right)=0\text{and}G\left(-1\right)=G\left(1\right)=0.$

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 $c=A^2$, since in this case the differential equations reduce to $F^{″}+A^{2}F=0$ and $G^{″}+A^{2}G=0$. The equation for G, together with its boundary conditions, then forces $A=n\pi /2$ for $n\in \mathbb{Z}$, and we obtain two sets of solutions,

$F\left(\mu \right)=c_{1}cos\left(\frac{\left(2n+1\right)\pi }{2}\mu \right)\mathrm{,}G\left(\nu \right)=c_{2}cos\left(\frac{\left(2n+1\right)\pi }{2}\nu \right)\text{for}n\in \mathbb{Z}\mathrm{,}$

or

$F\left(\mu \right)=c_{1}sin\left(n\pi \mu \right)\mathrm{,}G\left(\nu \right)=c_{2}sin\left(n\pi \nu \right)\text{for}n\in \mathbb{Z}\mathrm{.}$

In the general case without the assumption $c=A^2$, 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 $F^{\prime }\left(1\right)=1$, because other values simply lead to different scalings of F, and thus of $\psi$, which does not affect the stream surfaces of $\psi$. Contour plots for various values of the parameters are shown in Figures 6-9. In all of them $\alpha =0.01$ and $A^2$ 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 $\mu$ the solution for F in (21) resembles that of $F^{″}+A^{2}F=0$ and therefore exhibits near periodicity in $\mu$. Consequently, the contour plots of the corresponding stream function $\psi$ consist of repeating blocks along some of the curves shown in Figure 4(b) as indicated by the contour plots in Figures 6-9.

4.5. Oblate Spheroidal Coordinates

The oblate spheroidal coordinates are

$r=a\cosh \left(u\right)\cos \left(v\right),$

Figure 6. Contour plots of the stream functions $\psi$ obtained using prolate spheroidal coordinates. Results for the first six eigenmodes with $A^2=1$ are shown left to right and top to bottom.

$z=a\sinh \left(u\right)\sin \left(v\right),$

$u\ge \mathrm{0,}-\pi /2\le v\le \pi /2\mathrm{.}$

The curves along which u and v are constant are shown in Figure 4(c). The curves with constant u are the ellipses and the curves with constant v are the hyperbolas.

Figure 7. Contour plots of the stream functions $\psi$ obtained using prolate spheroidal coordinates. Results for the first six eigenmodes with $A^2=10$ are shown left to right and top to bottom.

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

$u=\frac{1}{2}ln\left(1-2q+2\sqrt{q^2-q}\right)\mathrm{,}v=(\begin{array}{ll}arcsin\left(\sqrt{p}\right)\hfill & \text{if}z\ge \mathrm{0,}\hfill \\ -arcsin\left(\sqrt{p}\right)\hfill & \text{if}z<\mathrm{0,}\hfill \end{array}$

Figure 8. Contour plots of the stream functions $\psi$ obtained using prolate spheroidal coordinates. Results for the first six eigenmodes with $A^2=50$ are shown left to right and top to bottom.

where

$p=\frac{-B+\sqrt{B^2+4a^2{z}^2}}{2a^2},q=\frac{-B-\sqrt{B^2+4a^2{z}^2}}{2a^2},B=r^2+z^2-a^2.$

In this coordinate system, Equation (12) becomes

$\frac{1}{a^2\left(sin{h}^2\left(u\right)+si{n}^2\left(v\right)\right)}\left(\frac{{\partial }^2\psi }{\partial u^2}+\frac{{\partial }^2\psi }{\partial v^2}-tanh\left(u\right)\frac{\partial \psi }{\partial u}+tan\left(v\right)\frac{\partial \psi }{\partial v}\right)=-{\alpha }^2\psi \mathrm{.}$

Figure 9. Contour plots of the stream functions $\psi$ obtained using prolate spheroidal coordinates. Results for the first six eigenmodes with $A^2=100$ are shown left to right and top to bottom.

Under the assumption of separability, $\psi \left(u\mathrm{,}v\right)=f\left(u\right)g\left(v\right)$, this equation becomes

$f^{″}g+f{g}^{″}-tanh\left(u\right)f^{\prime }g+tan\left(v\right)f{g}^{\prime }=-{\alpha }^2{a}^2\left(sin{h}^2\left(u\right)+si{n}^2\left(v\right)\right)fg\mathrm{,}$

which, with $A=\alpha a$, can be separated into the equations

$f^{″}-tanh\left(u\right)f^{\prime }+\left(A^{2}sin{h}^2\left(u\right)-c\right)f=\mathrm{0,}$

$g^{″}+tan\left(v\right)g^{\prime }+\left(A^{2}si{n}^2\left(v\right)+c\right)g=\mathrm{0,}$

for $c\in ℝ$. 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 $\psi \left(\mathrm{0,}z\right)=0$. The positive z-axis corresponds to $v=\pi /2$ and the negative z-axis corresponds to $v=-\pi /2$. Therefore, necessary boundary conditions are $g\left(-\pi /2\right)=g\left(\pi /2\right)=0$.

However, several more conditions have to be checked to ensure that $\psi$ is continuous and differentiable in $\left\{\left(r,z\right):r>0\right\}$. 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 $\left(\mathrm{0,}a\right]$ corresponding to $u=0$ and $\left[a\mathrm{,}\infty \right)$ corresponding to $v=0$. To ensure continuity of $\psi$ in $\left\{\left(r,z\right):r>0\right\}$, only continuity across the segment $\left(\mathrm{0,}a\right]$ needs to be addressed. If g is odd, there is a discontinuity in $\psi$ unless $f\left(0\right)=0$. If g is even, $\psi$ is continuous across the segment without a restriction on $f\left(0\right)$.

It can be verified that $\psi$ is differentiable in the z-direction across the segment $\left(\mathrm{0,}a\right]$ on the r-axis, so no new requirements arise there. However, in the case of an even g, the stream function $\psi$ is differentiable in the r-direction at the point $\left(a\mathrm{,0}\right)$ on the r-axis only when $f^{\prime }\left(0\right)=0$.

Therefore, the boundary conditions are $g\left(-\pi /2\right)=g\left(\pi /2\right)=0$ and either $f\left(0\right)=0$ when g is an odd function, or $f^{\prime }\left(0\right)=0$ 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 $F\left(\mu \right)$ and $G\left(\eta \right)$ with $\mu =\sinh u$ and $\eta =\sin v$ via $f\left(u\right)=F\left(\mu \right)$ and $g\left(v\right)=G\left(\eta \right)$, we can rewrite the system as

$F^{″}+\frac{A^2{\mu }^2-c}{1+{\mu }^2}F=0,G^{″}+\frac{A^2{\eta }^2+c}{1-{\eta }^2}G=0,$ (22)

to be solved for $\mu >0$ and $-1<\eta <1$, where the boundary conditions become $G\left(-1\right)=G\left(1\right)=0$, and $F\left(0\right)=0$ when G is odd, and $F^{\prime }\left(0\right)=0$ 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 $F^{\prime }\left(0\right)=1$ when $F\left(0\right)=0$ and by $F\left(0\right)=1$ when $F^{\prime }\left(0\right)=0$. Contour plots for various values of the parameters are shown in Figures 10-13. In all of them $\alpha =0.01$ and $A^2$ ranges from 1 to 100.

We again note that for large $\mu$ the solution for F in (22) resembles that of $F^{″}+A^{2}F=0$ and therefore exhibits near periodicity in $\mu$. Consequently, the contour plots of the corresponding stream function $\psi$ consist of repeating blocks along some of the curves shown in Figure 4(c) as indicated by the contour plots in Figures 10-13.

Figure 10. Contour plots of the stream functions $\psi$ obtained using oblate spheroidal coordinates. Results for the first six eigenmodes with $A^2=1$ are shown left to right and top to bottom.

4.6. 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.

Figure 11. Contour plots of the stream functions $\psi$ obtained using oblate spheroidal coordinates. Results for the first six eigenmodes with $A^2=10$ are shown left to right and top to bottom.

The Rayleigh criterion of stability with respect to axisymmetric perturbations states that if

$\frac{\partial {\left(rw\right)}^2}{\partial r}>0$

everywhere in the flow, then the flow is stable with respect to axisymmetric

Figure 12. Contour plots of the stream functions $\psi$ obtained using oblate spheroidal coordinates. Results for the first six eigenmodes with $A^2=50$ are shown left to right and top to bottom.

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) that in our considerations

${\left(rw\right)}^2=C^2={\alpha }^2{\psi }^2,$

Figure 13. Contour plots of the stream functions $\psi$ obtained using oblate spheroidal coordinates. Results for the first six eigenmodes with $A^2=100$ are shown left to right and top to bottom.

where C is the swirl, $\psi$ is the stream function, and $\alpha$ is the abnormality constant. Noting that for most of the solutions displayed above the level set $\psi =0$ divides the rz-plane into bounded regions in which ${\psi }^2>0$, it is clear the such solutions are unstable with respect to axisymmetric perturbations. They typically have $\partial {\left(rw\right)}^2/\partial r>0$ for smaller values of r in each such region and $\partial {\left(rw\right)}^2/\partial r<0$ for larger values of r.

5. 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 $z=0$ plane) as a stream surface, and therefore the part of the flow where $z>0$ 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 $\zeta \in ℝ$, 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(b), allowing the two focal points $\left(\mathrm{0,}\pm a\right)$ merge into one at the origin as $a\to 0$, and moving the merged point along the r-axis to transition to the right panel in the figure. In Figure 14, one can interpret the left two images, (a) and (b), as a two-cell vortex with a horizontal inflow and a vertical updraft near the corner with a central downdraft near the z-axis. The middle image, (c), representing the flow in spherical coordinates, corresponds to the flow in which the stagnation point reaches the ground.

Figure 14. Illustrations of flows that have the characteristics of a tornadic flow with increasing swirl. Prolate spheroidal coordinates with (a) $A^2=10$ and (b) $A^2=1$ ; (c) spherical coordinates; oblate spheroidal coordinates with (d) $A^2=1$ and (e) $A^2=10$. All plots correspond to $\alpha =0.01$ and to the fourth eigenmode.

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).

6. 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 $\alpha$ indicates that many other solutions can be found by exploring other three-dimensional coordinate systems and by allowing $\alpha$ to vary in space.

Acknowledgements

Sincere thanks to the Office of Undergraduate Research and Graduate Opportunity (URGO) at Augsburg University and to Dean and Amy Sundquist for providing the research opportunity and financial support for Ms. Xueqing Su. The authors would also like to thank Mr. Pierre-Henri Chavanis for his useful feedback and suggestions.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Marsh, G.E. (1996) Force-Free Magnetic Fields: Solutions, Topology and Applications. World Scientific Publishing Co. Pte Ltd., Singapore. https://doi.org/10.1142/2965
[2]	Dritschel, D.G. (1991) Generalized Helical Beltrami Flows in Hydrodynamics and Magnetohydrodynamics. Journal of Fluid Mechanics, 222, 525-541. https://doi.org/10.1017/S0022112091001209
[3]	Wang, C.Y. (1991) Exact Solutions of the Steady-State Navier-Stokes Equations. Annual Review Fluid Mechanics, 23, 159-177. https://doi.org/10.1146/annurev.fl.23.010191.001111
[4]	Constantin, P. and Majda, A. (1988) The Beltrami Spectrum for Incompressible Fluid Flows. Communications in Mathematical Physics, 115, 435-456. https://doi.org/10.1007/BF01218019
[5]	Majda, A.J. and Bertozzi, A. (2001) Vorticity and Incompressible Flows. Cambridge Texts in Applied Mathematics. Cambridge University Press, Cambridge.
[6]	Bluestein, H.B. (2013) Severe Convective Storms and Tornadoes. Observations and Dynamics. Springer, Berlin. https://doi.org/10.1007/978-3-642-05381-8
[7]	Noda, A.T. and Niino, H. (2010) A Numerical Investigation of a Super-Cell Tornado: Genesis and Vorticity Budget. Journal of the Meteorological Society of Japan, 88, 135-159. https://doi.org/10.2151/jmsj.2010-203
[8]	Weisman, M.L. and Rotunno, R. (2000) The Use of Vertical Wind Shear versus Helicity in Interpreting Supercell Dynamics. Journal of the Atmospheric Sciences, 57, 1452-1472. https://doi.org/10.1175/1520-0469(2000)057%3C1452:TUOVWS%3E2.0.CO;2
[9]	Sasaki, Y.K. (2014) Entropic Balance Theory and Variational Field Lagrangian Formalism: Tornadogenesis. Journal of the Atmospheric Sciences, 71, 2104-2113. https://doi.org/10.1175/JAS-D-13-0211.1
[10]	Davies-Jones, R.P. (2008) Can a Descending Rain Curtain in a Supercell Instigate Tornadogenesis Barotropically? Journal of the Atmospheric Sciences, 65, 2469-2497. https://doi.org/10.1175/2007JAS2516.1
[11]	Lilly, D.K. (1983) Dynamics of Rotating Thunderstorms. In: Lilly, D.K. and Gal-Chen, T., Eds., Mesoscale Meteorology—Theories, Observations and Models, NATO ASI Series (Series C: Mathematical and Physical Sciences), Vol. 114, Springer, Dordrecht, 531-544. https://doi.org/10.1007/978-94-017-2241-4_28
[12]	Lilly, D.K. (1986) The Structure, Energetics and Propagation of Rotating Convective Storms. Part II: Helicity and Storm Stabilization. Journal of the Atmospheric Sciences, 43, 126-140. https://doi.org/10.1175/1520-0469(1986)043%3C0126:TSEAPO%3E2.0.CO;2
[13]	Davies-Jones, R. and Richardson, Y.P. (2002) An Exact Anelastic Beltrami-Flow Solution for Use in Model Validation. Preprints, 19th Conference on Weather Analysis and Forecasting/15th Conference on Numerical Weather Prediction, San Antonio, TX, 11-16 August 2002, 43-46.
[14]	Shapiro, A. (1993) The Use of an Exact Solution of the Navier-Stokes Equations in a Validation Test of a Three-Dimensional Nonhydrostatic Numerical Model. Monthly Weather Review, 121, 2420-2425. https://doi.org/10.1175/1520-0493(1993)121%3C2420:TUOAES%3E2.0.CO;2
[15]	Skamarock, W.C., Klemp, J.B., Dudhia, J., Gill, D.O., Liu, Z., Berner, J., Wang, W., Powers, J.G., Duda, M.G., Barker, D.M. and Huang, X.Y. (2019) A Description of the Advanced Research WRF Version 4. NCAR, Tech. Rep., Tech. Note. https://opensky.ucar.edu/islandora/object/opensky:2898
[16]	Ahmad, N.N. and Proctor, F.H. (2011) Simulation of Benchmark Cases with the Terminal Area Simulation System (TASS). 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 4-7 January 2011, Orlando, Florida. https://doi.org/10.2514/6.2011-1005
[17]	Leprovost, N., Dubrulle, B. and Chavanis, P.-H. (2006) Dynamics and Thermodynamics of Axisymmetric Flows: Theory. Physical Review E, 73, Article ID: 046308. https://doi.org/10.1103/PhysRevE.73.046308
[18]	Monchaux, R., Ravelet, F., Dubrulle, B., Chiffaudel, A. and Daviaud, F. (2006) Properties of Steady States in Turbulent Axisymmetric Flows. Physical Review Letters, 96, Article ID: 124502. https://link.aps.org/doi/10.1103/PhysRevLett.96.124502 https://doi.org/10.1103/PhysRevLett.96.124502
[19]	Naso, A., Monchaux, R., Chavanis, P.-H. and Dubrulle, B. (2010) Statistical Mechanics of Beltrami Flows in Axisymmetric Geometry: Theory Reexamined. Physical Review E, 81, Article ID: 066318. https://doi.org/10.1103/PhysRevE.81.066318
[20]	Davies-Jones, R. (1985) Dynamical Interaction between an Isolated Convective Cell and a Veering Environmental Wind. 14th Conference on Severe Local Storms, Indianapolis, IN, American Meteorological Society, Boston, 28 October-1 November 1985, 216-219.
[21]	Marsh, G.E. (1992) Axially Symmetric Solutions to the Force-Free Magnetic-Field Equations in Spherical and Cylindrical Coordinates. Physical Review A, 45, 7520-7525. https://doi.org/10.1103/PhysRevA.45.7520
[22]	Elcrat, A.R., Fornberg, B. and Miller, K.G. (2008) Steady Axisymmetric Vortex Flows with Swirl and Shear. Journal of Fluid Mechanics, 613, 395-410. https://doi.org/10.1017/S002211200800342X
[23]	Wang, C.Y. (1966) On a Class of Exact Solutions of the Navier-Stokes Equations. Journal of Applied Mechanics, 33, 696-698. https://doi.org/10.1115/1.3625151
[24]	Dyck, N.J. and Straatman, A.G. (2020) Exact Solutions to the Three-Dimensional Navier-Stokes Equations Using the Extended Beltrami Method. Journal of Applied Mechanics, 87, Article ID: 011004. https://doi.org/10.1115/1.4044927
[25]	Olver, F.W.J., Lozier, D.W., Boisvert, R.F. and Clark, C.W. (2010) NIST Handbook of Mathematical Functions. Cambridge University Press, Cambridge. https://dlmf.nist.gov/
[26]	Davies-Jones, R., Trapp, R.J. and Bluestein, H.B. (2001) Tornadoes and Tornadic Storms. In: Doswell, C.A., Ed., Severe Convective Storms, American Meteorological Society, Boston, 167-221. https://doi.org/10.1007/978-1-935704-06-5_5
[27]	Church, C.R., Snow, J.T. and Agee, E.M. (1977) Tornado Vortex Simulation at Purdue University. Bulletin of the American Meteorological Society, 58, 900-908. https://doi.org/10.1175/1520-0477(1977)058%3C0900:TVSAPU%3E2.0.CO;2
[28]	Lakhtakia, A. (1994) Viktor Trkal, Beltrami Fields, and Trkalian Flows. Czechoslovak Journal of Physics, 44, 89-96. https://doi.org/10.1007/BF01701185
[29]	Trkal, V. (1994) A Note on the Hydrodynamics of Viscous Fluids. Czechoslovak Journal of Physics, 44, 97-106. https://doi.org/10.1007/BF01701186
[30]	Batchelor, G.K. (2000) An Introduction to Fluid Dynamics. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511800955
[31]	Moffatt, H.K. (1969) The Degree of Knottedness of Tangled Vortex Lines. Journal of Fluid Mechanics, 35, 117-129. https://doi.org/10.1017/S0022112069000991
[32]	Sun, C. (2017) Explicit Equations to Transform from Cartesian to Elliptic Coordinates. Mathematical Modelling and Applications, 2, 43-46. https://doi.org/10.11648/j.mma.20170204.12
[33]	Stakgold, I. and Holst, M. (2011) Green’s Functions and Boundary Value Problems. John Wiley & Sons, Inc., Hoboken. https://doi.org/10.1002/9780470906538
[34]	Drazin, P.G. (2002) Introduction to Hydrodynamic Stability. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511809064