1 Introduction↩︎

Tornadoes are among the most extreme and destructive weather phenomena and, due to the numerous casualties and substantial damage in property, their modeling has attracted considerable interest. At the core of mathematical models for the description of tornadoes is the concept of swirling flows. Assuming the tornado structure does not change significantly in time and that a vortex line resembles the tornado core, [1], a common model for tornado-like flows is the stationary incompressible axisymmetric Navier-Stokes equations. Namely, the following system in cylindrical coordinates is considered \[\label{axi-sys} \begin{align} u \frac{\partial u}{\partial r} + w \frac{\partial u}{\partial z} - \frac{v^2}{r} & = \nu \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial u}{\partial r}\Big)+ \frac{\partial^2 u}{\partial z^2} - \frac{u}{r^2} \Big] - \frac{\partial p}{\partial r}, \\ u \frac{\partial v}{\partial r} + w \frac{\partial v}{\partial z} + \frac{uv}{r} & = \nu \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial v}{\partial r}\Big)+ \frac{\partial^2 v}{\partial z^2} - \frac{v}{r^2} \Big], \\ u \frac{\partial w}{\partial r} + w \frac{\partial w}{\partial z} \qquad & = \nu \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial w}{\partial r}\Big)+ \frac{\partial^2 w}{\partial z^2} \qquad \Big] - \frac{\partial p}{\partial z}, \\ \frac{1}{r} \frac{\partial }{\partial r} (ru) + \frac{\partial w}{\partial z} & = 0, \end{align}\tag{1}\] where \(\vec{u}=(u,v,w): \mathbb{R}^3\times \mathbb{R}_+ \to \mathbb{R}^3\) is the velocity, \(p : \mathbb{R}^3 \times \mathbb{R}_+\to \mathbb{R}\) is the pressure and \(\nu \geq 0\) is the coefficient of kinematic viscosity. Modeling tornadoes necessitates the consideration of certain additional factors, like rotation induced from the cloud and buoyancy effects. The reader is referred to App. 9 for a quick presentation of models and references. Nevertheless, system 1 is considered as the base-model for describing their core structure, [2].

Long [3], [4] introduced a self-similar ansatz and reduced the stationary axisymmetric Navier-Stokes equations to a system of ordinary differential equations focusing on the boundary layer description. Independently, Goldshtik [5] examined the existence of self-similar solutions as the Reynolds number varies; he concluded his system is not solvable for all Reynolds numbers and characterized this loss of existence as a ‘paradox’. Few years later, Serrin [6] presented another self-similar ansatz and showed that only three types of solution can be associated with the interaction of an infinite vortex line with a boundary plane: the flow can be either ascending, or descending, or a combination of these two profiles, i.e. downward near the vortex line and inward near the r-axis. The latter is usually referred as a double-celled vortex and Serrin [6] presented extensive comparisons of such solutions to tornadoes. Many authors subsequently studied the system 1 on occasion considering more general families of self-similar solutions, for example [7]–[9], or more general geometries including for example conical flows, [9]–[14].

This work has two objectives: First, to consider the full system 1 and provide a complete study of the stationary self-similar solutions of axisymmetric Navier-Stokes and of axisymmetric Euler equations. Second, to compare those equations and study the emergence of Euler solutions from the Navier-Stokes solutions as the viscosity \(\nu\) goes to zero. One novel feature for the Euler equations is the consideration of solutions with a slip discontinuity in the velocities. Regarding the Navier-Stokes equations, we obtain an alternative derivation of the equations studied by Serrin [6] and study the emergence of stationary self-similar solutions of the Euler equations as zero-viscosity limits from corresponding solutions of the Navier-Stokes equations. Our analysis indicates that first there are no solutions with slip-discontinuity for the stationary, self-similar Euler equations, and second that the solutions of Navier-Stokes approach the single cell solution of the Euler equations in the zero-viscosity limit. Numerical calculations indicate the presence of double-celled solutions but this happens for finite (even relatively large viscosities) and does not persist as the viscosity decreases.

We now provide an outline of the work. We introduce the self-similar transformation \[\begin{align} u(r,z) = \frac{1}{r} U(\xi), \quad v(r,z) = \frac{1}{r} V(\xi), \quad w(r,z) =\frac{1}{r} W(\xi),\quad p(r,z) =\frac{1}{r^2} P(\xi) \, , \end{align}\] in \(\xi = \frac{z}{r}\) to 1 which captures the scale-invariance of Navier-Stokes equations. This leads, for the axisymmetric Euler equations, to a system of ordinary differential equations: \[\label{Eul1} \begin{align} \bigg[\frac{\theta^2}{2} + (1+\xi^2)P \bigg]' &= -\xi V^2, \\ V' \theta &= 0, \\ \bigg[\theta^2 - \xi \Big(\frac{\theta^2}{2}\Big)' + P \bigg]' &= 0, \end{align}\tag{2}\] for \((\theta, V, P)\), where \(\theta\) is a stream-function with \(U = - \theta^\prime\) and \(W = \theta - \xi \theta^\prime\). The same transformation leads for the axisymmetric Navier-Stokes equations to the system \[\label{NS1} \begin{align} \bigg[\frac{\theta^2}{2} + (1+\xi^2)P \bigg]' &= \nu \bigg[\xi\theta - (1+\xi^2) \theta' \bigg]' -\xi V^2, \\ V' \theta &= \nu \Big[3\xi V' + (1+\xi^2) V''\Big], \\ \bigg[\theta^2 - \xi \Big(\frac{\theta^2}{2}\Big)' + P \bigg]' &= \nu \Big[\xi\theta - \xi^2 \theta' - \xi (1+\xi^2) \theta'' \Big]' \, , \end{align}\tag{3}\] augmented with the physically relevant boundary conditions (see section 2 for the derivations).

First, we consider 2 . Imposing no-penetration conditions on the boundary and the vortex core we derive explicit solutions that correspond to either ascending or descending flows. To examine whether a double-celled vortex may occur for the Euler equations, we consider the existence of solutions with slip discontinuity in velocity. These would emerge from the Rankine-Hugoniot conditions for 2 , \[\label{intro:RH} \llbracket P \rrbracket = 0, \qquad \llbracket \theta V \rrbracket = 0, \qquad \llbracket \theta \, \theta' \rrbracket = 0, \qquad\tag{4}\] where \(\llbracket \cdot \rrbracket\) is the jump operator. We prove that there do not exist solutions of the self-similar Euler equations with discontinuities at a finite number of points that satisfy the jump conditions 4 . This nonexistence result is then extended to a class of flows where there is mass input or loss through the vortex line. For the aforementioned flows, the explicit continuous solutions are also derived. The non-existence result leads to conjecture that the double-celled vortices cannot persist in the limit as the viscosity goes to zero.

Next, we focus on the stationary self-similar axisymmetric Navier-Stokes system 3 in the context of a vortex core interacting with a boundary. This problem is reduced to a coupled integro-differential system for \((\theta(\xi), V(\xi))\) \[\label{NS2} \begin{align} \frac{\theta^2}{2} - \nu \bigg[(1+\xi^2) \theta' + \xi \theta \bigg] &= \mathcal{G}\left(\xi\right), \\ \nu (1+\xi^2) V'' + \Big(3\nu\xi - \theta\Big) V' &= 0, \end{align}\tag{5}\] with boundary conditions \[\theta = \theta' = V = 0 \,\, \textrm{ at } \; \xi = 0, \qquad V \to V_\infty , \,\, U \to 0 \, , \, \, \textrm{ as} \; \xi \to \infty.\] where the functional \(\mathcal{G}\) depending on \(V(\cdot)\) and \(\xi\) is defined via \[\label{intro-defG} \mathcal{G}\left(\xi\right) = \mathcal{G}\left(V(\cdot), \xi\right) = \xi \sqrt{1+\xi^2} \int_\xi^\infty \frac{1}{\zeta^2(1+\zeta^2)^\frac{3}{2}} \bigg( \int_0^\zeta s V^2(s) ds\bigg) d\zeta + \,\, {E_0} \phi(\xi).\tag{6}\] The system is recast into two equivalent integrodifferential formulations and provides a common framework encompassing the pioneering works on this problem of Goldshitk [5], Serrin [6] and Goldstick-Shtern [9]. Existence of solutions for 5 is provided under the conditions described in [6]. Moreover the system 5 is solved numerically using an iterative solver and a numerical bifurcation diagram is presented, see Figure 11.

We then focus on the zero-viscosity limit from Navier-Stokes to Euler in the self-similar stationary setting. This leads to study the zero-viscosity limit \(\nu \to 0\) for the integrodifferential system 5 . The problem is conveniently recasted as an autonomous system, \[\label{intro-ns} \begin{align} \nu \, \frac{d \Theta_\nu}{dx} &= \tfrac{1}{2} \Theta_\nu^2 - \mathcal{F}\left(V_\nu ; x\right), \\ \nu \frac{d^2 {V_\nu}}{dx^2} &= \Theta_\nu \, \frac{d {V_\nu}}{dx} , \end{align}\tag{7}\] where \(\mathcal{F}\left(V_\nu ; x\right)\) is the analog of 6 (and is defined in 105 see section 5.2). An analysis of the possible configurations of solutions of 7 (based on ideas of Goldstick [5]) leads to a-priori estimates for solutions and classifying the potential shapes of the stream function to three configurations. To investigate the zero-viscosity limit, we employ compactness methods and exploit techniques developed in the theory of zero-viscosity limits for Riemann problems of conservation laws, [15], [16]. The derivative \(dV\) is viewed as a limit of a family of probability measures capturing the averaging process as \(\nu\to 0\). We pursue two theories: the first is based on \(L^p\) estimates for the stream function and variation estimates for the velocity. It uses weak convergence methods and leads to the convergence Theorem 4. The latter is not fully satisfactory due to the weak bounds available for the stream function. In a second step, under more restrictive conditions on the data for the problem, we obtain uniform variation estimates and invoke Helly’s theorem to conclude almost everywhere convergence (along subsequences) from the Navier-Stokes solutions to the Euler solutions. The convergence results are stated in Theorem 5, Proposition 8, and Corollary 3. They permit to characterize admissibility restrictions for the stationary Euler solutions emerging from Navier-Stokes in this setting.

Finally, we study the boundary layer in the case of a model problem for 7 that captures the essential behavior of our system. We use methods of asymptotic analysis, namely inner and outer solutions and matched asymptotic expansions. We deduce that the boundary layer is of order \(\nu^{2/3}\) and provide an asymptotic description of the stream function \(\Theta_\nu\) (see 145 ) and of \(V_\nu\) (see 150 ).

The work is organized as follows: In Section 2, we introduce the self-similar ansatz and derive the systems 2 and 3 that are used in the rest of this work. Sections 3 and 4 are devoted to stationary self-similar axisymmetric Euler equations. In section 3, an explicit continuous solution of 2 is derived, while in section 4 we investigate the existence of discontinuous solutions with finite number of discontinuities. In section 4, Euler solutions are also studied for flows with input of mass through the vortex line. In Section 5, we study self-similar stationary solutions for axisymmetric Navier-Stokes. After presenting some of their properties, we examine their limiting behavior as \(\nu\to0\) using compactness methods and analytical ideas from the theory of Borel measures. The results are stated in Theorem 4, Theorem 5 and Corollary 3. In Section 6, there is an asymptotic analysis of the boundary layer for a model problem, where inner and outer expansions leading to an explicit form of an asymptotic solution. In Section 7, a numerical scheme is presented for solving the stationary, self-similar system 7 . We present indicative profiles of typical flows that appear for various values of the parameters, and calculate a bifurcation of solutions in terms of the dimensionless parameters \((\frac{E_0}{V_\infty} , \frac{\nu}{V_\infty})\). The diagram is computed by solving the system 7 numerically and characterizing solutions according to the zone they belong, see Sections 5.3 and 7.4.

In Appendix 8 we list Navier-Stokes in cylindrical coordinates. As our study is motivated from the study of tornadoes, we present for the convenience of the reader in Appendix 9 a quick review of models used in the literature for tornado modeling.

2 Preliminaries↩︎

We consider the equations of motion for an incompressible homogeneous viscous fluid, with constant density \(\rho=1\) : \[\label{NS} \begin{align} \vec{u}_t + (\vec{u} \cdot \nabla) \vec{u} &= - \nabla p + \nu \, \Delta \vec{u} , \\ \nabla\cdot \vec{u} & = 0, \end{align}\tag{8}\] where \(\vec{u} : \mathbb{R}^3\times \mathbb{R}_+ \to \mathbb{R}^3\) is the velocity vector of the fluid, \(p : \mathbb{R}^3 \times \mathbb{R}_+\to \mathbb{R}\) is the pressure and \(\nu \geq 0\) is the coefficient of kinematic viscosity. The first equation represents the conservation of momentum while the second is the incompressibility condition and may be interpreted as describing conservation of mass. For \(\nu>0\), the system 8 is the so-called Navier-Stokes equations, while when viscocity effects are omitted the corresponding system 8 with \(\nu =0\) is the (incompressible) Euler equations.

2.1 Axisymmetric Navier - Stokes equations↩︎

In three space dimensions, introducing cylindrical coordinates \((r,\vartheta,z)\): \(x_1 = r \,\cos\vartheta\), \(x_2 = r\,\sin\vartheta\), \(x_3 = z\), and the attached unit vector system \[\vec{e}_r = (\cos\vartheta,\sin\vartheta,0), \,\; \vec{e}_\vartheta = (-\sin\vartheta,\cos\vartheta,0), \,\;\vec{e}_z = (0, 0, 1) \, ,\] we express the velocity vector \(\vec{u}\) in cylindrical coordinates as \[\vec{u} = u(r,\vartheta,z,t) \vec{e}_r + v(r,\vartheta,z,t) \vec{e}_{\vartheta} + w(r,\vartheta,z,t) \vec{e}_z.\] A flow is called axisymmetric if the velocity \(\vec{u}\) does not depend on the azimuthal angle \(\vartheta\), i.e. \[\vec{u} = u(r,z,t) \vec{e}_r + v(r,z,t) \vec{e}_{\vartheta} + w(r,z,t) \vec{e}_z.\] The velocity component \(v\) in the direction of \(\vec{e}_\vartheta\) is called the swirl velocity. If the swirl is everywhere equal to zero, i.e. \(v\equiv 0\), we name such flows as flows without swirl. Otherwise we call them axisymmetric flows with swirl.

Using that the velocity \(\vec{u}\) does not depend on \(\vartheta\), we derive the axisymmetric Navier-Stokes equations which have the following form (see App 8) \[\begin{align} \tag{9} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial r} + w \frac{\partial u}{\partial z} - \frac{v^2}{r} & = \nu \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial u}{\partial r}\Big)+ \frac{\partial^2 u}{\partial z^2} - \frac{u}{r^2} \Big] - \frac{\partial p}{\partial r}, \\ \tag{10} \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial r} + w \frac{\partial v}{\partial z} + \frac{uv}{r} & = \nu \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial v}{\partial r}\Big)+ \frac{\partial^2 v}{\partial z^2} - \frac{v}{r^2} \Big], \\ \tag{11} \frac{\partial w}{\partial t} + u \frac{\partial w}{\partial r} + w \frac{\partial w}{\partial z} \qquad & = \nu \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial w}{\partial r}\Big)+ \frac{\partial^2 w}{\partial z^2} \qquad \Big] - \frac{\partial p}{\partial z}, \\ \tag{12} \frac{1}{r} \frac{\partial }{\partial r} (ru) + \frac{\partial w}{\partial z} & = 0. \end{align}\]

For general axisymmetric flows, the vorticity vector \(\vec{\omega} = \nabla \times \vec{u}\) is expressed as \[\begin{align} \vec{\omega} &= \big({\omega}_r, {\omega}_{\vartheta}, {\omega}_z\big) = \bigg(-\frac{\partial v}{\partial z}, \frac{\partial u}{\partial z}-\frac{\partial w}{\partial r}, \, \frac{1}{r}\frac{\partial}{\partial r} (rv) \bigg) \, , \end{align}\] where \(\omega_\vartheta := \frac{\partial u}{\partial z}-\frac{\partial w}{\partial r}\) is the component of vorticity in the \(\vec{e}_\vartheta\)-direction. Note that for flows without swirl only the vorticity in the direction of \(\vec{e}_\vartheta\) survives: \(\vec{\omega} = \omega_\vartheta \vec{e}_\vartheta\).

The continuity equation 12 may be integrated using the axisymmetric stream function \(\psi(r,z,t)\) by expressing the velocity components \(u\) and \(w\) in terms of \(\psi\) as follows \[u = - \frac{1}{r} \frac{\partial \psi}{\partial z} \quad \text{and} \quad w= \frac{1}{r} \frac{\partial \psi}{\partial r}.\] The stream function \(\psi\) and the vorticity component \(\omega_\vartheta\) are independent of the swirl velocity \(v\); hence, an equivalent formulation of \(\eqref{tNS1} - \eqref{tNS4}\) is given by \[\begin{align} \tag{13} \frac{\partial {\omega}_{\vartheta}}{\partial t} + u \frac{\partial {\omega}_{\vartheta}}{\partial r} + w \frac{\partial {\omega}_{\vartheta}}{\partial z} - {\omega}_{\vartheta} \frac{u}{r} & = \nu \Big[\Delta {\omega}_{\vartheta} - \frac{ {\omega}_{\vartheta}}{r^2} \Big] + \frac{\partial }{\partial z} \Big(\frac{v^2}{r} \Big),\\ \tag{14} \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial r} + w \frac{\partial v}{\partial z} + v \frac{u}{r} & = \nu \Big[\Delta v - \frac{v}{r^2} \Big], \\ \tag{15} {\omega}_{\vartheta} & = - \frac{1}{r} \Delta \psi + \frac{2}{r^2} \frac{\partial \psi}{\partial r}, \\ \tag{16} \Big( u,w \Big) & = \frac{1}{r} \nabla^{\perp} \psi, \end{align}\] in terms of \(\omega_\theta\), \(v\) and \(\psi\), where \[\Delta f = \frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial f}{\partial r}\Big)+ \frac{\partial^2 f}{\partial z^2} \qquad \text{and} \qquad \nabla^{\perp} f = \Bigg(- \frac{\partial f}{\partial z} , \frac{\partial f}{\partial r} \Bigg).\] This formulation is equivalent with the three-dimensional Navier-Stokes equations, [17].

2.2 Self-Similar Formulation↩︎

In the analysis of partial differential equations it is often beneficial to seek solutions that conform with symmetry properties of the problem, for instance invariance under rotations, dilations, etc. The invariance of the Navier-Stokes equations under the scaling \[\vec{u}(t,r,z) = \lambda \, \vec{u}(\lambda^2 t,\lambda r, \lambda z) \quad \textrm{and} \quad p(t,r,z) = \lambda^2\, p(\lambda^2 t,\lambda r, \lambda z),\] for \(\lambda>0\), suggests to look for self-similar solutions of 9 12 that follow the ansatz \[\label{ansatz-time} \begin{align} & u(t,r,z) = \frac{1}{r} U(s,\xi), \quad v(t,r,z) = \frac{1}{r} V(s,\xi), \quad w(t,r,z) =\frac{1}{r} W(s,\xi) , \\[5pt] & p(t,r,z) =\frac{1}{r^2} P(s,\xi), \quad \psi(t,r,z) = r \Psi(s,\xi), \quad \omega_\vartheta(t,r,z) =\frac{1}{r^2} \Omega(s,\xi), \end{align}\tag{17}\] in the variables \(\xi = \frac{z}{r}\) and \(s= \frac{r}{\sqrt{t}}\). Such an ansatz induces a singularity at \(r=0\). In the literature related to tornadoes the singularity is considered to represent the line vortex resembling the tornado core, [1]. Using the ansatz 17 in 9 12 , a tedious but straightforward calculation yields the system of partial differential equations for \((U,V,W)\), \[\begin{align} s\Big(U-\frac{s^2}{2}\Big) U_s + \Big(W-\xi U\Big) U_{\xi} - \Big(U^2+V^2\Big) &= \nu \Big [ \mathcal{H}U - U \Big ] + 2P + \xi P_{\xi} - sP_s, \\ s\Big(U-\frac{s^2}{2}\Big) V_s + \Big(W-\xi U\Big) V_{\xi}&= \nu \Big[ \mathcal{H}V - V \Big], \\ s\Big(U-\frac{s^2}{2}\Big) W_s + \Big(W-\xi U\Big) W_{\xi}- \,\, UW &= \nu \Big[ \mathcal{H}W \qquad \Big] - P_{\xi}, \\ sU_s + W_{\xi} - \xi U_{\xi} &= 0 , \end{align}\] where the operator \(\mathcal{H}\) is defined as \[\label{defH} \mathcal{H}f = \bigg[ -\xi f + \Big((1+\xi^2)f\Big)_\xi \bigg]_\xi - sf_s - 2 s \xi f_{s\xi} + s^2 f_{ss}.\tag{18}\] The same ansatz transforms the velocity-vorticity equations \(\eqref{V1} - \eqref{V4}\) into the form \[\begin{align} s\Big(U-\frac{s^2}{2}\Big) \Omega_s + \big(W-\xi U\big) \Omega_{\xi} - 3U\Omega -2VV_{\xi} &= \nu \big[ \mathcal{D}\Omega - \Omega \big], \\ s\Big(U-\frac{s^2}{2}\Big) V_s + \big(W-\xi U\big) V_{\xi}&= \nu \big[ \mathcal{H}V - V \big], \\ \Psi -\xi \Psi_{\xi} - s\Psi_s + 2s\xi \Psi_{s\xi} - s^2 \Psi_{ss} - (1+\xi^2)\Psi_{\xi\xi} &= \quad \Omega, \\ \Big(- \Psi_{\xi} , \Psi - \xi \Psi_{\xi} + s \Psi_s \Big) & = \big( U,W \big), \end{align}\] where \(\mathcal{H}\) is defined in 18 and \(\mathcal{D}\) is given by \[\label{defD} \mathcal{D}f = 4f +5 \xi f_{\xi} - 3 sf_s - 2 s \xi f_{s\xi} + s^2 f_{ss} +(1+\xi^2)f_{\xi\xi}.\tag{19}\] The choice of ansatz is not unique and various self-similar transformations can be used for the Navier-Stokes equations, for instance, one may consider an ansatz in variables \(\xi = \frac{z}{r}\) and \(\tau= \frac{t}{r^2}\).

2.3 Stationary Self-Similar Formulation↩︎

The tornadic funnel typically moves slowly compared to the internal swirling velocities, [2]. This motivates to represent the core of the tornado via a vortex singularity, following [3], [5], and to study stationary axi-symmetric self-similar solutions for the Navier-Stokes equations. Letting \(\xi = \frac{z}{r}\), we seek for self-similar stationary solutions of 9 12 in the form \[\label{ansatz} \begin{align} & u(r,z) = \frac{1}{r} U(\xi), \quad v(r,z) = \frac{1}{r} V(\xi), \quad w(r,z) =\frac{1}{r} W(\xi) , \\ & p(r,z) =\frac{1}{r^2} P(\xi), \quad \psi(r,z) = r \Psi(\xi), \quad \omega_\vartheta(r,z) =\frac{1}{r^2} \Omega(\xi). \end{align}\tag{20}\] Such solutions satisfy the system of ordinary differential equations \[\begin{align} \tag{21} -U \Big(\xi U\Big)' + U' W - V^2 &= \nu \Big[\mathcal{L} U - U \Big] + 2P + \xi P' ,\\ -U \Big(\xi V\Big)' + V' W + UV &= \nu \Big[\mathcal{L} V - V \Big], \\ \tag{22} -U \Big(\xi W\Big)' + W' W &= \nu \mathcal{L} W - P', \\ W' - \xi U' &= 0 \end{align}\] where \(\xi \in \mathbb{R}_+\) and the operator \(\mathcal{L}\) defined as \[\label{defL} \mathcal{L}f = \bigg[ -\xi f + \Big((1+\xi^2)f\Big)' \bigg] ' \, .\tag{23}\] If we multiply \(\eqref{eq-1}\) by \(\xi\) and subtract it from 22 , we can rewrite the system as \[\begin{align*} \bigg[\frac{1}{2}\Big(W - \xi U\Big)^2 + (1+\xi^2)P \bigg]' &= \nu \big\{ \mathcal{L}W - \xi \mathcal{L}U + \xi U \big\} - \xi V^2 ,\\ V' \Big(W - \xi U\Big) &= \nu \big\{ \mathcal{L}V - V \big\}, \\ \bigg[W \Big(W - \xi U\Big) + P \bigg]' &= \nu \mathcal{L}W, \\ \Big(W - \xi U\Big)' &= -U . \end{align*}\]

For convenience, we introduce a new function \(\theta(\xi)\) with \(\theta = W -\xi U\). One checks from 16 and 20 that \(\theta (\xi) = \Psi (\xi)\) is the self-similar version of the stream function. Now, \(U\) and \(W\) are expressed in terms of \(\theta\) by \[\label{stream} U = - \theta' \, , \qquad W = \theta - \xi \theta^\prime \, .\tag{24}\] A tedious but straightforward computation renders [sss-ns] to a non-linear system of ordinary differential equations for \((\theta, V, P)\), \[\tag{25} \begin{align} \bigg[\frac{\theta^2}{2} + (1+\xi^2)P \bigg]' &= \nu \bigg[\xi\theta - (1+\xi^2) \theta' \bigg]' -\xi V^2, \\ V' \theta &= \nu \Big[3\xi V' + (1+\xi^2) V''\Big], \tag{26} \\ \bigg[\theta^2 - \xi \Big(\frac{\theta^2}{2}\Big)' + P \bigg]' &= \nu \Big[\xi\theta - \xi^2 \theta' - \xi (1+\xi^2) \theta'' \Big]', \end{align}\] where the viscosity \(\nu\) is a parameter and \((U,W)\) are determined by 24 . The aim is to study 25 in both the viscous \(\nu>0\) as well as the inviscid \(\nu=0\) cases, and to investigate the limiting relationship. Sections 3 and 4 are devoted to the inviscid case \(\nu=0\): In Section 3 an explicit solution for 25 (with \(\nu=0\)) is obtained. The existence of solutions with slip discontinuities is examined in Section 4. In Section 5, we express the system 25 for \(\nu>0\) into an equivalent integrodifferential formulation and study its limit as \(\nu\to0\). Numerical solutions are illustrated in Section 7.

3 A stationary self-similar solution for the axisymmetric Euler equations↩︎

Consider first the Euler equations 25 with \(\nu = 0\): \[\begin{align} \tag{27} \bigg[\frac{\theta^2}{2} + (1+\xi^2)P \bigg]' &= -\xi V^2 \\ \tag{28} V' \theta &= 0 \\ \tag{29} \bigg[\theta^2 - \xi \Big(\frac{\theta^2}{2}\Big)' + P \bigg]' &= 0\\ \tag{30} \theta' &= - U \end{align}\] Equation \(\eqref{v-eq}\) implies that if \(\theta(\xi) \ne 0\) then \(V(\xi)\) is a constant function. So, assume first that \(\theta \neq 0\) for \(\xi \in (0, \infty)\) and that \[\label{eqnV} V \equiv V_{0},\tag{31}\] where \(V_{0}\) is a constant proportional to the circulation around the axis. (The case where \(\theta\) vanishes and where \(V\) may exhibit discontinuities is examined in the next section.) Substituting 31 yields \[\begin{align} \label{sseuler-eq1} \xi (1+\xi^2) \Big(\frac{\theta^2}{2}\Big)' - (1+2 \xi^2) \Big(\frac{\theta^2}{2}\Big) &= -\frac{\xi^2}{2} {V_0}^2 - E_0 (1+ \xi^2)+ A_0 \\ P &= \xi \Big(\frac{\theta^2}{2}\Big)' - \theta^2 + E_0 \end{align}\tag{32}\] where \(E_{0}\) and \(A_{0}\) are integration constants. Dividing equation 32 by \(\xi^2 (1+\xi^2)^{\frac{3}{2}}\) gives \[\Bigg[ \frac{1}{\xi \sqrt{1+\xi^2}} \frac{\theta^2}{2} \Bigg] ' = -\frac{{V_0}^2}{2} \Bigg[ \frac{\xi}{\sqrt{1+\xi^2}} \Bigg] ' + {E_0} \Bigg[ \frac{\sqrt{1+\xi^2}}{\xi} \Bigg] ' - A_0 \Bigg[ \frac{1+2\xi^2}{\xi \sqrt{1+\xi^2}} \Bigg] '\] After an integration, we eventually obtain an explicit form for the stream function \(\theta\) and an explicit solution for the Euler system 27 - 30 \[\tag{33} \begin{align} \frac{\theta^2}{2} &= - \Big(\frac{V_0^2}{2} + 2{A_0} - {E_0}\Big) \xi^2 + {k_0} \xi \sqrt{1+\xi^2} + \Big({E_0} - {A_0}\Big) \tag{34} \\ U &= \frac{1}{\theta} \Bigg[ 2\Big(\frac{V_0^2}{2} + 2{A_0} - {E_0}\Big) \xi - k_0 \frac{1+2\xi^2}{\sqrt{1+\xi^2}} \Bigg] \tag{35} \\ W &= \frac{1}{\theta} \Bigg[k_0 \frac{\xi}{\sqrt{1+\xi^2}} + 2\Big({E_0} - {A_0}\Big)\Bigg] \\ P &= - k_0 \frac{\xi}{\sqrt{1+\xi^2}} + 2{A_0} - {E_0} \end{align}\] where \(V(\xi)\) is given by 31 and \(E_0\), \(A_{0}\) and \(k_0\) are integration constants.

We supplement system 27 30 with suitable boundary conditions. Recall that our objective is to model a vortex line that interacts with a boundary surface. We impose no-penetration boundary condition at the boundary, namely \(w = 0\) at \(z = 0\). Moreover, we assume that \(u = 0\) as \(r \to 0\) which ensures that no mass is added or lost through the vortex line. Expressing these conditions in terms of the self-similar functions 20 leads to \[W(0) = 0 \quad \textrm{and} \quad U(\xi) \to 0 \quad \textrm{as} \quad \xi \to \infty,\] or in terms of \(\theta\) \[\label{bc-euler} \theta(0) = 0 \quad \textrm{and} \quad \theta'(\xi) \to 0 \quad \textrm{as} \quad \xi \to \infty \, .\tag{36}\] The boundary conditions applied to 33 provide the relations \[\begin{align} \tag{37} {E_0} &= {A_0} \\ \tag{38} k_0 &= \frac{{V_0}^2}{2} + {E_0} \end{align}\] The derivation of 37 is direct. To derive relation 38 , note that 34 , 35 imply \[\begin{align} U = \frac{1}{\theta} \Bigg[ 2\Big(\frac{V_0^2}{2} + {E_0}\Big) \xi - k_0 \frac{1+2\xi^2}{\sqrt{1+\xi^2}} \Bigg] &=\frac{ 2\Big(\frac{V_0^2}{2} + {E_0}\Big) - k_0 \frac{1+2\xi^2}{\xi\sqrt{1+\xi^2}}}{\pm\sqrt{- 2\Big(\frac{V_0^2}{2} + {E_0}\Big) + 2{k_0}\frac{\sqrt{1+\xi^2}}{\xi}}} \\[8pt] &\xrightarrow{\xi\to\infty} \mp {\sqrt{- 2\Big(\frac{V_0^2}{2} + {E_0}\Big) + 2{k_0}}} \end{align}\] Thus, as \(\xi\to\infty\) \[U (\infty) = 0 =\mp{\sqrt{- 2\Big(\frac{V_0^2}{2} + {E_0}\Big) + 2{k_0}}} \quad\] implies 38 .

Substituting 37 38 into 33 , we obtain the explicit formula \[\label{euler-sol-theta} \theta^2 = 2{k_0} \;\phi(\xi)\tag{39}\] with \(\phi\) the function \[\label{defphi} \phi(\xi) := \xi \big( \sqrt{1+\xi^2} - \xi \big)\tag{40}\]

The sign of \(\phi\) provides a constraint on the existence of \(\theta\). From 39 , \(\theta(\xi)\) is well-defined iff \(k_0 = E_0 + \frac{V_0^2}{2}>0\). Under this restriction, we derive an explicit family of solutions depending on two parameters \(V_0\), \(E_0\) \[\label{eulersolution} \begin{align} &\theta = \pm \sqrt{2 {k_0} \phi(\xi)} \quad &&U = \mp \sqrt{\frac{k_0}{2\phi(\xi)}} \Bigg[ \frac{1-2\phi(\xi)}{\sqrt{1+\xi^2}} \Bigg] \qquad &&&V = V_{0} \\ &W = \pm \sqrt{\frac{k_0}{2\phi(\xi)}} \frac{\xi}{\sqrt{1+\xi^2}} \quad &&P = \bigg(k_0 - \frac{V_0^2}{2}\bigg) - k_0 \frac{\xi}{\sqrt{1+\xi^2}} \quad &&& \end{align}\tag{41}\] Observe that \(\theta(\xi)\) can be either positive or negative. If \(\theta\) is positive, then the flow is directed inward near the plane \(z = 0\) and upward near the vortex line. Conversely, if \(\theta\) is negative, the flow is directed outward near the plane \(z = 0\) and downward near the vortex line, see Figure 2.

Note that \(P(\xi=0) = E_0\) and \(P(\xi=\infty)= - \frac{V_0^2}{2}\) and, as \(P(\xi)\) is decreasing, the values of \(P\) range from \(-\frac{V_0^2}{2}\) to \(E_0\). The constraint \(k_0 >0\) may be interpreted as \[k_0 = E_0 + \frac{V_0^2}{2} = P(\xi=0) - P(\xi=\infty)>0 \, .\] As \(p = \frac{P}{r^2}\) and \(v = \frac{V}{r}\) the constraint may be thought as presenting a difference between pressure induced forces at the top of the vortex core and the intersection of the vortex core with the boundary. \(V_0\) is a measure of the circulation around the vortex core.

Returning to the original variables through 20 we obtain \[\begin{align} &u(r,z) = \mp \sqrt{\frac{k_0}{2}} \frac{ (\sqrt{r^2 + z^2} - z )^2 }{ r \sqrt{ z (r^2 + z^2) (\sqrt{r^2 + z^2} -z) }} \quad &&v(r,z) = \frac{V_0}{r} \\ &w(r,z) = \pm \sqrt{\frac{k_0}{2}} \frac{ z }{ \sqrt{ z (r^2 + z^2) (\sqrt{r^2 + z^2} -z) }} \quad &&p(r,z) = \frac{1}{r^2} \Big ( \big (k_0 - \frac{V_0^2}{2}\big) - k_0 \frac{z}{\sqrt{r^2 + z^2}} \Big ) \end{align}\] The streamlines of the vector field \((u,w)\) are sketched in Figure 2, showing the two type of flows. In Figure 3 contour plots of the pressure are shown for four different choices of the parameters \(E_0, \;V_0\).

4 Do the axisymmetric Euler equations admit self-similar discontinuous solutions?↩︎

We consider the Euler equations 27 - 30 and focus on solutions where \(\theta\) vanishes and \(V\) exhibits discontinuities. We study weak solutions for the system 27 - 30 and investigate the existence of discontinuous solutions subject to the boundary conditions used before.

4.1 Discontinuous solutions - Weak Formulation↩︎

Let \(\varphi \in C^{\infty}_c ((0,\infty))\). Multiplying the equations 27 - 30 by \(\varphi\) and integrating by parts, we obtain \[\tag{42} \begin{align} \tag{43} -\int_0^\infty \bigg(\frac{\theta^2}{2} + (1+\xi^2)P \bigg) \varphi'(\xi) d\xi &= -\int_0^\infty \xi V^2(\xi) \varphi(\xi) d\xi, \\ \tag{44} -\int_0^\infty \bigg(\theta(\xi) V(\xi)\bigg) \varphi'(\xi) d\xi &= -\int_0^\infty U(\xi) V(\xi) \varphi(\xi) d\xi, \;\\ \tag{45} -\int_0^\infty \bigg(\theta^2 - \xi \Big(\frac{\theta^2}{2}\Big)' + P \bigg) \varphi'(\xi) d\xi &= 0,\\ \tag{46} -\int_0^\infty \theta \varphi'(\xi) d\xi &= -\int_0^\infty U(\xi) \varphi(\xi) d\xi. \end{align}\]

In the definition 1, \(BV\) denotes the space of functions of bounded variation. This property ensures that right and left limits of \((\theta, V, P)\) exist at any point \(\xi\). Since \(\theta\) belongs to \(W^{1,1}((0,\infty))\), then \(\theta\) is absolutely continuous and of bounded variation. The derivative of \(\theta\), i.e. \(\theta'=-U\), exists almost everywhere and is Lebesgue integrable on \((0,\infty)\).

Suppose \((\theta,V,P)\) is a weak solution of 27 - 30 . Using the theory of Sobolev spaces, [18], [19], one may explore properties of system 42 . Recalling that \(\theta\) and as a consequence \(\theta^2\) are absolutely continuous, we get that \(\theta^2 + \xi \Big(\frac{\theta^2}{2}\Big)' \in L^1_{loc}((0,\infty))\). Also, \(P \in L^{\infty}((0,\infty)) \subset L^1_{loc}((0,\infty))\). Then, using [18], 45 implies there exists a constant \(A\) such that \[\label{f1} \theta^2 - \xi \Big(\frac{\theta^2}{2}\Big)' + P = A,\tag{47}\] almost everywhere (a.e). The same lemma applied to 43 , 44 and 46 as follows. Using [18] we rewrite 43 , 44 and 46 in the form \[\begin{align} &\quad -\int_0^\infty \bigg(\frac{\theta^2}{2} + (1+\xi^2)P + \int_0^\xi \zeta V^2(\zeta) d\zeta\bigg) \varphi'(\xi) d\xi = 0. \\ &\quad -\int_0^\infty \bigg(\theta(\xi) V(\xi) + \int_0^\xi U(\zeta) V(\zeta) d\zeta \bigg) \varphi'(\xi) d\xi = 0. \\ &\quad-\int_0^\infty \bigg(\theta + \int_0^\xi U(\zeta) d\zeta\bigg) \varphi'(\xi) d\xi = 0 \quad \forall \, \varphi. \end{align}\] The integrals \(\int_0^\xi \zeta V^2(\zeta) d\zeta\), \(\int_0^\xi U(\zeta) V(\zeta) d\zeta\) and \(\int_0^\xi U(\zeta) d\zeta\) for \(U=-\theta' \in L^1\) are well-defined. Using [18], there are constants \(B\), \(C\) and \(D\) such that \[\begin{align} \tag{48} \frac{\theta^2}{2} + (1+\xi^2)P + \int_0^\xi \zeta V^2(\zeta) d\zeta = B \quad\textrm{a.e.}, \\ \tag{49} \theta(\xi) V(\xi) + \int_0^\xi U(\zeta) V(\zeta) d\zeta = C \quad\textrm{a.e.}, \\ \tag{50} \theta(\xi) + \int_0^\xi U(\zeta) d\zeta = D \quad\textrm{a.e.}. \end{align}\] Letting \(\xi\to a_-\) and \(\xi\to b_+\) for \(0<a,b<\infty\) in 47 and 48 50 , one may conclude that \((\theta,V,P)\) a weak solution of 27 - 29 satisfies for \(a<b\) \[\label{weak2} \begin{align} \bigg(\frac{\theta^2(\xi)}{2} + (1+\xi^2)P(\xi) \bigg)\Bigg|_{a-}^{b_+} &= - \int_a^b \zeta V^2(\zeta) d\zeta, \\ \bigg(\theta(\xi) V(\xi)\bigg) \Bigg|_{a-}^{b_+} &= - \int_a^b U(\xi) V(\xi) d\xi, \;\\ \bigg(\theta^2 - \xi \Big(\frac{\theta^2}{2}\Big)' + P(\xi) \bigg)\Bigg|_{a-}^{b_+} &= 0, \\ \theta(\xi) \Big|_{a-}^{b_+} &= - \int_a^b U(\xi) d\xi, \end{align}\tag{51}\]

Jump Discontinuities↩︎

Let \((\theta,V,P)\) be a weak solution of 27 - 30 defined by definition 1. Utilizing 51 , we will compute the Rankine-Hugoniot jump conditions associated with the system. Suppose there exists a discontinuity at some point \(\xi = \sigma\), \(0<\sigma<\infty\) as in Figure 4. Letting \(a,b \to \sigma\), 51 reduces to \[\tag{52} \begin{align} \tag{53} \frac{1}{2}\Big(\theta_+^2 - \theta_-^2\Big) + (1+\sigma^2)\Big(P_+ - P_-\Big) &= 0, \\ \tag{54} \theta_+ V_+ - \theta_- V_- &= 0, \\ \tag{55} \bigg(\theta_+^2 - \sigma \Big(\frac{\theta_+^2}{2}\Big)' - \Big(\theta_-^2 - \sigma \Big(\frac{\theta_-^2}{2}\Big)'\Big) \bigg) + \bigg(P_+ - P_- \bigg) &= 0, \\ \tag{56} \theta_+ - \theta_- &= 0, \end{align}\] where \(\theta_\pm = \theta(\sigma\pm)\), \(V_\pm = V(\sigma\pm)\) and \(P_\pm = P(\sigma\pm)\) denote the one-sided limits at \(\xi =\sigma\). Note that all one-sided limits exist and are finite since \((\theta,V,P)\) are of bounded variation.

The last equation in 52 implies that \(\theta\) must be continuous for any \(\xi \in (0, \infty)\), that is to say \[\theta_- =\theta_+,\] and thus, 53 55 provide the Rankine-Hugoniot jump conditions \[\tag{57} \begin{align} \tag{58} \llbracket P \rrbracket = 0, \\ \llbracket \theta V \rrbracket = 0, \\ \tag{59} \llbracket \theta \, \theta' \rrbracket = 0, \end{align}\] where \(\llbracket F \rrbracket = F_+ - F_-\) is the jump operator. It is easy to see that \(P(\xi)\) must also be continuous for any \(\xi \in (0, \infty)\). On the other hand, for \({\theta}(\sigma)=0\) 57 implies that \(V(\xi)\) and \(\theta'(\xi)\) have a jump discontinuity at \(\xi = \sigma\). (The case where \(\theta(\sigma) \neq 0\) implies that \(V\) is continuous and leads to solutions described in Section 3.)

4.1.1 Non-Existence of Discontinuous Solutions↩︎

Consider \((\theta,V,P)\) a weak solution of 27 - 29 given by definition [weak-def] with a discontinuity at a fixed point \(\xi = \sigma\), \(\sigma \in (0,\infty)\). Recall that \(\theta'\) and \(V\) may exhibit discontinuities, it suggests looking for solutions in the form \[\begin{align*} \theta(\xi) = \left\{ \begin{align} & \theta_{-}(\xi), \\ &\theta_{+}(\xi) \end{align} \right. \quad \text{and} \quad V(\xi) = \left\{ \begin{align} & V_-(\xi) \,\, , \quad &&\xi \in (0,\sigma) \\ & V_+ (\xi)\,\, , \quad &&\xi \in (\sigma,\infty) \end{align} \right. \end{align*}\] We solve the system 27 - 29 on \((0,\sigma)\) and \((\sigma,\infty)\) independently. Similar calculations to those described in Section 3 lead again to a solution \(\theta\) in form 34 on each domain. Note that \(V(\xi)\) is again a constant function with \(V_-\) and \(V_+\) constants. Since discontinuous solutions occur when \({\theta}(\sigma)=0\), the relation for \(\theta\) becomes \[\label{th-disc1} \frac{\theta_\pm^2}{2} = {k_\pm} \bigg(\phi(\xi) - \phi(\sigma)\bigg) + \Big({k_\pm}-\frac{V_\pm^2}{2} - 2{A_\pm} + {E_\pm}\Big) \big(\xi^2 - \sigma^2\big).\tag{60}\] where \(k_\pm, A_\pm\) and \(E_\pm\) are integration constants. We next impose the boundary conditions 36 , that is \[{\theta}_-(0) = 0 \quad \textrm{and} \quad {\theta'_+}(\xi) \to 0 \quad \textrm{as} \quad \xi \to \infty,\] which imply \[\begin{align} \tag{61} {k_-}-\frac{V_-^2}{2} - 2{A_-} + {E_-} &= - {k_-} \frac{1}{\sigma^2} \phi(\sigma), \\ \tag{62} k_+ - \frac{{V_+}^2}{2} - 2{A_+} + {E_+} &= 0. \end{align}\] The derivation of 61 is direct. To derive relation 62 , observe that 60 yields \[\begin{align} \theta'_+ &= \frac{1}{\theta} \bigg[k_+ \phi'(\xi) + 2\bigg(k_+ - \frac{V_+^2}{2}- 2A_+ + E_+\bigg) \xi \bigg] \\ &= \frac{k_+ \frac{\phi'(\xi)}{\xi} + 2\bigg(k_+ - \frac{V_+^2}{2}- 2A_+ + E_+\bigg)}{\pm\sqrt{2k_+ \frac{\phi(\xi)-\phi(\sigma)}{\xi^2} + 2\bigg(k_+ - \frac{V_+^2}{2}- 2A_+ + E_+\bigg) \frac{\xi^2-\sigma^2}{\xi^2}}}\\ &\xrightarrow []{\xi\to\infty} \pm \sqrt{2\bigg(k_+ - \frac{V_+^2}{2}- 2A_+ + E_+\bigg)} , \end{align}\] Therefore, as \(\xi \to \infty\) \[\theta'(\infty) = 0 = k_+ - \frac{{V_+}^2}{2} - 2{A_+} + {E_+},\] provides 62 . Substituting now 61 62 into 60 , we obtain the explicit formula \[\begin{align} \label{th-disc2} \frac{\theta^2}{2}(\xi) = \left\{ \begin{aligned} &{k_-} \Bigg[\phi(\xi) - \phi(\sigma) - \frac{\xi^2-\sigma^2}{\sigma^2}\phi(\sigma)\Bigg] \,\, , \quad &&\xi \in (0,\sigma) \\ &{k_+} \bigg[ \phi(\xi) - \phi(\sigma) \bigg] \,\, , \quad &&\xi \in (\sigma,\infty) \end{aligned} \right. \end{align}\tag{63}\] From 33 and 63 , we then derive an explicit family of discontinuous solutions depending on parameters \(V_\pm, E_\pm\) and \(\sigma\) \[\label{euler-disc} \begin{align} &U = \left\{ \begin{aligned} &\frac{k_-}{\theta_-} \,\, \Bigg[\frac{2\phi(\sigma)}{\sigma^2} \xi - \phi'(\xi)\Bigg] \\ & - \frac{k_+}{\theta_+} \,\, \phi'(\xi) \end{aligned} \right.\, \, && \qquad W = \left\{ \begin{align} & \frac{k_-}{\theta_-} \;\frac{\xi}{\sqrt{1+\xi^2}} \,\, , \quad &&\xi \in (0,\sigma) \\ &\frac{k_+}{\theta_+} \Bigg[\frac{\xi}{\sqrt{1+\xi^2}} - 2 \phi(\sigma)\Bigg] \,\, , \quad &&\xi \in (\sigma,\infty) \end{align} \right.\\ &V = \left\{ \begin{align} & V_- \\ & V_+ \end{align} \right.\qquad && \qquad P = \left\{ \begin{aligned} & E_- - k_- \frac{\xi}{\sqrt{1+\xi^2}} \,\, , \quad &&\xi \in (0,\sigma) \\ & E_+ - {k_+} \Bigg(\frac{\xi}{\sqrt{1+\xi^2}} - 2\phi(\sigma) \Bigg) \,\, , \quad &&\xi \in (\sigma,\infty) \end{aligned} \right. \end{align}\tag{64}\] where \(E_- = k_-\big(1 + \frac{\phi(\sigma)}{\sigma^2}\big) - \frac{V_-^2}{2}\) and \(E_+ = k_+ ( 1-2\phi(\sigma)) -\frac{V_+^2}{2}\). Although solutions \(U(\xi)\) and \(W(\xi)\) tend to infinity near the discontinuity at \(\xi=\sigma\), they achieve the required regularity in definition [weak-def]. In particular, examining the denominator yields \[\begin{align*} U(\xi) \sim W(\xi) \sim \left\{ \begin{align} &\frac{k_-}{\theta_-} \sim \sqrt{\frac{k_-}{\Big(\phi(\xi) - \phi(\sigma)\Big) - \frac{\phi(\sigma)}{\sigma^2} \Big(\xi^2 - \sigma^2\Big)}}\sim \frac{\sqrt{k_-}}{\sqrt{\lvert \xi -\sigma \rvert}} \,\, , \quad &&\textrm{as} \quad \xi \to \sigma_- \\ & \frac{k_+}{\theta_+} \sim \sqrt{\frac{k_+}{\Big(\phi(\xi) - \phi(\sigma)\Big)}} \sim \frac{\sqrt{k_+}}{\sqrt{\lvert \xi -\sigma \rvert}} \,\, , \quad &&\textrm{as} \quad \xi \to \sigma_+ \end{align} \right. \end{align*}\] and therefore the singularities \(\xi=\sigma\) are integrable, i.e. \(U,W \in L^1_{loc}((0,\infty))\). Analysis of these singularities leads to the following theorem.

4.1.2 Non-Existence of Multiple Discontinuities↩︎

Let \((\theta,V,P)\) be a weak solution of 27 - 30 with discontinuities at multiple points \(\xi = \sigma_i\), \(0<\sigma_i <\infty\), \(i>0\). We consider first the case where the solution exhibits discontinuities at two points \(\sigma_1\) and \(\sigma_2\), \(\sigma_1<\sigma_2\). This suggests seeking solutions in the form \[\begin{align*} \theta (\xi) = \left\{ \begin{align} & \theta_{1}(\xi) \\ &\theta_{2}(\xi) \\ &\theta_{3}(\xi) \end{align} \right. \quad \text{and} \quad V (\xi)= \left\{ \begin{align} & V_1 (\xi)\,\, , \quad &&\xi \in (0,\sigma_1) \\ & V_2 (\xi)\,\, , \quad &&\xi \in (\sigma_1,\sigma_2) \\ & V_3 (\xi) \,\, , \quad &&\xi \in (\sigma_2,\infty) \end{align} \right. \end{align*}\] To examine whether such solutions exist, it is sufficient to derive the solution \(\theta\) by solving system 27 - 29 on each domain independently. Same calculations as those described in Section 3 lead to functions \(\theta_1, \theta_2\) and \(\theta_3\) of the form 34 , that is \[\label{eq6} \frac{\theta_j^2}{2} = {k_j} \phi(\xi) + \Big({k_j}-\frac{V_j^2}{2} - 2{A_j} + {E_j}\Big) \xi^2 + \Big({E_j} - {A_j}\Big), \quad j=1,2,3\tag{67}\] where \(A_j,E_j, k_j\) are integration constants. Note that \(V(\xi)\) is again a constant function with \(V_1, V_2\) and \(V_3\) constants. Recall that discontinuous solutions occur when \(\theta(\sigma_1)=\theta(\sigma_2)=0\), it implies \[\theta_1(\sigma_1) = \theta_2(\sigma_1) = \theta_2(\sigma_2) = \theta_3(\sigma_2) = 0,\] and thus, \(\theta\) takes the following form \[\label{eq7} \frac{\theta^2}{2} = \left\{ \begin{align} &\frac{\theta_{1}^2}{2} = {k_1} \Bigg[\phi(\xi) - \phi(\sigma_1) + \Big({k_1}-\frac{V_1^2}{2} - 2{A_1} + {E_1}\Big) (\xi^2-\sigma_1^2)\Bigg] \,\, , \quad &&\xi \in (0,\sigma_1) \\ &\frac{\theta_{2}^2}{2} = {k_2} \Bigg[\phi(\xi) - \phi(\sigma_1) - \frac{\phi(\sigma_2)-\phi(\sigma_1)}{\sigma_2^2 - \sigma_1^2} (\xi^2-\sigma_1^2) \Bigg] \,\, , \quad &&\xi \in (\sigma_1,\sigma_2) \\ &\frac{\theta_{3}^2}{2}= {k_3} \bigg[ \phi(\xi) - \phi(\sigma_2) + \Big({k_3}-\frac{V_3^2}{2} - 2{A_3} + {E_3}\Big) (\xi^2-\sigma_2^2) \bigg] \,\, , \quad &&\xi \in (\sigma_2,\infty) \end{align} \right.\tag{68}\]

Imposing now the boundary conditions 36 into 68 yields the explicit formula of discontinuous solution \(\theta\) \[\label{disc-sol3} \frac{\theta^2}{2} = \left\{ \begin{align} &{k_1} \Bigg[\phi(\xi) - \phi(\sigma_1) - \frac{\xi^2-\sigma_1^2}{\sigma_1^2}\phi(\sigma_1)\Bigg] \,\, , \quad &&\xi \in (0,\sigma_1) \\ &{k_2} \Bigg[\phi(\xi) - \phi(\sigma_1) - \frac{\phi(\sigma_2)-\phi(\sigma_1)}{\sigma_2^2 - \sigma_1^2} (\xi^2-\sigma_1^2) \Bigg] \,\, , \quad &&\xi \in (\sigma_1,\sigma_2) \\ &{k_3} \bigg[ \phi(\xi) - \phi(\sigma_2) \bigg] \,\, , \quad &&\xi \in (\sigma_2,\infty) \end{align} \right.\tag{69}\] Note that functions \(\theta_1\) and \(\theta_3\) have the same form as \(\theta_-\) and \(\theta_+\) in 63 for \(\sigma=\sigma_1\) and \(\sigma=\sigma_2\) respectively. Likewise \(k_-\) and \(k_+\), we conclude that \(k_1>0\) and \(k_3>0\) using the same reasoning. To define the sign of \(k_2\) in the remaining case, we consider the following lemma.

Lemma [lemma1] implies that \(k_2>0\). Considering the signs of constants \(k_1,k_2\) and \(k_3\), we proceed to the main theorem about the non-existence of solutions with discontinuities at two points.

Recall theorem 1 and theorem 2, we observe that jump condition 59 is not satisfied neither near the boundary nor across the vortex line, despite the existence of an intermediate region. Therefore, a similar outcome will also be obtained if a finite number of intermediate layers are considered. This leads to the following corollary.

4.2 Interaction of vortex with boundary↩︎

Motivated by the study of Euler equations presented in the previous sections, we are interested in extending it to a class of flows where there is mass input or loss through the vortex line. In other words, we assume that the vortex line can be either a source or a sink of the fluid motion. Since that the vortex line resembles the tornado core, this assumption is not unnatural for a tornado-like flow. In terms of the self-similar functions 20 , the assumption leads to \[U(\xi) \to U_\infty \quad \textrm{as} \quad \xi \to \infty,\] or in terms of \(\theta\) \[\label{bc-new} \theta'(\xi) \to -U_\infty \quad \textrm{as} \quad \xi \to \infty,\tag{74}\] where \(U_\infty\) is a non-zero constant.

4.2.1 Continuous solutions↩︎

Consider first the case where solutions are continuous. Same calculations as those described in Section 3 lead again to solutions in form 33 . The boundary condition 74 on the vortex line applied to 34 yields \[\label{bc-new1} \frac{1}{2} \, U_\infty^2 = k_0 - \frac{V_0^2}{2}- 2A_0 + E_0.\tag{75}\] To derive this, note that 34 implies \[\begin{align} \theta' = \frac{1}{\theta} \bigg[k_0 \phi'(\xi) + 2\bigg(k_0 - \frac{V_0^2}{2}- 2A_0 + E_0\bigg) \xi\bigg] \xrightarrow []{\xi\to\infty} {\pm\sqrt{2\bigg(k_0 - \frac{V_0^2}{2}- 2A_0 + E_0\bigg)}}. \end{align}\] Thus, as \(\xi\to \infty\) \[\theta'(\infty) = -U_\infty = {\pm\sqrt{2\bigg(k_0 - \frac{V_0^2}{2}- 2A_0 + E_0\bigg)}}.\] In addition, we impose a no-penetration boundary condition at \(\xi=0\), i.e. \(\theta(0)=0\). Applying this to 34 provides \[\begin{align} \label{bc-new2} A_0 = E_0. \end{align}\tag{76}\] As a result the relation 75 reduces to \[\label{bc-new3} \frac{1}{2} U_\infty^2 = k_0 - \frac{V_0^2}{2}- E_0.\tag{77}\] Substituting 76 and 77 into 33 , we obtain an explicit family solutions of \(\eqref{theta-eq} - \eqref{p-eq}\) depending on parameters \(U_\infty, V_0\) and \(E_0\) \[\begin{align} &\frac{\theta^2}{2} = k_0 \phi(\xi) + \frac{1}{2} U_\infty^2 \xi^2, \qquad U = - \frac{1}{\theta} \bigg[\frac{k_0}{\sqrt{1+\xi^2}} \;\big(1-2\phi(\xi)\big) + U_\infty^2 \xi \bigg], \\ &V = V_{0}, \qquad W = \frac{k_0}{\theta} \frac{\xi}{\sqrt{1+\xi^2}}, \qquad P = E_0 - k_0 \frac{\xi}{\sqrt{1+\xi^2}}, \end{align}\] where \[k_0= \frac{1}{2} U_\infty^2 + \frac{V_0^2}{2}+E_0.\] We next investigate the existence of discontinuous solutions in this framework.

4.2.2 Nonexistence of Discontinuous Solutions↩︎

Let \((\theta,V,P)\) be a weak solution of 27 - 30 with a discontinuity at a fixed point \(\xi = \sigma\), \(\sigma \in (0,\infty)\), which satisfies the boundary condition 74 . This suggests seeking solutions in the form \[\begin{align*} \theta (\xi) = \left\{ \begin{align} & \theta_{-} (\xi), \\ &\theta_{+} (\xi) \end{align} \right. \quad \text{and} \quad V (\xi) = \left\{ \begin{align} & V_- (\xi) \,\, , \quad &&\xi \in (0,\sigma) \\ & V_+ (\xi)\,\, , \quad &&\xi \in (\sigma,\infty) \end{align} \right. \end{align*}\] Repeating calculations described in Section 3, we obtain again solutions in form 34 on \((0,\sigma)\) and \((\sigma,\infty)\) respectively. The system 27 - 29 is solved on each domain independently. Note that \(V(\xi)\) is again a constant function with \(V_-\) and \(V_+\) constants. Recall that discontinuous solutions occur when \({\theta}(\sigma)=0\), the relation for \(\theta\) becomes \[\label{th-disc3} \frac{\theta_\pm^2}{2} = {k_\pm} \bigg(\phi(\xi) - \phi(\sigma)\bigg) + \Big({k_\pm}-\frac{V_\pm^2}{2} - 2{A_\pm} + {E_\pm}\Big) \big(\xi^2 - \sigma^2\big).\tag{78}\] where \(k_\pm, A_\pm\) and \(E_\pm\) are integration constants.

We impose a no-penetration boundary condition at the boundary and the condition 74 at the vortex line, namely \[\theta_-(0) = 0, \quad \textrm{and} \quad \theta'_+(\xi) \to -U_\infty \quad \textrm{as} \quad \xi \to \infty,\] which imply \[\begin{align} \tag{79} {k_-}-\frac{V_-^2}{2} - 2{A_-} + {E_-} &= - {k_-} \frac{1}{\sigma^2} \\ \tag{80} k_+ - \frac{V_+^2}{2}- 2A_+ + E_+ &= \frac{1}{2} U_\infty^2 \end{align}\] To derive 80 , note that 78 provides \[\theta'_+ \xrightarrow []{\xi\to\infty} {\pm\sqrt{2\bigg(k_+ - \frac{V_+^2}{2}- 2A_+ + E_+\bigg)}},\] and thus, \[\theta'_+(\infty) = - U_\infty = \pm\sqrt{2\bigg(k_+ - \frac{V_+^2}{2}- 2A_+ + E_+\bigg)}.\] Substituting 79 and 80 into 78 , we obtain the explicit formula \[\begin{align} \label{disc-sol2} \frac{\theta^2}{2} = \left\{ \begin{aligned} &{k_-} \Bigg[\phi(\xi) - \phi(\sigma) - \frac{\xi^2-\sigma^2}{\sigma^2}\phi(\sigma)\Bigg] \,\, , \quad &&\xi \in (0,\sigma) \\ &{k_+} \bigg[ \phi(\xi) - \phi(\sigma) \bigg] + \frac{1}{2} U_\infty^2 (\xi^2-\sigma^2) \,\, , \quad &&\xi \in (\sigma,\infty) \end{aligned} \right. \end{align}\tag{81}\] where \(k_+, k_-\) are constants. From 33 and 81 , one may derive the explicit formulas of discontinuous solutions \(U, W\) and \(P\).

Examining now whether discontinuous solutions in the form 81 are feasible leads to the following theorem.

5 Stationary self-similar axisymmetric Navier-Stokes equations↩︎

In this section we study self-similar solutions for the stationary axisymmetric Navier-Stokes equations 1 . The ansatz 20 leads to the system 25 in the variables \((\theta, V, P)\), where the velocities \((U, W)\) are determined from the stream function \(\theta\) via 24 . We first provide a convenient reformulation of 25 as a coupled integrodifferential system and then study the limiting behavior as the viscosity \(\nu \to 0\), the form of the boundary layer, and identify conditions on admissible solutions of the Euler equations.

5.1 Formulation via an integro-differential system↩︎

The system 25 for \(\nu > 0\) gives after an integration \[\begin{align} \tag{83} \frac{\theta^2}{2} + (1 + \xi^2) P &= \nu \Big [ \xi \theta - (1 + \xi^2) \theta^\prime \Big ] - \int_{0}^{\xi} s V^2 \,ds + A_0 \\ \tag{84} \theta^2 - \xi \Big(\frac{\theta^2}{2}\Big)' + P &= \nu \Big[\xi\theta - \xi^2 \theta' - \xi (1+\xi^2) \theta'' \Big] + {E_0} , \end{align}\] where \(A_0\), \(E_0\) are integration constants.

We impose a no-slip boundary condition \(\vec{u} = 0\) at the boundary \(z=0\), which implies \(u = v = w = 0\) at \(z = 0\). Expressed in terms of self-similar functions it yields \[\label{bd1} U(0) = V(0) = W(0) = 0, \quad \textrm{and thus} \quad \theta(0) = \theta'(0) = 0.\tag{85}\] We also require that no mass is added or lost through the vortex line, that is \(u (2 \pi r) \to 0\) as \(r \to 0\). In the self-similar setting the condition at the vortex core becomes \[\label{bd2} U = - {\theta}' \to 0 \quadas \; \xi \to \infty\tag{86}\] Since 26 is second order, an additional condition is required for \(V(\xi)\). We request the vortex line to have constant swirl: \[\label{bd3} V \to V_\infty, \,\, \textrm{as} \quad \xi \to \infty.\tag{87}\]

Then 85 together with 83 and 84 imply that \(E_0 = A_0 = P(0)\). Eliminating \(P\) from 83 and 84 gives \[\begin{align} \xi (1+\xi^2) \Big(\frac{\theta^2}{2}\Big)' - (1+2\xi^2) \frac{\theta^2}{2}& + \nu \bigg[\xi^3 \theta + (1-\xi^4) \theta' - \xi(1+\xi^2)^2 \theta'' \bigg] \nonumber\\ &= - \int_{0}^{\xi} s V^2 \,ds - {E_0} \xi^2 \end{align}\] Upon dividing by \(\xi^2 (1+\xi^2)^{\frac{3}{2}}\) and using the calculus identities \[\begin{align} \tag{88} &\frac{d}{d\xi} \frac{1}{\xi (1+\xi^2)^{\frac{1}{2}} } = - \frac{1+2\xi^2}{\xi^2 (1+\xi^2)^{\frac{3}{2}}} \qquad \qquad &&\frac{d}{d\xi} \frac{ (1+\xi^2)^{\frac{1}{2}} }{\xi} = - \frac{1}{ \xi^2 (1+\xi^2)^{\frac{1}{2}}} \\ \tag{89} &\frac{d}{d\xi} \frac{1}{ (1+\xi^2)^{\frac{1}{2}} } = - \frac{\xi}{ (1+\xi^2)^{\frac{3}{2}}} \qquad \qquad &&\frac{d}{d\xi} \Big ( \frac{1}{\xi \sqrt{1+ \xi^2} \big (\xi + \sqrt{1+\xi^2} \big )^2 } \Big ) = - \frac{1}{ \xi^2 (1 + \xi^2)^{\frac{3}{2}}} \end{align}\] we obtain after some rearrangements the equation \[\label{diffeq10} \frac{d}{d\xi} \left ( \frac{1}{\xi\sqrt{1+\xi^2}} \frac{\theta^2}{2} - \nu \Big[\frac{1}{\sqrt{1+\xi^2}} \theta + \frac{\sqrt{1+\xi^2}}{\xi} \theta' \Big] \right ) = - \frac{1}{\xi^2 (1+\xi^2)^{\frac{3}{2}}} \int_{0}^{\xi} s V^2 \,ds - \frac{d}{d\xi} \Big ( \frac{E_0 \xi }{ (1+\xi^2)^{\frac{1}{2}}} \Big ) \, .\tag{90}\]

Using the boundary condition 86 , \(\theta^\prime (\infty) = 0\), which also implies that \(\lim_{\xi \to \infty} \frac{\theta(\xi)}{\xi} = 0\), we integrate 90 over the interval \((\xi, \infty)\) to obtain \[\frac{\theta^2}{2} - \nu \Big ( 1+\xi^2) \theta' + \xi \theta \Big ) = \xi(1+\xi^2)^{\frac{1}{2}} \int_\xi^\infty \frac{1}{\zeta^2(1+\zeta^2)^\frac{3}{2}} \bigg( \int_0^\zeta s V^2(s) ds\bigg) d\zeta + \,\, {E_0} \Big ( \xi (1+\xi^2)^{\frac{1}{2}} - \xi^2 \Big ) \, .\] The latter is an integrodifferential equation depending on \(\mathcal{G} (V ; \xi)\) a functional on \(V\), \[\label{defG} \begin{align} \mathcal{G}\left(V ; \xi\right) &:= \xi \sqrt{1+\xi^2} \int_\xi^\infty \frac{1}{\zeta^2(1+\zeta^2)^\frac{3}{2}} \bigg( \int_0^\zeta s V^2(s) ds\bigg) d\zeta \\ &= \xi \sqrt{1+\xi^2} \int_\xi^\infty \frac{1}{s \sqrt{1+s^2} \big (s + \sqrt{1+s^2} \big )^2 } s V^2 ds + \frac{1}{ (\xi + \sqrt{1 + \xi^2} )^2} \int_0^\xi s V^2 ds \, , \end{align}\tag{91}\] where the last equality follows via integration by parts using 89 .

In summary, the system 25 with boundary conditions 85 87 is transformed to a coupled integrodifferential system for \((\theta, V)\), \[\tag{92} \begin{align} \tag{93} \frac{\theta^2}{2} - \nu \bigg[(1+\xi^2) \theta' + \xi \theta \bigg] &= \mathcal{G}\left( V ; \xi\right) + E_0 \phi (\xi) \\ \tag{94} \nu (1+\xi^2) V'' + \Big(3\nu\xi - \theta\Big) V' &= 0 \, , \end{align}\] where \(\mathcal{G} (V ; \xi)\) is the functional in 91 and \(\phi(\xi)\) the function in 40 that played a prominent role in the Euler equations. The system is supplemented with the boundary conditions \[\label{bc95xi} \begin{align} \theta (0) = V (0) = 0 , \,\, &\textrm{ at } \quad \xi = 0 \\ V \to V_\infty , \,\, &\textrm{ as} \quad \xi \to \infty \end{align}\tag{95}\] and satisfies \(\theta^\prime(0) = \theta^\prime(\infty) = 0\). The parameter \(E_0= P(0)\) reflects on the pressure \(p = \frac{1}{r^2} P(\xi)\) and \(P\) may be computed by 84 .

5.2 Alternative equivalent formulations↩︎

Next, we derive two more equivalent formulations of the system 92 with simpler structure. This will benefit the analysis and provides a different perspective on the problem.

First we set \[\label{firsttrans} \hat{\theta} (\xi) = \sqrt{1+\xi^2} \, \theta(\xi) \quad \textrm{and} \quad \hat{V} (\xi) =V(\xi), \quad \xi\in[0,\infty),\tag{96}\] and obtain \[\label{eq-ns-v2} \begin{align} \frac{\hat{\theta}^2}{2} - \nu (1+\xi^2)^\frac{3}{2} \, \hat{\theta}' &= (1+\xi^2) \Big ( \mathcal{G} \big (\hat{V}(\xi) ; \xi \big) + E_0 \phi (\xi) \Big ) , \\ \nu\Big ( (1+\xi^2)^\frac{3}{2} \hat{V}'\Big )^\prime &= \hat{\theta} \hat{V}', \end{align}\tag{97}\] with boundary conditions 95 .

An alternative form of 92 is derived by introducing a change of variable \(\xi \to x\) such that \(\frac{d}{dx} = (1 + \xi^2)^{\frac{3}{2}} \frac{d}{d\xi}\). This is achieved by defining \(x\) by \[\label{def95x} x = \frac{\xi}{\sqrt{1+ \xi^2}} \quad\textrm{with inverse thansformation} \quad \xi = \frac{x}{\sqrt{1- x^2}} \, .\tag{98}\] The change of variables \(\xi \to x\) is a surjective map \(T: [0, \infty) \to [0, 1)\) and has an interpretation as a change from cylindrical \((r,\vartheta,z)\) to spherical \((\rho,\vartheta,\varphi)\) coordinates, see Figure 5. Since \(\xi = \frac{z}{r}\), it follows \(\xi = \frac{\cos \varphi}{\sin \varphi} = \cot \varphi\) and setting \(x = \cos\varphi\) yields \[\xi = \frac{\cos \varphi}{\sin \varphi} = \frac{x}{\sqrt{1- x^2}}.\] The functions \(\hat{\theta}(\xi)\) and \(\hat{V}(\xi)\) are expressed in terms of the variable \(x\) as follows \[\label{defvar} \Theta (x) = \hat{\theta}(\xi) \quad \textrm{and} \quad V (x) = \hat{V}(\xi),\tag{99}\] where \(\xi, \, x\) are related via 98 .

Using the transformation 98 , 99 , the system 92 reduces to the equivalent form \[\tag{100} \begin{align} \tag{101} \nu \, \frac{d \Theta}{dx} &= \frac{\Theta^2 (x)}{2} - \mathcal{F}\left(V ; x\right), \\ \tag{102} \nu \frac{d^2 {V}}{d^2x} &= \Theta(x) \, \frac{d {V}}{dx} , \\ \tag{103} &\qquad \Theta(0) = {V}(0) = 0 , \\ \tag{104} &\qquad {V} \to V_\infty, \,\,\textrm{ as} \,\;x \to 1, \end{align}\] where \(\mathcal{F}\) is \[\begin{align} \mathcal{F}\left(V ; x\right) &:= (1+\xi^2) \Big ( \mathcal{G} \big (\hat{V}(\xi) ; \xi \big) + E_0 \phi (\xi) \Big ) \Bigg |_{\xi = \frac{x}{\sqrt{1-x^2}}} \nonumber \\ &= \frac{x}{(1-x^2)^2}\int_x^1 \frac{1-t^2}{t^2} \bigg(\displaystyle \int_0^\frac{t}{\sqrt{1-t^2}} s \hat{V}^2(s) ds \bigg) \, dt + \, {E_0} \,\, \frac{x - x^2}{(1-x^2)^2} \label{defF} \\ &= \frac{x}{(1-x^2)^2}\Bigg[\int_x^1 \frac{1-t^2}{t^2} \bigg(\displaystyle \int_0^t \frac{\sigma}{(1-\sigma^2)^2} V^2(\sigma) d\sigma \bigg) \, dt + \, {E_0} \,\, (1 - x)\Bigg] \nonumber \\ &= \frac{1}{(1-x^2)^2}\Bigg[ x \int_x^1 \frac{1}{ (t+1)^2} V^2 (t) dt + (1-x)^2 \int_0^x \frac{t}{(1 - t^2)^2} V^2 (t) dt + \, {E_0} \, x (1 - x)\Bigg] \nonumber \end{align}\tag{105}\] The first expression in 105 is the definition of \(\mathcal{F}\left(V ; x\right)\) and the remaining equalities are obtained using 40 , 91 , the change of variables \(\xi = \frac{x}{\sqrt{1-x^2}}\) and the formula \(V(x) = \hat{V} \Big ( \frac{x}{ \sqrt{1-x^2}} \Big )\) connecting \(V(x)\) defined on \([0,1)\) with \(\hat{V}(\xi)\) defined on \([0,\infty)\).

Variants of 100 appear in [5], [6], [8], [13], [20] originating from different authors who initiate their studies by considering either spherical or cylindrical coordinates. Serrin [6] provided a detailed analysis of existence of solutions for 100 and a (indicative) bifurcation diagram for solutions. We refer to Section 7 for numerical computations leading to a computational bifurcation diagram, see Figure 11. As \(\Theta\) satisfies a Ricatti-type equation, there are regions of \(\nu\) for which \(\Theta (x)\) blows up for \(x \in [0,1)\) leading to non-existence of solutions for such \(\nu\). This is explained in the following section.

5.3 Properties and a-priori estimates↩︎

In the sequel we are interested in the limit \(\nu \to 0\) and the convergence of 100 to solutions of the Euler equations. To fix ideas we restrict to the case \(V_\infty > 0\). Then 102 , 104 imply that \(V(x)\) is strictly increasing and yields the representation formula \[\label{V-sol} V(x) = \displaystyle {V_\infty \, \frac{ \int_0^x e^{ \frac{1}{\nu} \int_0^t \Theta(s) ds} dt}{ \int_0^1 e^{ \frac{1}{\nu} \int_0^t \Theta(s) ds} dt} }\tag{106}\] leading to the following lemma.

We turn next to the equation 101 , \[\label{eqtheta} \nu \, \frac{d \Theta}{dx} = \frac{\Theta^2 (x)}{2} - \frac{x}{(1-x^2)^2}\Big(H(x) + \, {E_0} \,\, (1 - x)\Big)\tag{107}\] and proceed to establish properties for \(\Theta\). The functional \(H(x)\) in 105 is given by \[\label{H} \begin{align} H(x) = \int_x^1 \frac{1-t^2}{t^2} \bigg(\displaystyle \int_0^\frac{t}{\sqrt{1-t^2}} s \hat{V}^2(s) ds \bigg) \, dt \, , \end{align}\tag{108}\] here expressed in terms of the velocity \(\hat{V} (\xi)\) defined on \((0,\infty)\). The derivatives of \(H(x)\) are \[\begin{align} \frac{dH}{dx} (x) &=-\frac{1-x^2}{x^2} \int_0^\frac{x}{\sqrt{1-x^2}} s {\hat{V}}^2(s) ds < 0, \quad x \in [0,1), \tag{109} \\ \frac{d^2 H}{dx^2} (x) &= -\frac{2}{x^3} \int_0^\frac{x}{\sqrt{1-x^2}} s {\hat{V}}(s) \frac{d {\hat{V}}}{d s}(s) ds < 0, \quad x \in [0,1), \tag{110} \end{align}\] and \(H(x)\) satisfies \[\label{defh0} 0 < H(0) = \int_0^1 \frac{1-t^2}{t^2} \bigg(\displaystyle \int_0^\frac{t}{\sqrt{1-t^2}} s \hat{V}^2(s) ds \bigg) \, dt < \frac{V_\infty^2}{2} \, ,\tag{111}\] \(H(1) = 0\), \(\lim_{x\to0} \frac{dH}{dx}= 0\) and \(\lim_{x\to0} \frac{d^2H}{dx^2} = 0\). Hence, \(H(x)\) is decreasing and concave. A sketch of \(H(x)\) is presented in Figure [fig:H40x41]. If we define \[\label{defalpha} \frac{dH}{dx} (1) = - \lim_{\xi \to \infty} \frac{1}{\xi^2} \int_0^\xi s {\hat{V}}^2 (s) ds = : - \beta < 0\tag{112}\] we observe that due to the concavity of \(H(x)\) we have the bounds \[H(0) (1 - x) \le H(x) \le \beta (1 -x) \qquad x \in [0,1]\] with \(H(0)\), \(\beta\) defined in 111 , 112 and satisfying \(H(0) < \beta\).

As a consequence, the function \(F = H + E_0 (1-x)\) in the Ricatti-type equation 107 obeys the bounds \[\label{boundH} (H(0) + E_0) (1-x) \le F(x) = H(x) + E_0 (1-x) \le (\beta + E_0) (1-x)\tag{113}\] and resembles one of the following configurations, see Fig. [fig:F40x41] :

1. \(F\) in negative on \([ 0,1]\), which occurs when \(\beta + E_0 < 0\)

2. \(F\) is first negative, it becomes zero at a point \(0<x_1<1\) and then \(F\) is positive. This occurs when \(H(0) + E_0 < 0 < \beta + E_0\).

3. \(F\) in positive for \([0,1)\), occurring when \(0 < H(0) + E_0\).

This restricts considerably the possible shapes of \(\Theta\). A computation using 107 gives \[\nu \, \frac{d^2 \Theta}{dx^2} (0) = - (H(0) + E_0 ) = - F(0) \, .\] Since \(\Theta(0) = \frac{d\Theta}{dx} (0) =0\) the value of \(F(0)\) determines in which half-plane \(\Theta (x)\) initially starts. Consider first the case where \(F(0) = H(0) + E_0>0\). Then \(F(x)>0\) for all \(x\in(0,1)\), that is the setting of case \({(c)}\), illustrated in Figure [fig:F95c]. \(\Theta(x)\) starts in the lower half-plane. If it crosses the x-axis at some point \(x_1<1\), then \[0 \le \nu \, \frac{d \Theta}{dx} (x_1) = - \frac{x_1}{(1-x_1^2)^2} \, {F}(x_1)\] which contradicts that \(F(x)>0\). Therefore, when \(H(0) + \, {E_0} >0\), the function \(\Theta(x)\) remains negative for all \(x\in (0,1)\).

Next, consider the case \(H(0) + E_0<0\) which corresponds to the cases \({(a)}\) or \({(b)}\). We now have \[\begin{align} \nu \, \frac{d^2 \Theta}{dx^2} (0) &= - \Big ( H(0) + \, {E_0} \Big ) > 0 \, \end{align}\] and \(\Theta(x)\) starts at the upper half-plane. If \(\Theta\) crosses the x-axis at some point \(x_1<1\) then \[\label{eq5} 0 \ge \, \nu \, \frac{d \Theta}{dx} (x_1) = - \frac{x_1}{(1-x_1^2)^2}\, {F}(x_1)\tag{114}\] It is possible to cross going downwards but it is not possible to cross a second time going upward again. We conclude that when case \((b)\) happens then \(\Theta(x)\) can possibly cross the axis but it cannot cross a second time. By contrast when case \((a)\) happens, 114 contradicts \(F(x)<0\), \(x\in (0,1)\), i.e. \(\Theta(x)\) cannot cross the axis. In this regime, the differential inequality \[\nu \, \frac{d^2 \Theta}{dx^2} > \frac{1}{2} \Theta^2, \quad x \in [0,1)\] implies that the solution \(\Theta(x)\) blows up in finite time, and if \(\nu\) is sufficiently small, this will happen within the interval \([0,1]\).

Summarizing there are the following possibilities: When \(F(x) > 0\) as in Fig. [fig:F95c] then \(\Theta (x) < 0\); when \(F(x) < 0\) as in Fig. 7 then \(\Theta (x) > 0\). When \(F(x)\) changes sign as in Fig. [fig:F95b], then \(\Theta\) starts at the upper half plane, and either \(\Theta\) stays there afterwards or it might cross to the lower half plane. Accordingly, there are the following configurations: Either \(\Theta (x) < 0\) which is called Zone A; or it starts with \(\Theta (x) >0\) but crosses to the negative half-plane and stays there thereafter, called Zone B; or it starts and stays \(\Theta (x) >0\), called Zone C. In Zone C the solution \(\Theta\) lies in the upper half-plane and if \(\nu\) is sufficiently small it will blow up before reaching \(x=1\).

A-priori estimates for \(\Theta(x)\) are derived below.

When the bounds of Lemma 4 hold, the solution \(\Theta\) of 107 cannot blow up in \((0,1)\). The hypothesis \(\beta + E_0 > 0\) and thus item (ii) of Lemma 4 is of theoretical interest, as the parameter \(\beta\) cannot be a-priori determined. Nevertheless, using the bound \(H(x) < \frac{1}{2} V_\infty^2\), we can prove a variant of (ii) providing a bound from below.

5.4 Convergence as \(\nu \to 0\)↩︎

Let \(\big\{ (\Theta_\nu, V_\nu) \big\}_{\nu >0}\) be a family of solutions of 100 , \[\label{eql} \begin{align} \nu \, \frac{d \Theta_\nu}{dx} &= \tfrac{1}{2} \Theta_\nu^2 - \mathcal{F}\left(V_\nu ; x\right), \\ \nu \frac{d^2 {V_\nu}}{d^2x} &= \Theta_\nu \, \frac{d {V_\nu}}{dx} , \end{align}\tag{116}\] with \(\nu > 0\) and we are interested in the limit \(\nu \to 0\). We assume that \(\kappa := 2E_0+ V_\infty^2>0\) and that the family \((\Theta_\nu, V_\nu)\) satisfies the uniform bound \[\label{hypothesis} -\sqrt{\kappa} \, \frac{\sqrt{x(1-x)}}{1-x^2} < \Theta_\nu(x) < 0, \quad \forall \,\, x \in (0,1).\qquad{(1)}\] The assumption \(\kappa > 0\) is natural since only then solutions of stationary self-similar axisymmetric Euler equations exist, see sections 3 and 4.

Regarding the uniform bounds ?? : The left hand inequality is guaranteed by Corollary 2 while the right hand inequality follows from Lemma 4 in the more restrictive range \(E_0>0\). The reader should note, that, as expected from numerical computations in the parameter range \(\kappa > 0\), the solution of 116 enters Zone A - the region that \(\Theta(x) < 0\) - as \(\nu\) decreases and stays there. This is corroborated by the bifurcation diagram sketched in Figure 11.

Without loss of generality, we will examine the case \(V_\infty>0\). Then 102 implies \[\begin{align} \label{dV} \frac{d V_\nu}{dx} &= V_\infty \frac{e^\mathlarger{{\frac{1}{\nu} \int_0^x \Theta_\nu(s) ds}}}{\int_0^1 e^{\frac{1}{\nu} \int_0^t \Theta_\nu(s) ds} dt} \, > 0 \\ 0 &< V_{\nu} (x) < V_\infty \, , \quad x \in (0,1) \, . \end{align}\tag{117}\] The sequence \(\{ V_\nu \}\) is of bounded variation and by Helly’s theorem it admits a convergent subsequence (which will be still denoted by \(V_\nu\)) such that \[\label{convV} V_{\nu}(x) \rightarrow V(x) \quad a.e \; \; x \in (0,1) \, .\tag{118}\]

For the sequence \(\{ \Theta_\nu \}\) the uniform bound imply that for \(p \in [1,2)\) \[\int_0^1 |\Theta_\nu (x)|^p dx \le C \int_0^1 (1-x)^{-\frac{p}{2}} dx = : K < \infty\] Using weak convergence techniques [18], there exists a subsequence (still denoted by \(\Theta_\nu\)) and a function \(\Theta \in L^p (0,1)\) such that \[\label{wkconvTheta} \begin{align} \Theta_\nu &\rightharpoonup \Theta \qquadweakly in L^p(0,1) \\ that is \qquad \int \Theta_\nu \psi dx &\to \int \Theta \psi dx \qquad for \psi \in L^{p^\prime}(0,1) \end{align}\tag{119}\] with \(p^\prime\) the dual exponent of \(p\). Setting \(\psi (y) = \chi_{[0,x]} (y)\), the characteristic function of \([0,x]\), we deduce \[\label{limit-G} G_\nu(x) = \int_0^x \Theta_\nu (s) ds \rightarrow G(x) := \int_0^x \Theta(s) \, ds \quad\tag{120}\]

Next, we employ ideas developed in the context of zero-viscosity limits for Riemann problems of conservation laws [15], [16], [21]. We set \[\label{measure} \mathop{\mathlarger{\pi_\nu \mathrel{\vcenter{:}}= \frac{d \, V_\nu}{dx}= V_\infty \, \, \frac{e^{\frac{1}{\nu} G_\nu(x)}}{\int_0^1 e^{\frac{1}{\nu} G_\nu(t)} dt}}},\tag{121}\] where \(\mathlarger{G_\nu(x) = \int_0^x \Theta_\nu(s) \, ds}\) and view \(\{ \pi_\nu \}\) as a sequence of probability measures. By 121 , the sequence \(\{ \pi_\nu \}\) is uniformly bounded in measures and there exists a subsequence \(\pi_\nu\) and a probability measure \(\pi\) such that \[\label{convprob} \pi_\nu \stackrel{\ast}{\rightharpoonup} \pi, \quad \textrm{in measures } \quad \mathcal{M } [0,1] \, .\tag{122}\] Due to the correspondence between functions of bounded variations and Borel measures [19] one sees that the distribution function of the measure \(\pi\) is precisely the function \(V(x+)\) and 122 reflects 118 .

We are now ready to state the first convergence theorem.

The inequality ?? is due to the fact that only weak convergence for \(\Theta\) is available under ?? . It can be improved if we have pointwise convergence for \(\Theta\), namely if \[\label{hypconvergence} \Theta_\nu \to \Theta \quad a.e., \,\, as \quad \nu\to 0. \tag{129}\] We will later justify 129 under the hypothesis \(E_0 > 0\). Using 129 , one may pass to the limit in 107 in the sense of distributions and deduce \[\label{limit-Th} \frac{\Theta^2 (x)}{2} = \mathcal{F}(V,x).\tag{130}\] Let us compute \(\Theta(x)\) for Case \(2\). Since \(V(x)\) is given by 128 , 130 yields \[\begin{align} \frac{\Theta^2 (x)}{2} = \mathcal{F}(V_\infty,x) = \big(\frac{1}{2}V_\infty^2 + \, {E_0}\big) \,\, \frac{x(1 - x)}{(1-x^2)^2} \end{align}\] and thus, since \(\Theta(x)<0\), \[\Theta(x) = -\sqrt{V_\infty+2E_0}\,\, \frac{\sqrt{x(1 - x)}}{(1-x^2)} \, .\]

For Case \(1\), we recall that \(V(x)\) is given by 127 and \(\Theta(x)=0\) a.e. for \(x\in[0,a]\), \(a>0\). Then 130 with 105 imply \[\begin{align} E_0 \,\, (1-x) = - \int_x^1 \frac{1-t^2}{t^2} \bigg(\displaystyle \int_0^t \frac{\sigma}{(1-\sigma^2)^2} V^2(\sigma) d\sigma \bigg) \, dt \, , \quad for x \in [0,a]. \end{align}\] If we differentiate this, we get \[\begin{align} E_0 \,\, \frac{x^2}{1-x^2} = - \displaystyle \int_0^x \frac{\sigma}{(1-\sigma^2)^2} V^2(\sigma) d\sigma. \end{align}\] Differentiating once more yields \[\begin{align} V^2(x) = - 2\,\, E_0 \quad for \,\, x\in(0,a) \end{align}\] which contradicts the assumption \(supp \, \pi = [0,a]\) with \(a> 0\) . We conclude that only Case 2 can happen and \((\Theta_\nu, V_\nu)\) converges almost everywhere as \(\nu\to0\) to a solution \((\Theta, V)\) of the Euler equations. This provides a criterion to select the type of Euler solution that occurs and clearly only solutions with solutions with \(\Theta<0\) are admissible.

In conclusion, we have the following theorem

If the solution takes values in the range where \(\Theta(x) < 0\) then it lies in Zone A. The analysis of section 5.3 shows that \(F\) in 113 then takes values \(F(x) > 0\) on \((0,1)\). The reader should note that numerical computations suggest that as \(\nu\) decreases the solution of 100 enters Zone A, that is the region that \(\Theta(x) < 0\), and stays there. No oscillations in \(\Theta\) are observed numerically.

We next restrict in the range \(E_0 > 0\) and justify 129 .

Combining Proposition 8 with Theorem 5 gives

Note that when \(V_\infty > 0\) among the two solutions of the Euler equations 41 the one selected at the zero-viscosity limit corresponds to the negative sign, see Fig. 2 (b). One easily checks that when \(V_\infty < 0\) again the negative sign is selected.

6 Boundary layer analysis for a model problem↩︎

In this section, we investigate the asymptotic behavior of solutions of system 100 as \(\nu\to0\). We consider a model problem which is a simplification of the initial equations, with the objective to understand the boundary layer. For small viscosities, \(V(x)\) is approximated by setting \(V(x) \equiv V_\infty\), leading to \[\mathcal{F}(V,x) = \mathcal{F}(V_\infty,x) = \big(\frac{V_\infty^2}{2} + \, {E_0}\big) \frac{x - x^2}{(1-x^2)^2}.\] This reduces system 100 to a simpler form which will be referred to as the model problem, namely, \[\tag{135} \begin{align} \tag{136} \frac{\bar{\Theta}^2 (x)}{2} - \nu \, \frac{d \bar{\Theta}}{dx} (x) &= K \frac{x - x^2}{(1-x^2)^2},\\ \tag{137} \nu \frac{d^2 \bar{V} }{d^2x} & = \bar{\Theta} \,\, \frac{d \bar{V}}{dx}, \\[5pt] \bar{\Theta}(0)=0,\quad \bar{V}(0)&=0, \quad \bar{V}(1)=V_\infty, \end{align}\] where \(K= \frac{V_\infty^2}{2} + \, {E_0}>0\). Note that equation 136 for \(\bar{\Theta}\) is now independent of \(\bar{V}\), while 137 is still coupled.

We use the method of matched asymptotic expansions from singular perturbation theory. The method aims to to construct an asymptotic approximation of the solution inside the boundary layer and a solution valid away from the boundary layer, and then combine them through a matching process. We refer to solutions within the boundary layer as inner solutions, and to solutions away of the layer as the outer solutions, [22], [23].

We consider first the equation 136 for \(\bar{\Theta}\) \[\begin{align} \tag{138} \frac{\bar{\Theta}^2 (x)}{2} - \nu \, \frac{d \bar{\Theta}}{dx} (x) &= K \frac{x - x^2}{(1-x^2)^2},\\ \tag{139} \bar{\Theta}(0)&=0, \end{align}\] and apply the method of matched asymptotic expansions. To construct the outer solution, assume that \(\bar{\Theta}\) can be written as a power series with powers of \(\nu\), \[\bar{\Theta}(x) \approx \bar{\Theta}_0(x) + \nu \bar{\Theta}_1(x) + \mathcal{O}(\nu^2),\] and substitute it back to 138 . If we focus on the leading terms, i.e terms of order \(\nu^0\), we obtain the equation \[\bar{\Theta}_0(x) = \pm \sqrt{ 2 \, K \frac{x - x^2}{(1-x^2)^2}}.\] The choice for the sign \(\bar{\Theta}_0\) will depend on \(\bar{\Theta}\). Recall that \(\bar{\Theta}\) solve the equation 138 , we have \[\bar{\Theta}(0) = 0, \quad \frac{d \bar{\Theta}}{dx}(0) = 0, \quad \textrm{and} \quad \frac{d^2 \bar{\Theta}}{dx^2}(0)<0,\] which implies \(\bar{\Theta}\) should be negative for all \(0<x<1\). Hence, we choose \(\bar{\Theta}_0\) to be negative. Note that the boundary condition at \(x=0\) is automatically satisfied.

We proceed now with the inner solution. Expecting the boundary layer to locate at \(x=0\), we introduce the stretched variable \(\eta = \nu^{-\frac{2}{3}} \, x\). Setting \(\Phi(\eta) = \nu^{-\frac{1}{3}}{\bar{\Theta}(x)}\), the problem 138 139 takes the form \[\begin{align} \label{asy-Phi} & \frac{\Phi^2}{2} - \frac{d \Phi}{d\eta} = K \, \frac{\eta}{(1-\nu^{\frac{2}{3}} \,\eta) (1+\nu^{\frac{2}{3}} \, \eta)^2} \\ & \Phi(0)=0, \end{align}\tag{140}\] If the approximation of \(\Phi\) in powers of \(\nu\) is given by \[\Phi(\eta) \approx \Phi_0(\eta) + \nu \Phi_1(\eta) + \mathcal{O}(\nu^2),\] and using the Taylor expansion of the right-hand side of 140 , \[\frac{\eta}{(1-\nu^{\frac{2}{3}} \,\eta) (1+\nu^{\frac{2}{3}} \, \eta)^2} \approx \eta + \mathcal{O}(\nu^\frac{2}{3}),\] we obtain for the leading term \[\tag{141} \begin{align} \tag{142} & \frac{\Phi_0^2}{2} - \frac{d \Phi_0}{d\eta} = K \, \eta, \\ \tag{143} & \Phi_0(0)=0. \end{align}\]

Inspired by [24], equation 142 can be transformed into the well-known Airy equation \[y''(t) - t \,y(t)=0.\] Using the transformation \[\label{Airy} \Phi_0(\eta) = -\frac{2}{\mathcal{U}(\eta)} \frac{d \mathcal{U}}{d\eta},\tag{144}\] equation 142 reduces to \[\frac{d^2 \mathcal{U}}{d\eta^2} = \frac{K}{2} \eta \, \mathcal{U}(\eta).\] and its solutions are given as a linear combination of special functions \(Ai\), the Airy function of the first kind, and \(Bi\), the Airy function of the second kind, see [23]. Namely, \[\mathcal{U}(\eta) = \text{Ai}\left(\big(\frac{K}{2}\big)^{1/3} \eta \right) c_1 +\text{Bi} \left(\big(\frac{K}{2}\big)^{1/3} \eta \right) c_2,\] where \(c_1, c_2\) are integration constants. Using 144 , we express \(\Phi_0\) as a linear combination of Airy functions \(Ai, Bi\) and their derivatives \(Ai', Bi'\). Hence, \(\Phi_0\) takes the form \[\Phi_0(\eta) = -2 \bigg(\frac{K}{2}\bigg)^{1/3} \, \frac{\text{Ai}'\left(\big(\frac{K}{2}\big)^{1/3} \eta \right) c_1 +\text{Bi}'\left(\big(\frac{K}{2}\big)^{1/3} \eta \right) c_2}{\text{Ai}\left(\big(\frac{K}{2}\big)^{1/3} \eta \right) c_1 +\text{Bi} \left(\big(\frac{K}{2}\big)^{1/3} \eta \right) c_2},\] or equivalently, \[\Phi_0(\eta) = -2 \bigg(\frac{K}{2}\bigg)^{1/3} \, \frac{\text{Ai}'\left(\big(\frac{K}{2}\big)^{1/3} \eta \right) C +\text{Bi}'\left(\big(\frac{K}{2}\big)^{1/3} \eta \right)}{\text{Ai}\left(\big(\frac{K}{2}\big)^{1/3} \eta \right) C +\text{Bi} \left(\big(\frac{K}{2}\big)^{1/3} \eta \right) },\] where \(C= \frac{c_1}{c_2}\) is a constant. Imposing now the boundary condition 143 , we get \[\begin{align} 0 = \frac{\text{Ai}'\left(0\right) C +\text{Bi}'\left(0 \right)}{\text{Ai}\left(0\right) C +\text{Bi} \left(0 \right)} = \frac{-\frac{C}{3^{1/3} \, \Gamma \left(\frac{1}{3}\right)} + \frac{3^{1/6}}{\Gamma \left(\frac{1}{3}\right)}}{\frac{C}{3^{2/3} \Gamma \left(\frac{2}{3}\right)}+\frac{1}{3^{1/6} \Gamma \left(\frac{2}{3}\right)}}= \frac{3^{1/3} \, \Gamma(\frac{2}{3}) (-C+ \sqrt{3})}{\Gamma(\frac{1}{3})(C+ \sqrt{3})}, \end{align}\] where \(\Gamma\) denotes the Gamma function. Thus, the constant \(C\) is determined as \[C = \sqrt{3}.\] Consequently, the leading term \(\Phi_0\) of the inner solution becomes \[\Phi_0(\eta) = -2 \bigg(\frac{K}{2}\bigg)^{1/3} \frac{\sqrt{3}\, \text{Ai}'\left(\big(\frac{K}{2}\big)^{1/3} \eta \right)+\text{Bi}'\left(\big(\frac{K}{2}\big)^{1/3} \eta \right)}{\sqrt{3} \, \text{Ai}\left(\big(\frac{K}{2}\big)^{1/3} \eta \right)+\text{Bi} \left(\big(\frac{K}{2}\big)^{1/3} \eta \right)}.\]

The last step of the method of matched asymptotic expansions is to derive a uniform expansion of the solution of 138 over the whole domain \([0,1)\). Combining the approximations of inner and outer solution, we conclude that the asymptotic expansion of \(\bar{\Theta}\) as \(\nu\to0\) is \[\label{asymd-theta} \bar{\Theta}(x) \approx \left\{ \begin{align} &-\sqrt{ 2 \, K \frac{x - x^2}{(1-x^2)^2}}, \quad && A\nu^{\frac{2}{3}}<x<1\\ &\Phi_0(\nu^{-\frac{2}{3}} \, x), \quad && 0\leq x \leq A \nu^{\frac{2}{3}} \end{align} \right.\tag{145}\] where \(A\) is a positive constant and \[\Phi_0(\nu^{-\frac{2}{3}} \, x) = - 2 \nu^{{1}/{3}} \bigg(\frac{K}{2}\bigg)^{1/3} \frac{\sqrt{3} \text{Ai}'\left(\big(\frac{K}{2}\big)^{1/3} \nu^{-\frac{2}{3}} \, x \right)+\text{Bi}'\left(\big(\frac{K}{2}\big)^{1/3} \nu^{-\frac{2}{3}} \, x \right)}{\sqrt{3} \text{Ai}\left(\big(\frac{K}{2}\big)^{1/3} \nu^{-\frac{2}{3}} \, x \right)+\text{Bi} \left(\big(\frac{K}{2}\big)^{1/3} \nu^{-\frac{2}{3}} \, x \right)}.\] Moreover, the boundary layer is formed at \(x=0\) and its size is of order \(\nu^{\frac{2}{3}}\).

Let us consider now the equation 137 for \(\bar{V}\). Motivated by the asymptotic expansion of \(\bar{\Theta}\), we introduce the following simplified problem \[\tag{146} \begin{align} \tag{147} & \nu \frac{d^2 \bar{V} }{d^2x}(x) = -\sqrt{2\,K\;\frac{x - x^2}{(1-x^2)^2}} \,\, \frac{d \bar{V}}{dx} (x), \\ \tag{148} & \bar{V}(0)=0, \quad \bar{V}(1)=V_\infty, \end{align}\] and apply the method of matched asymptotic expansions. As before, we expect the boundary layer to be located at \(x=0\).

First, we derive the approximation of the outer solution. Assuming that \(\bar{V}\) can be expressed as a power series with powers of \(\nu\), i.e. \[\bar{V}(x) \approx \bar{V}_0(x) + \nu \bar{V}_1(x) + \mathcal{O}(\nu^2),\] the equation 147 yields the following differential equation for the lead term \(\bar{V}_0\) \[-\sqrt{2\,K\,\frac{x - x^2}{(1-x^2)^2}} \,\, \frac{d \bar{V}_0}{dx} =0 \quad \textrm{which implies} \quad \bar{V}_0 (x) \equiv c,\] where \(c\) is a constant. To determine this constant, we consider the boundary condition away from the boundary layer, i.e at \(x=1\) and thus, 148 implies \[\bar{V}_0(x) \equiv V_\infty.\] For the inner solution we introduce the variable \(\eta = \nu^{-2/3}\,x\) and set \(\bar{V}(x) = Y(\eta)\). The problem 146 reduces to \[\begin{align} \label{asy-Y} \frac{d^2 Y}{d^2 \eta} + \sqrt{2\,K} \, \sqrt{\frac{\eta}{(1-\nu^{\frac{2}{3}} \,\eta) (1+\nu^{\frac{2}{3}} \, \eta)^2}} \,\, \frac{d Y}{d\eta} = 0 \quad \textrm{with} \quad Y(0) =0. \end{align}\tag{149}\] If the approximation of \(Y\) in powers of \(\nu\) is expressed as \[Y(\eta) \approx Y_0(\eta) + \nu Y_1(\eta) + \mathcal{O}(\nu^2),\] and the Taylor expansion of the right-hand side of 149 is \[\sqrt{\frac{\eta}{(1-\nu^{\frac{2}{3}} \,\eta) (1+\nu^{\frac{2}{3}} \, \eta)^2}} \approx \sqrt{\eta} + \mathcal{O}(\nu^\frac{2}{3}),\] then we have to solve the following problem for leading order term \(Y_0\) \[\begin{align} \frac{d^2 Y_0}{d^2 \eta} + \sqrt{2\,K\,\eta}\, \frac{d Y_0}{d\eta} &= 0,\\ Y_0(0) &= 0. \end{align}\] Its solution is \[\begin{align} Y_0(\eta) = C \int_0^\eta \mathlarger{\mathlarger{e^ {-\frac{2 \sqrt{2\,K}}{3} \, s^{3/2}}}} ds, \end{align}\] where \(C\) is an integration constant. Since both boundary conditions 148 been have taken into consideration, the unknown \(C\) is determined by matching the outer and the inner expansions. Considering that both the inner and outer solutions approximate the same function in different regions, we impose that the two solutions are equal in a transition area close to the boundary layer, [23]. Therefore, we require \[\lim_{x\to 0} \bar{V}_0(x) = \lim_{\eta \to \infty} Y_0(\eta),\] that implies \[\begin{align} C = \mathop{ \frac{V_\infty}{\int_0^\infty \mathlarger{\mathlarger{e^ {-\frac{2 \sqrt{2\,K}}{3} \, s^{3/2}}}} ds}}, \end{align}\] and completes the derivation of the inner solution.

As the last step of the method of matched asymptotic expansions, we derive a uniform expansion of the solution of 146 over the whole domain \([0,1)\). Therefore, as \(\nu\to0\) the asymptotic expansion of \(\bar{V}(x)\) takes the form \[\begin{align} \label{V-asympt} \bar{V}(x) \approx {V_\infty} \,\, \mathlarger{\frac{\bigints_0^{\nu^{-2/3}\,x} e^ {-\frac{2 \sqrt{2\,K}}{3} \, s^{3/2}} ds}{\bigints_0^\infty e^ {-\frac{2 \sqrt{2\,K}}{3} \, s^{3/2}} ds}} + \mathcal{O}(\nu). \end{align}\tag{150}\] and the boundary layer at \(x=0\) is of the same size as for \(\bar{\Theta}(x)\).

7 Stationary Navier-Stokes - Numerical Results↩︎

In this section we construct a numerical scheme to solve the the coupled system of ordinary differential equations 92 using an iterative algorithm. Moreover, we illustrate some numerical experiments for different values of the parameters \(\nu, \;E_0\) and \(V_\infty\).

Let us recall the system 92 . For convenience we use the equivalent formulation of the problem, system 101 104 in variable \(x\), where \(0<x<1\) \[\begin{align} \tag{151} \frac{\Theta^2 (x)}{2} - \nu \, \frac{d \Theta}{dx} (x) &= \mathcal{F}(V,x),\\ \tag{152} \nu \frac{d^2 V}{d^2x}(x) &= \Theta(x) \, \frac{d V}{dx} (x), \end{align}\] equipped with the boundary conditions \[\begin{align} \tag{153} \Theta = V = 0 , \,\, &\textrm{ at } \quad x = 0, \\ \tag{154} V \to V_\infty , \,\, &\textrm{ as} \quad x \to 1, \end{align}\] where the functional \(\mathcal{F}\left(V,x\right)\) is defined by \[\begin{align} \label{defFE} \mathcal{F}\left(x\right) = \mathcal{F}\left(V,x; E_0 \right) & = \frac{x}{(1-x^2)^2}\Bigg[\int_x^1 \frac{1-t^2}{t^2} \bigg(\displaystyle \int_0^t \frac{\sigma}{(1-\sigma^2)^2} V^2(\sigma) d\sigma \bigg) \, dt + \, {E_0} \,\, (1 - x)\Bigg]. \end{align}\tag{155}\]

7.1 Discretization↩︎

To obtain numerical approximations to system 151 154 we use an iterative algorithm: Given \(V\) we solve numerically the initial value problem 151 by a Runge-Kutta method and update \(\Theta\) which is then used to solve the two-point boundary value problem 152 using finite differences. This process is repeated until convergence. The algorithm is described in detail in Section 7.2.

The discretization of 152 is based on finite differences. On \([0,1]\) we introduce a uniform mesh of fixed width \(\Delta x\). Given \(N\in\mathbb{N}\), set \(\Delta x = \frac{1}{N}\) and define the discrete points \(x_j = j\Delta x\), \(j= 0,1,2,\dotsc, N\). We denote by \(\Theta^j \approx \Theta(x_j)\) and \(V^j \approx V(x_j)\) and use central finite differences to discretize 152 and obtain \[\label{FM1} \begin{align} \nu \,\frac{V^{j+1} - 2V^j + V^{j-1}}{\Delta x^2} & = \Theta^j \, \frac{V^{j+1} - V^{j-1}}{2\Delta x}, \quad j=1,\dotsc, N-1, \\ V^0 = 0, & \;V^N = V_{\infty}, \end{align}\tag{156}\] or equivalently \[\begin{align} \label{FM2} V^{j-1} \bigg(2\nu + {\Delta x} \,\Theta^j \bigg) - 4 \nu\, V^j & + V^{j+1} \bigg(2 \nu - {\Delta x} \, \Theta^j \bigg) = 0 , \qquad j=1,\dotsc, N-1, \\ V^0 = 0, & \;V^N = V_{\infty}, \end{align}\tag{157}\] which forms a tri-diagonal linear system solved by a direct method.

To discretize the nonlinear initial value problem 151 we use a Runge-Kutta method. The Runge-Kutta methods are multistage methods that compute approximations to the solution at intermediate points which are later combined to advance the solution at the next discretization point. Equation 151 can be rewritten as \[\nu \frac{d \Theta}{dx} (x) = \frac{\Theta^2 (x)}{2} - \mathcal{F}(x).\] For simplification of the presentation, we denote the right hand side of the above relation as follows \(f(x, V, \Theta) = \tfrac{1}{2} \Theta^2 - \mathcal{F}(x)\). The Runge-Kutta method for 151 is obtained by discretizing the derivative using the following method:

1. Set \(\Theta^{j,1} =\Theta^j\), where \(\Theta^j\) is the approximate solution at grid point \(x_j\).

2. Compute the \(\mathcal{K}\) intermediate stages for \(m =1,\dots,\mathcal{K}\) \[\label{RK1} \nu \Theta^{j,m} = \nu \Theta^j + {\Delta x} \sum_{s=1}^{\mathcal{K}} \alpha_{m,s} f \left(x_j+c_s {\Delta x},V^{j,s},\Theta^{j,s}\right),\tag{158}\] where \(V^{j,s} \approx V(x_j+c_s \Delta x)\) and \(\Theta^{j,s} \approx \Theta(x_j+c_s \Delta x)\). The coefficients \(\alpha_{m,s}\) are constants with sum for every row equals to \(c_m\), i.e. \[\sum_{s=1}^{\mathcal{K}} \alpha_{m,s} = c_m, \quad m =1,\dots,\mathcal{K}.\]

3. Set \[\nu\Theta^{j+1} = \nu \Theta^j + {\Delta x} \sum_{s=1}^{\mathcal{K}} \beta_s f \left(x_j+c_s {\Delta x},V^{j,s},\Theta^{j,s}\right)\] the solution at the next point. Here the coefficients \(\beta_s\) are constants with total sum equal to 1 and values of \(V^{j,s}\) are computed using interpolation.

The coefficients \(\alpha_s, \beta_s\) and \(c_s\) of the associated with the Runge-Kutta method are usually presented in a matrix form referred as Butcher tableau. That is, \[\renewcommand\arraystretch{1.2} \begin{array} {c|cccc} c_1 & \alpha_{1,1} &\alpha_{1,2} &\cdots &\alpha_{1,\mathcal{K}}\\ c_2 & \alpha_{2,1} &\alpha_{2,2} &\cdots &\alpha_{2,\mathcal{K}} \\ \vdots & \vdots & \vdots & & \vdots\\ c_\mathcal{K} & \alpha_{\mathcal{K},1}& \alpha_{\mathcal{K},2} &\cdots &\alpha_{\mathcal{K},\mathcal{K}}\\ \hline & \beta_1 &\beta_2 &\cdots &\beta_\mathcal{K} \end{array}\] In our numerical experiments we use the well known Runge-Kutta-Fehlberg (RKF4(5)) method which allows adaptive stepsize control. For small values of \(\nu\) the solution \(\Theta\) of the i.v.p 151 might blow up, see the discussion in Section 5.3. Furthermore, as \(x\to 1\) we expect the system to be singular since the point \(x=1\) corresponds to the vortex line. In that respect, choosing the stepsize adaptively allows us to capture the correct behaviour of the solution close to the singularity.

7.2 Implementation Details↩︎

To compute the numerical approximation of the solution for the coupled system 151 - 152 , we construct an iterative algorithm and each equation is solved separately using information from the previous iterations. In particular, we follow the following algorithm:

1. Set \(V_0^j \equiv V_\infty\), for all \(j\)

2. For every iterative step \(i=1,\dotsc\):

   1. Given \(V_{i-1}\), compute \(\Theta_i^j\), \(\forall \, j\), using the RKF-method in \([0,1]\)

   2. Given \(\Theta_i\), compute \(V_i^j\), \(\forall \, j\), by solving the linear system 157

   3. Solution \((\Theta,V)\) is obtained when the error between two consecutive iterations is small, i.e. \[\bigg\{\Delta x \sum_j |V_{i+1}^j - V_i^j | ^2 \bigg\}^\frac{1}{2} < \epsilon \quad \textrm{and}\quad \bigg\{\Delta x \sum_j |\Theta_{i+1}^j - \Theta_i^j | ^2 \bigg\}^\frac{1}{2} < \epsilon\] for some \(\epsilon>0\) small.

Note that for the approximation of the integrals in \(\mathcal{F}\) we use the composite Simpson’s rule.

7.3 Numerical Tests↩︎

In this section we exhibit numerical approximation of the solution for 151 - 152 . We illustrate the results for a different combinations of parameters \(\nu\) and \(E_0\) while the parameter \(V_\infty\) is fixed. For the following examples we take \(V_\infty=1\).

Consider first the case for \(\nu = 0.2\) and \(E_0 = -0.5\). Then \(\Theta(x)\) is positive, and the flow is directed inward near the plane \(z = 0\) and upward near the vortex line, see Figure 8.

The converse behavior, i.e. \(\Theta\) is negative, occurs for \(\nu = 0.05\) and \(E_0 = 0.15\). The flow now has the reverse direction, it is directed outward near the plane \(z = 0\) and downward near the vortex line, see Figure 9.

For the parameters \(\nu = 0.1\) and \(E_0 = -0.2\), the stream function \(\Theta\) is first positive and then becomes negative; the flow is directed inward near the plane \(z = 0\) and downward near the vortex line, see Figure 10.

The solid lines in Figures 9, 8 and 10 are the contour lines of the streamfunction \(\Theta\).

7.4 Bifurcation diagram↩︎

We present now a \(E_0 - \nu\) bifurcation diagram for 100 . This diagram is computed using the methodology detailed in Section 7.2. To take into account the effect of the parameter \(V_\infty\) we consider the following scaled variables: \[\phi = \frac{V}{V_\infty} \qquad \vartheta = \frac{\Theta}{V_\infty}, \qquad \mu = \frac{\nu}{V_\infty}, \qquad p_0 = \frac{E_0}{V_\infty^2}\] Then system 151 154 becomes \[\label{scaledsys} \begin{align} \mu \frac{d\vartheta}{dx} & = \frac{1}{2}\vartheta^2 - \mathcal{F}(x, \phi, p_0), \\ \mu \frac{d^2 \phi}{dx^2} & = \vartheta\frac{d\phi}{dx}, \\ \phi = \vartheta = 0 & \;\text{ at } \;x=0, \\ \phi \to 1, & \;\text{ as } \;x\to 1 , \end{align}\tag{159}\] where \(\mathcal{F}\) is given by 155 . The new scaled system exhibits the same behaviour as the original one, so we proceed in identifying numerically the four different zones, see Section 5.3:

• Zone A: \(\vartheta\) is negative

• Zone B: \(\vartheta\) starts positive and then changes sign

• Zone C: \(\vartheta\) is positive

• No Solution: \(\vartheta\) is positive and for small \(\mu\) blows up

We take \(V_\infty=1\) and we consider values of \((E_0, \nu) \leftrightarrow (p_0, \;\mu)\) in \(\mathcal{B}= [-2, 2]\times[4\cdot10^{-5}, 0.6]\). The set \(\mathcal{B}\) is covered initially by \(256=16\times 16\) patches \(\mathcal{B}_{k,\ell}\) each of size \(\delta p_0\times \delta\mu\). Along the horizontal axis, we take a uniform size \(\delta p_0 = 0.25\), while \(\delta\mu\) is non-uniform, with a finer grid around \(\mu=4\cdot10^{-5}\) and gradually increasing towards \(\mu=0.6\), with an average \(\delta\mu \sim 3.75\cdot 10^{-2}\). In each \(\mathcal{B}_{k,\ell}\) we consider further a \(51\times 51\) uniform grid of values \(\{{p_0}_{i,j}^{k,\ell}, \mu_{i,j}^{k,\ell}\}, \;i,j=1,\dotsc,51\). For all patches and for all values in the patch we identify in which zone the solution belongs to, by solving numerically the scaled system 159 . The total computational cost of such process is rather small since the work in each patch can be computed independently.

The results of this process are depicted in Figure 11. On the left graph, the four zones are clearly marked by their boundaries. To distinguish further each zone we “separate" them by shifting slightly each zone along the horizontal axis and Zone C along the vertical axis. On the right graph, the full structure of each zone is now revealed. The boundaries of all zones meet at the point \((p_0, \mu) = (E_0, \nu) = (-\frac{1}{2}, 4\cdot10^{-5})=(-\beta,4\cdot10^{-5})\), where \(\beta = \frac{dH}{dx}(1)\) is the constant appearing in 112 . For any \(E_0 < -\beta\) on the line \(\mu=\nu= 4\cdot10^{-5}\) the solution ceases to exist. For larger values of \(\nu\), where the solution does exist, the line \(E_0=-\beta\) defines the border between Zone C and Zone B as predicted theoretically, see Section 5.3.

A similar diagram appears in Serrin [6]. It uses a slightly different (but equivalent) formulation of the problem, but it is derived through a completely different process. Our bifurcation diagram in Figure 11 is based on computation of the numerical solution of 100 as detailed in the previous paragraph. The diagram of Serrin is a consequence of a series of theoretical estimates and bounds that provide necessary and/or sufficient conditions on the underlying parameters, so that the solutions to the corresponding initial and boundary value problem belong to one of the desired zones; see [6]. Qualitatively the two diagrams are quite similar, in terms of the overall structure and shape of the borders between the various zones. However, quantitively there are some differences, due to the different sets of variables and constants used in the two approaches. In our case solutions can exist when the pressure constant \(E_0\) is negative, a region which is not addressed in Serrin’s work, where solutions are obtained only for positive values of pressure. The value \(P=1\) in [6] is the point where all zones converge and beyond this point along the axis \(\nu = 0\), solution ceases to exist. In the present work this value corresponds to \(E_0=-\frac{1}{2}\) and the whole diagram is mirrored along this line when compared to that in [6]. In both diagrams the two vertical lines define the border between zones B and C. However, it’s worth noticing that, apart from having different(opposite) orientations, the curved boundary that separates the NoSolution zone from the other zones is very similar in both works.

8 Navier-Stokes equations in Cylindrical Coordinates↩︎

Cylindrical coordinates \((r,\vartheta,z)\) are connected with rectangular coordinates through the transformation \(x_1 = r \,\cos\vartheta\), \(x_2 = r\,\sin\vartheta\), \(x_3 = z\), see Fig. [fig:cyl-coord-3d], while the associated orthonormal system attached to \(P\) is \[\vec{e}_r = (\cos\vartheta,\sin\vartheta,0), \,\; \vec{e}_\vartheta = (-\sin\vartheta,\cos\vartheta,0), \,\;\vec{e}_z = (0, 0, 1).\]

The velocity vector \(\vec{u}\) is expressed in cylindrical coordinates as \[\vec{u} = u(r,\vartheta,z,t) \vec{e}_r + v(r,\vartheta,z,t) \vec{e}_{\vartheta} + w(r,\vartheta,z,t) \vec{e}_z,\] while the vorticity \(\vec{\omega}\) takes the form \[\vec{\omega} = \bigg(\frac{1}{r}\frac{\partial u}{\partial \vartheta} -\frac{\partial v}{\partial z}\bigg) \vec{e}_r + \bigg( \frac{\partial u}{\partial z}-\frac{\partial w}{\partial r}\bigg) \vec{e}_\vartheta + \bigg( \frac{1}{r}\frac{\partial}{\partial r} (rv)- \frac{1}{r}\frac{\partial w}{\partial \vartheta} \bigg) \vec{e}_z.\] A change of variables for the Navier - Stokes equations 8 yields the following representation in cylindrical coordinates, see [25], \[\begin{align} \tag{160} \frac{\partial u}{\partial t} + (\vec{u} \cdot \nabla) u - \frac{v^2}{r} & = -\frac{\partial p}{\partial r} + \nu \Big[\Delta u - \frac{u}{r^2} -\frac{2}{r^2}\frac{\partial v}{\partial \vartheta} \Big], \\ \tag{161} \frac{\partial v}{\partial t} + (\vec{u} \cdot \nabla) v + \frac{uv}{r} & = -\frac{1}{r}\frac{\partial p}{\partial \vartheta}+\nu \Big[\Delta v - \frac{v}{r^2} + \frac{2}{r^2}\frac{\partial u}{\partial \vartheta} \Big], \\ \tag{162} \frac{\partial w}{\partial t} + (\vec{u} \cdot \nabla) w \qquad & = - \frac{\partial p}{\partial z} + \nu \Delta w, \\ \tag{163} \frac{1}{r} \frac{\partial }{\partial r} (ru) +\frac{1}{r}\frac{\partial v}{\partial \vartheta} + \frac{\partial w}{\partial z} & = 0. \end{align}\] where the operators \(\vec{u} \cdot \nabla\) and \(\Delta\) are \[(\vec{u} \cdot \nabla) = u \frac{\partial }{\partial r} + \frac{1}{r} v \frac{\partial }{\partial \vartheta} + w \frac{\partial }{\partial z} \qquad \text{and} \qquad \Delta = \frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial }{\partial r}\Big)+\frac{1}{r^2} \frac{\partial^2 }{\partial \vartheta^2} + \frac{\partial^2 }{\partial z^2}\]

9 Tornadoes↩︎

Tornadoes are considered among the most extreme and violent weather phenomena on Earth. They can occur under appropriate circumstances in all continents expect Antarctica, at various seasonal times and can be hazardous causing loss of human lives and extensive property damage. According to meteorologists, a tornado is defined as a rapidly rotating mass of air that extends downward from a cumuliform cloud, that is a cloud formed due to vertical motion of air parcels, to the ground. There exists several types of tornadoes, such as landspouts and waterspouts, but the majority of the most destructive tornadoes are known as supercells since their generation takes place within supercell thunderstorms.

Although there is a high interest in forecasting such hazardous tornadoes, it remains a challenging task for researchers to predict when a supercell thunderstorm will lead to a tornado. It is observed that not all supercells are tornadic since a combination of atmospheric instability (caused by the storm) with a wind shear, i.e. a variation of wind speed and direction with altitude, is required for tornado formation. These two ingredients are important for tornado formation, however the process by which tornadoes are formed is still not fully understood and thus, difficult to predict. For a detailed presentation on the subject of tornadoes and tornado formation, we refer to [26] and [27].

9.1 Modeling Tornadoes↩︎

Due to the complexity of tornadoes, the current knowledge about them comes mainly from laboratory experiments and numerical models of idealized supercell thunderstorms, [2]. In 1972, Ward [28] conducted a pioneering laboratory experiment reproducing a tornado-like flow considering a fluid with constant density. The idea was to create a flow using a fan that passes through a hole of radius \(r_0\) and is placed above a rotating plate in some distance \(h\), under the assumption that the ratio of \(h/r\) is small. Based on this work, several experimental and numerical simulations have taken place, referred to as Ward-type simulations, [2]. It was shown that the vortex form changes as the rotation increases, from single-celled (centerline updraft) to single-celled below to double-celled above, to double-celled (central downdraft surrounded by updraft) to multiple vortices, [2]. Also, this structural change is largely independent of the Reynolds number.

Fiedler [29] proposed an idealization of a tornado-like flow that is defined on a closed domain and is for theoretical analysis. Here, the buoyancy force is taken into consideration. The behavior of such flows can be analyzed numerically using the axisymmetric, incompressible Navier-Stokes equations in cylindrical coordinates. Hence, the model takes the form \[\begin{align} \frac{Du}{Dt} &= \frac{v^2}{r} + 2\Omega v + \frac{1}{Re} \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial u}{\partial r}\Big)+ \frac{\partial^2 u}{\partial z^2} - \frac{u}{r^2} \Big] - \frac{\partial p}{\partial r}, \\ \frac{Dv}{Dt} & = - \frac{uv}{r} - 2\Omega u + \frac{1}{Re} \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial v}{\partial r}\Big)+ \frac{\partial^2 v}{\partial z^2} - \frac{v}{r^2} \Big], \\ \frac{Dw}{Dt} & = b + \frac{1}{Re} \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial w}{\partial r}\Big)+ \frac{\partial^2 w}{\partial z^2} \Big] - \frac{\partial p}{\partial z}, \\ 0&=\frac{1}{r} \frac{\partial }{\partial r} (ru) + \frac{\partial w}{\partial z} , \end{align}\] where \(\frac{D}{Dt} = \frac{\partial }{\partial t} + u \frac{\partial }{\partial r}+ w\frac{\partial }{\partial z}\) stands for the material derivative, \(b\) is the buoyancy and \(\Omega\) is the non-dimensional swirl ratio which depends on both the angular momentum and the buoyancy. Numerical experiments of this model produce results analogous to Ward-type experiments for different values of \(\Omega\) and \(Re\), [2].

In addition, various analytical models have been introduced to describe a tornado-like flow behavior. Assuming that a vortex line resembles the tornado core, we consider again the incompressible axisymmetric Euler and Navier - Stokes equations.

Let us review some widely used vortex models. A detailed presentation can be found in [30] and [31] and in references therein.

The Rankine vortex model is considered as the simplest one. Here the flow is assumed to be one-dimensional, steady, inviscid and all body forces are omitted. Hence, the model takes the form \[\begin{align} \frac{d p(r)}{d r} =\rho \frac{v^2}{r}, \end{align}\] where \(\rho\) is the density. Also, it is assumed that the velocity component is discontinuous and is written as \[\begin{align} \bar{v}(\bar{r}) = \left\{ \begin{aligned} \bar{r} &\quad for \,\, \bar{r}<1,\\ \frac{1}{\bar{r}} &\quad \,\, for \,\, \bar{r}>1. \end{aligned} \right. \end{align}\] where \(\bar{v} = \frac{v}{v_{max}}\) is the normalized velocity and \(\bar{r}=\frac{r}{R}\) is the normalized distance for \(R\) the radius of the core vortex. Sometimes, a modified version of velocity is used that is \[\bar{v} = \frac{2\bar{r}}{(1+\bar{r}^2)}.\] If the discontinuous velocity is considered, solving the differential equation yields the normalized pressure \(\bar{p}(\bar{r}) = \frac{p(r)}{\rho v_{max}^2}\) that is \[\begin{align} \bar{p}(\bar{r}) = \left\{ \begin{aligned} \bar{p}(0) + \frac{1}{2} \bar{r}^2 &\quad for \,\, \bar{r}<1,\\ \bar{p}|_{r\to\infty} -\frac{1}{\bar{r}^2} &\quad \,\, for \,\, \bar{r}>1. \end{aligned} \right. \end{align}\]

Another vortex model is the Burgers-Rott, where the flow is assumed to be steady, with constant viscosity and zero body forces. Moreover, it is assumed that \(u=u(r)\), \(v=v(r)\), \(w=w(z)\) and \(p=p(r,z)\). The model then has the following form \[\label{s1} \begin{align} u \frac{\partial u}{\partial r} - \frac{v^2}{r} & = \mu \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial u}{\partial r}\Big) - \frac{u}{r^2} \Big] - \frac{1}{\rho}\frac{\partial p}{\partial r}, \\ u \frac{\partial v}{\partial r} + \frac{uv}{r} & = \mu \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial v}{\partial r}\Big) - \frac{v}{r^2} \Big], \\ w \frac{\partial w}{\partial z} \qquad & = - \frac{1}{\rho}\frac{\partial p}{\partial z}, \\ \frac{1}{r} \frac{\partial }{\partial r} (ru) + \frac{\partial w}{\partial z} & = 0, \end{align}\tag{164}\] where \(\rho\) is the density and \(\mu\) the dynamic viscosity. It is also assumed that \[\begin{align} \bar{w}(\bar{z}) &= 2\,\alpha \,\bar{z}, \\ \bar{u}(\bar{r}) &= -\alpha \,\bar{r}, \end{align}\] where \(\bar{z}=\frac{z}{R}\) is the normalized vertical height, \(\bar{u}=\frac{u}{v_{max}}\) and \(\bar{w}=\frac{w}{v_{max}}\) are the normalized velocities and \(\alpha = \frac{2\mu}{v_{max}}\). Under these assumptions, solving the system 164 implies the following \[\begin{align} \bar{v}(\bar{r}) &= \frac{1}{\bar{r}} (1-exp(-\bar{r}^2)), \quad \textrm{and} \\ \bar{p}(\bar{r}, \bar{z}) &= \bar{p}(0,0) + \int_0^{\bar{r}} \frac{\bar{v}^2(s)}{s} ds - \frac{\bar{\alpha}}{2} (\bar{r}^2 + 4 \bar{z}^2) \end{align}\]

The Sullivan vortex model has also been used widely to model tornado-like flows. As in the case of the Burgers-Rott model, we consider a flow that is stationary, with constant viscosity and zero body forces. In addition, it is considered that velocity components are given in the form \(u=u(r)\), \(v=v(r)\), \(w=w(r,z)\), while pressure is of the form \(p=p(r,z)\). One may conclude to the following \[\begin{align} \bar{u}(\bar{r}) &= -\bar{\alpha}\bar{r} +\frac{2b\bar{v}}{\bar{r}} (1-e^{-\bar{r}^2}),\\ \bar{v}(\bar{r}) &= \frac{1}{\bar{r}} \frac{H(x)}{H(\infty)}, \\ \bar{w}(\bar{r},\bar{z}) &= 2 \bar{\alpha}\bar{z} (1-b\,e^{-\bar{r}^2}), \end{align}\] where \(H(x) = \int_0^{x} e^{-s + 3 \int_0^s \frac{1}{\sigma} (1-e^{-\sigma^2}) d\sigma} ds\) for \(x=\bar{r}^2\). It is worth mentioning that although the Sullivan and the Burgers-Rott models have some similarities, the Sullivan model allows the generation of a double-celled vortex while Burgers-Rott model does not.

9.2 Mathematical Approach↩︎

A theoretically sound approach towards study of tornadoes was introduced by Long [3], [4]. Assuming the tornado core is modeled by a semi-infinite vortex line in a fluid interacting with a plane boundary surface, he presented the reduction of incompressible Navier-Stokes equations to a system of differential equations motivated by boundary layer theory. Several subsequent studies [14], [20], [32], [33] took a similar direction and studied the formation of a boundary layer considering the near-axis boundary layer approximation to the incompressible axisymmetric Navier-Stokes equations. This is usually referred as quasi-cylindrical approximation and leads to the system \[\begin{align*} \frac{v^2}{r} & = \frac{\partial p}{\partial r}, \\ u \frac{\partial v}{\partial r} + w \frac{\partial v}{\partial z} + \frac{uv}{r} & = \nu \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial v}{\partial r}\Big) - \frac{v}{r^2} \Big], \\ u \frac{\partial w}{\partial r} + w \frac{\partial w}{\partial z} \qquad & = \nu \Big[\frac{1}{r}\frac{\partial }{\partial r} \Big(r \frac{\partial w}{\partial r}\Big) \qquad \Big] - \frac{\partial p}{\partial z}, \\ \frac{1}{r} \frac{\partial }{\partial r} (ru) + \frac{\partial w}{\partial z} & = 0, \end{align*}\] The boundary layer is associated in the literature with the vortex breakdown, that is the change of direction of the flow near the boundary as \(\nu\to0\).

Independently, Goldshtik (1960) showed that a similar reduction of incompressible axisymmetric Navier-Stokes equations to a system of differential equations leads to a ‘paradoxical’ exact solution that vanishes for some values of Reynolds number, [5]. Serrin (1972) broadened this class of solutions and described the existence of three different solution profiles depending on an arbitrary parameter and the kinematic viscosity, [6]. Following this work, several authors have extended the study to the generalized case of conical flows, [20], [13], [14]. The ideas of Long and Goldshtik have been applied to investigate the formation of a boundary layer and the loss of existence of such solutions using different boundary conditions or a modified self-similar ansatz, [9], [10], [32]–[34]. This line of research is systematized in the present manuscript by studying stationary solutions of the axisymmetric Navier-Stokes equations.

The above References concern the interaction of a vortex with a boundary when the flow is assumed stationary. Models have been devised concerning the motion of the vortex core structure and its coupling with environmental flows, and the reader is referred to [35], [36] and references therein regarding that subject.

Acknowledgment. Research partially supported by King Abdullah University of Science and Technology (KAUST) baseline funds.

Conflicts of interest. The authors have no conflicts of interest to report.

References↩︎

---

1. Department of Mathematics and Applied Mathematics, University of Crete, Heraklion, 70013 Crete, Greece. Email: thodoros.katsaounis@uoc.gr↩︎

2. Institute of Applied and Computational Mathematics, FORTH, Heraklion, Greece.↩︎

3. Computer, Electrical, Mathematical Sciences & Engineering Division, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia. Email: ioanna.mousikou@kaust.edu.sa↩︎

4. Computer, Electrical, Mathematical Sciences & Engineering Division, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia. Email: athanasios.tzavaras@kaust.edu.sa↩︎