A fluid mechanic's analysis of the teacup singularity

TL;DR

This paper investigates the formation of singularities in inviscid, wall-bounded fluid flows using numerical simulations of the Euler equations, emphasizing the role of pressure and flow confinement in the blowup process.

Contribution

It introduces a pressure decomposition method and a primitive-variables model to analyze the singularity mechanism in axisymmetric swirling flows.

Findings

Numerical simulations support the existence of finite-time singularities.

Pressure confinement within the cylinder is crucial for velocity gradient blowup.

The proposed model captures key mechanics of the singularity formation.

Abstract

The mechanism for singularity formation in an inviscid wall-bounded fluid flow is investigated. The incompressible Euler equations are numerically simulated in a cylindrical container. The flow is axisymmetric with swirl. The simulations reproduce and corroborate aspects of prior studies reporting strong evidence for a finite-time singularity. The analysis here focuses on the interplay between inertia and pressure, rather than on vorticity. Linearity of the pressure Poisson equation is exploited to decompose the pressure field into independent contributions arising from the meridional flow and from the swirl, and enforcing incompressibility and enforcing flow confinement. The key pressure field driving the blowup of velocity gradients is that confining the fluid within the cylinder walls. A model is presented based on a primitive-variables formulation of the Euler equations on the cylinder wall, with closure coming from how pressure is determined from velocity. The model captures key features in the mechanics of the blowup scenario.

Tables

Table 1. Table 1: Tabulated principal pressure curvatures and other quantities at t = 0.0031 𝑡 0.0031 t=0.0031 .

Curvature	value		Curvature	value		Quantity	value
$P$	$-1.9877\times {10}^7$		$Q$	$-5.1872\times {10}^7$		${{\partial }_r{p}_a|}_c$	$-8.6991\times {10}^3$
$P_{2D}$	$-2.8556\times {10}^7$		$Q_{2D}$	$-4.3192\times {10}^7$		${{\partial }_r{p}_c|}_c$	$8.6991\times {10}^3$
$P_{\mathrm{swirl}}$	$8.6797\times {10}^6$		$Q_{\mathrm{swirl}}$	$-8.6797\times {10}^6$		$V$	$5.9895\times {10}^3$
$P_a$	$9.5265\times {10}^6$		$Q_a$	$-9.5178\times {10}^6$		$W$	$-5.9895\times {10}^3$
$P_b$	$-8.4676\times {10}^5$		$Q_b$	$5.7998\times {10}^5$		$\Omega$	$1.5428\times {10}^5$
$P_c$	$0$		$Q_c$	$2.5808\times {10}^5$

Equations

$${\partial }_{t}u+u\cdot \nabla u$$

$$\nabla \cdot u$$

$$u(r,z,t)=u_r(r,z,t){\hat{e}}_r+u_{\theta }(r,z,t){\hat{e}}_{\theta }+u_z(r,z,t){\hat{e}}_z,$$

$${u_r|}_{r=1}=0.$$

$$u(r,z,t=0)=100r{e}^{-30(1-r^2{)}^4}\sin \left(\frac{2\pi }{L}z\right){\hat{e}}_{\theta }.$$

$$\left(u_r,u_{\theta },u_z\right)(r,z,t)=\left(u_r,-u_{\theta },-u_z\right)(r,-z,t).$$

$${\partial }_t\left(\nabla \cdot u\right)+\nabla \cdot \left(u\cdot \nabla u\right)=-{\nabla }^{2}p.$$

$${\nabla }^{2}p=-\nabla \cdot \left(u\cdot \nabla u\right).$$

$${\nabla }^{2}p=-\nabla \cdot \left(u\cdot \nabla u\right)\equiv S,{{\partial }_{r}p|}_{r=1}={u_{\theta }^2|}_{r=1}\equiv b,$$

$$W\equiv {{\partial }_z{u}_z|}_c,\Omega \equiv {{\partial }_z{u}_{\theta }|}_c,V\equiv {{\partial }_r{u}_r|}_c,$$

$$P\equiv {{\partial }_{zz}p|}_c,Q\equiv {{\partial }_{rr}p|}_c,$$

$$\dot{W}+W^2=-P,\dot{\Omega }+W\Omega =0,\dot{V}+V^2=-Q.$$

$$\dot{W}+W^2=-P,$$

$$\dot{\Omega }+W\Omega =0,$$

$$P+Q=-2W^2$$

$$a^2\equiv Q/P,$$

$$\dot{W}=-\frac{W^2}{\gamma },\dot{\Omega }=-W\Omega$$

$$\gamma \equiv \frac{a^2+1}{a^2-1}=\frac{Q+P}{Q-P}$$

$$W(t)=-\frac{\gamma }{T-t}\sim (T-t{)}^{-1},\Omega (t)=\frac{{\Omega }_0{T}^{\gamma }}{(T-t{)}^{\gamma }}\sim (T-t{)}^{-\gamma },$$

$${\partial }_t{u}_z+u_z{\partial }_z{u}_z=0,{\partial }_t{u}_{\theta }+u_z{\partial }_z{u}_{\theta }=0.$$

$$\dot{W}+W^2=0,\dot{\Omega }+W\Omega =0.$$

$$W(t)=-\frac{1}{T-t},\Omega (t)=\frac{{\Omega }_{0}T}{T-t},$$

$$\dot{W}+W^2=W^2,\dot{\Omega }+W\Omega =0.$$

$$W(t)=W_0,\Omega (t)={\Omega }_0\exp (-W_{0}t).$$

$${\nabla }^2{p}_{2D}$$

$${\nabla }^2{p}_{\mathrm{swirl}}$$

$$P=P_{2D}+P_{\mathrm{swirl}},Q=Q_{2D}+Q_{\mathrm{swirl}},$$

$$P_{2D}+Q_{2D}=-2W^2,P_{\mathrm{swirl}}+Q_{\mathrm{swirl}}=0.$$

$$|P|=|P_{2D}+P_{\mathrm{swirl}}|<|Q_{2D}+Q_{\mathrm{swirl}}|=|Q|.$$

$${\nabla }^2{p}_a=0,{{\partial }_r{p}_a|}_{r=1}$$

$${\nabla }^2{p}_b={\tilde{S}}_{\mathrm{swirl}},{{\partial }_r{p}_b|}_{r=1}$$

$${\nabla }^2{p}_c=⟨S_{\mathrm{swirl}}⟩,{{\partial }_r{p}_c|}_{r=1}$$

$${\partial }_t\omega +u{\partial }_z\omega ={\partial }_z\theta ,{\partial }_t\theta +u{\partial }_z\theta =0,$$

$${\partial }_{z}u=H(\omega ).$$

$$H(\theta )=H(b)=H(⟨b⟩+\tilde{b})=H(\tilde{b})=H({{\partial }_r{p}_a|}_{r=1})=-{{\partial }_z{p}_a|}_{r=1}.$$

Full text

Euler equations, finite-time singularity, pressure, swirl, primitive variables, wall bounded

I Introduction

In 1926 Einstein published a short paper explaining the meandering of rivers Einstein . He famously began the paper by discussing the secondary flow generated in a stirred teacup – the flow now widely known to be responsible for the collection of tea leaves at the centre of a stirred cup of tea. In 2014, Luo and Hou presented detailed numerical evidence of a finite-time singularity at the boundary of a rotating, incompressible, inviscid flow LH_PNAS ; LH_MMS . The key to generating this singularity is the teacup effect. The present work is not aimed at proving the existence of a singularity for this flow, nor is it aimed at generating more highly resolved numerical evidence for the singularity than already exists. Rather, I assume that the flow simulated by Luo and Hou genuinely develops a singularity in finite time. My goal is to understand, from a fluid-mechanics perspective, why.

II Preliminaries

II.1 Problem statement

The flow under investigation is depicted in Fig. 1. The system is initialised with a pure azimuthal flow (swirl) having a sinusoidal dependence on the axial coordinate $z$. A pressure field is instantaneously generated to provide the radially inward force necessary to keep fluid parcels moving along circular paths. This results in high pressure at the cylinder wall where the circulation is largest ($z=\pm L/4$) and low pressure where there is no azimuthal flow ($z=0$ and $z=\pm L/2$). Necessarily, then, there is a vertical variation in the pressure at the cylinder wall and this drives a secondary meridional flow. This is the teacup effect – the portion of the fluid just from $z=0$ to $z=L/4$ corresponds to a cup of tea. (In an actual cup of tea the variation in swirl with $z$ is due to a boundary layer at the bottom of the cup. Here we disregard viscous effects even though they play a role in the flow of real tea in a teacup.)

We consider inviscid fluid flow governed by the incompressible Euler equations

[TABLE]

The initial condition employed by Luo and Hou, and reproduced here, is a pure swirl

[TABLE]

This initial condition possesses symmetries that are preserved under evolution of (1). The most important is centro symmetry about $z=0$

[TABLE]

The full set of symmetry planes is $z_{j}=jL/4$, $j=0,\pm 1,\pm 2$; $u_{r}$ is even and $u_{z}$ is odd about these all planes; $u_{\theta}$ is odd about planes $z_{o}$, $z_{\pm 2}$ and is even about planes $z_{\pm 1}$. The pressure $p$ is even about all four planes.

Extensive analysis of finely resolved numerical simulations of the Euler equations indicates that starting from the above initial condition, the flow evolves to form a singularity on the critical ring, $(r=1,z=0)$, at time $T\simeq 0.0035056$  LH_PNAS ; LH_MMS ; LH_Review . In the present work, simulations are well resolved to time $t=0.0031$. Details of the simulations are given in Appendix D. I rely heavily on the studies of Luo and Hou (hereafter referred to as LH), to know that the flow at $t=0.0031$ is indicative of the flow all the way to $t=0.003505$, extremely close to the singularity time. To be clear, the simulations presented here are not aimed at numerically establishing a singularity (LH have already done this), but instead at understanding the physical mechanisms at work, and for this purpose they are adequate.

II.2 Mechanics

Pressure is the only stress acting within an inviscid fluid and it is the only means to provide force to, and thereby accelerate, the flow. It is at the heart of the teacup effect and it is therefore natural to investigate its role in the singularity. In general, stress is a tensor field $\tau$ whose divergence gives the net force acting on infinitesimal fluid parcels. Pressure is isotropic and so for inviscid flow $\tau$ has only diagonal components: $\tau_{ij}=-P\delta_{ij}$, where $P$ is the pressure field, and $\delta_{ij}$ is the Kronecker delta. Hence $\nabla\cdot\tau=-\nabla P$ is the force per volume acting within the fluid. For incompressible flow, the fluid’s mass density $\rho$ is constant and we define $p\equiv P/\rho$, so that $-\nabla p$ sets the fluid acceleration and thus appears on the right-hand-side of the momentum equation (1a). Subsequently the symbol $P$ will be used for another quantity and we will refer to $p$ simply as pressure. Pressure gradients with $-\nabla p$ anti-parallel to velocity $u$ are known as adverse pressure gradients and result in flow deceleration (decrease in fluid speed).

The role of pressure in incompressible flow is seen by taking the divergence of (1a)

[TABLE]

This equation governs the evolution of the flow divergence. Given a divergence-free velocity field $u$ satisfying (1b), in general $\nabla\cdot\left(u\cdot\nabla u\right)$ will not be zero, meaning that nonlinearity acting alone does not maintain incompressibility. A pressure field is generated within the fluid (simultaneously everywhere) to accelerate the flow exactly so as to counterbalance this effect of nonlinearity. From (3), the relationship between pressure and velocity required to maintain a divergence-free flow is the Poisson equation

[TABLE]

This is not the full story, however. The flow of interest is wall bounded and this puts a condition on the stress field within the fluid. The initial velocity field satisfies (1c) and thus has no radial component at the cylinder wall. From the $\hat{e}_{r}$ component of the momentum equation at the wall, this will be maintained as long as $\left.{\partial_{r}p}\right|_{r=1}=\left.{u_{\theta}^{2}}\right|_{r=1}$. Thus, pressure is determined by a Poisson equation together with its boundary condition

[TABLE]

where these expressions define the source term ${S}$ and the boundary term ${b}$. As long as $p$ satisfies (4), the flow evolving under (1a) will remain incompressible and confined within the cylinder. An important focus of this work will be disentangling the contributions to the stress associated with incompressibility from those associated with flow confinement.

III Basics of the singularity mechanism

Figure 2 presents a quantitative overview of the flow dynamics with visualisation of the initial and final states from the numerical simulation. We see the teacup effect: the initial primary flow (pure swirl) results in high pressure on the cylinder wall at $z=\pm L/4$, where the swirl is largest. This in turn results in a vertical component of the pressure gradient that produces a secondary flow driving fluid on the cylinder wall toward the midplane $z=0$, (and by symmetry also toward $z=\pm L/2$). Swirl of opposite signs from above and below the midplane is thus transported towards $z=0$, resulting in intense radial vorticity $\omega_{r}=-\partial_{z}u_{\theta}$ on the critical ring. This has been described clearly by LH LH_PNAS ; LH_MMS ; LH_Review , who show very strong evidence that the flow continues to develop a singularity in a nearly, (but not exactly  chae_Tsai_2015 ; Sperone_2017 ), self-similar way. The final state of the present simulations shown in Fig 2 is representative of the flow as it approaches the singularity.

Figure 3(a) shows an enlargement of the final flow in a very small region around the critical ring. Analysis of this state will be a major focus of the paper. The critical ring is a saddle point of the meridional flow. A local pressure maximum exists on the critical ring (barely visible in Fig. 2). This secondary local maximum in the pressure field first appears at time $t\simeq 0.002$ and accounts for the stress required to accelerate the flow around the saddle point – bending incoming axial velocity near the cylinder wall to radially inward velocity near the midplane. This pressure field is similar to that reported by LH at $t=0.003505$, very close to the singularity time $T\simeq 0.0035056$. (See Ref. LH_MMS , but note that its Fig. 17 has a distorted aspect ratio.) LH emphasise that the pressure maximum on the critical ring means that there is locally an adverse axial pressure gradient decelerating the incoming axial flow on the cylinder wall. This is an important point. However, it does not mean that the pressure maximum inhibits the singularity. On the contrary, a pressure maximum like that in Fig. 3(a) can drive a singularity. This fact is central to this work.

To understand the mechanics of this particular situation, we turn to the velocity-gradient dynamics on the critical ring. Differentiating velocity gives the velocity-gradient tensor $\nabla u$ and differentiating the pressure gradient gives the pressure Hessian $\nabla(\nabla p)$. Symmetries dictate that on the critical ring the only non-zero derivatives entering these are

[TABLE]

where $\left.{}\right|_{c}$ means evaluated on the critical ring. $P$ and $Q$ will be especially important in what follows. I refer to these as pressure curvatures since they are the principal curvatures of a graph of the pressure $p(r,z)$. It is to be understood that I always mean “on the critical ring” when referring to these curvatures. Straightforward differentiation of (1a) gives

[TABLE]

By incompressibility on the critical ring: $V+W=0$. Thus, $V$ can be eliminated, giving the velocity-gradient dynamics

[TABLE]

The meaning associated with each equation is indicated. The equations are exact, and while they are not closed ((6c) is insufficient to determine $P$ and $Q$ separately), they are extremely useful in examining what transpires in singularity formation. For this flow, $-\left.{\omega_{r}}\right|_{c}=\left.{\partial_{z}u_{\theta}}\right|_{c}=\Omega$ is the absolute vorticity maximum LH_PNAS ; LH_MMS , so $\Omega=\|\omega\|_{\infty}=\left.{\partial_{z}u_{\theta}}\right|_{c}$ as indicated in Fig. 2. Equation (6c) is the pressure Poisson equation evaluated on the critical ring, but equally it is the geometrical statement that sum of principal curvatures is twice the mean curvature, where the mean curvature of $p$ is $-W^{2}$. (See works by D. Chae and collaborators chae2008incompressible ; chae2008blow ; constantin2008singular ; chae2010lagrangian ; Chae_etal_2012 for more general treatments of the Euler equations in the velocity-gradient formulation, including several blowup scenarios.)

From Fig. 3(a) we see that the principal pressure curvatures, $P$ and $Q$, are both negative (a pressure maximum occurs on the critical ring), but that they are not equal. The axial curvature is smaller in magnitude than the radial curvature, that is $|P|<|Q|$. This can be seen in the ratio of axial to radial length scales in the pressure contours. This leads us to define $a\geq 0$ by

[TABLE]

so that $a$ is this ratio of length scales. The case $a>1$ corresponds to $|P|<|Q|$ and is seen in Fig. 3(a).

To understand the importance of $|P|<|Q|$ to blowup, we proceed as follows. Using (7) to eliminate $Q$ from (6c) gives $P=-2W^{2}/(a^{2}+1)$, which can then be used to eliminate $P$ from (6a). The velocity-gradient equations (6) then become

[TABLE]

where

[TABLE]

The case of interest $\infty>a>1$ corresponds to $1<\gamma<\infty$.

The solution to Eqs. (8) with $\gamma$ constant is simple and captures the essence of the blowup described by these equations. Effectively we set $\gamma$ to its limiting value, assumed to be finite, as $t\to T$. We are interested in a saddle point flow in the meridional plane, with fluid converging towards the critical ring in the axial direction. We are only ever interested in this situation and always assume $W(t=0)=W_{0}<0$. The solution to Eqs. (8) is then

[TABLE]

where $T=-\gamma/W_{0}>0$ is the singularity time and $\Omega_{0}=\Omega(0)$. These are the known divergences as $t\to T$ LH_PNAS ; LH_MMS . In particular, the vorticity $\Omega=\|\omega\|_{\infty}$ diverges with exponent $-\gamma$. All other divergences associated with the singularity follow from invariances of the Euler equations and the value of $\gamma$. By treating $\gamma$ as a constant, we obtain the scaling of the blowup as a simple exact solution to (8). This will be useful in what follows.

We know from LH that the vorticity diverges with exponent $\gamma\simeq 2.46$, corresponding to $a\simeq 1.54$. The corresponding ratio of length scales is indicated in Fig. 3(a). The contours do not exactly manifest this ratio of scales, in part because contours are a finite distance from the critical ring and in part because the flow is seen at a time a finite distance from the singularity time. From the data at $t=0.0031$, $\sqrt{Q/P}\simeq 1.62$. (See data in Table 1.)

The fundamental point is the following. Incompressibility locks axial contraction and radial expansion together such that it is not the signs of principal pressure curvatures $P$ and $Q$ that are important for singularity formation; it is their inequality. A persistent inequality in the pressure curvatures on the critical ring can drive the flow to a singularity. Of interest here is flow converging axially toward the critical ring with $|P|<|Q|$, so the axial curvature is smaller than the radial curvature in magnitude. The pressure contours in Fig. 3(a) are the signature of this simple mechanism. If the flow evolves such that this situation persists, (that is such that $\inf_{t\geq 0}Q/P>1$), then the solution will blow up. One can deduce from the results of LH that a ratio of pressure curvatures of approximately the same amount as is seen in Fig. 3(a) exists as close to the singularity time as they could resolve (Fig. 17 of Ref. LH_MMS ).

IV Illustrative cases

Before continuing to a detailed analysis of the pressure field, I consider the the velocity-gradient dynamics (6) in two limiting cases. These cases will appear again later in the paper (see Fig. 7), and they are very useful for understanding the interplay between inertia and pressure in the mechanics of the singularity.

Consider simply dropping pressure and the incompressibility constraint from the Euler equations, and also dropping the radial component for the momentum equation. On the cylinder wall, the remaining two components of the momentum equation become Burger’s equation and advection of swirl as a passive scalar:

[TABLE]

The velocity-gradient dynamics on the critical ring become

[TABLE]

(Pressure does not appear and (6)(c) is dropped.) Starting with a flow converging towards $z=0$, $W_{0}<0$, these equation have blowup given by

[TABLE]

where $T=-1/W_{0}>0$ is the singularity time. This is just a special case of Eq. (10) with $\gamma=1$. In the absence of stresses, there is no deceleration of the fluid parcels converging towards $z=0$, resulting in the well-known blowup of Burger’s equation. This illustrates how inertia, or equivalently the associated advective nonlinearity in Eulerian coordinates, itself can easily lead to a finite-time singularity.

Consider now the case in which principal pressure curvatures are equal at all times: $P=Q$. This would correspond to pressure contours locally forming circular arcs about the critical ring in the meridional plane (similar to what is seen in Fig. 3(b), although those contours are not perfectly circular). With this assumption, Eqs. (6) are closed because both pressure curvatures equal the mean curvature: $P=Q=-W^{2}$. With this, the velocity-gradient dynamics on the critical ring become

[TABLE]

(This case corresponds to $\gamma=\infty$.) The system does not develop a singularity and instead has solution

[TABLE]

Since $W_{0}<0$, the vorticity grows exponentially, but only exponentially in time. This illustrates what is observed to be the normal situation for incompressible inviscid flow – the stress that develops within the flow to maintain incompressibility is such as to accelerate the fluid sufficiently to prevent blowup that would come from inertia acting alone. Algebraically, the term $W^{2}$ from inertia on the left-hand side of (13) is exactly balanced by the term $W^{2}$ from pressure on the right-hand side.

The case of interest, where $|P|<|Q|$ and hence $1<\gamma<\infty$, falls between the two extremes just considered. We have a flow configuration evolving under the full incompressible Euler equations, with the pressure stress acting, but such that the axial pressure curvature $P$ is too small to compensate inertia. As a result, fluid parcels converging towards $z=0$ are not sufficiently decelerated and a singularity occurs.

V Analysis of pressure

In this section I will analyse in depth the structure of the pressure field near the critical ring and show how it is dictated by specific aspects of the fluid flow. I will then use this information in Sec. VI to gain further insights into the blowup scenario.

V.1 Meridional and swirl pressure fields

We exploit the linearity of the Poisson equation (4) to separate pressure into contributions from distinct effects. To begin, the source term for the equation can be decomposed as ${S}={{S}}_{\rm 2D}+{{S}}_{\rm swirl}$, where ${{S}}_{\rm 2D}$ depends only on the meridional (2D) velocity components $(u_{r},u_{z})$ and ${{S}}_{\rm swirl}$ depends only on the swirl velocity $u_{\theta}$. (See Appendix A for details.) The boundary term ${b}$ in (4) also depends only on $u_{\theta}$. Thus, the pressure $p$ can be written as a linear superposition $p=p_{\rm 2D}+p_{\rm swirl}$, where

[TABLE]

These pressure fields are plotted in Fig. 3(b). Contours of $p_{\rm 2D}$ are nearly circular arcs indicating approximate rotational symmetry locally about the critical ring within the meridional plane. Contours of $p_{\rm swirl}$ are those of a hyperbolic saddle with the expected high pressure along the cylinder wall where the swirl is largest. Since we will be especially concerned with the axial momentum balance, these fields are plotted along the cylinder wall in Fig. 4(a).

Let

[TABLE]

where $P_{\rm 2D}=\left.{\partial^{2}_{z}p_{\rm 2D}}\right|_{c}$, $Q_{\rm 2D}=\left.{\partial^{2}_{r}p_{\rm 2D}}\right|_{c}$, etc, are the principal curvatures of the component pressure fields. (See Table 1.) Since $\left.{{{S}}_{\rm swirl}}\right|_{c}=0$, we have from (15b) and (6c)

[TABLE]

Hence the mean curvature of the pressure field $p$ is contained entirely in the component field $p_{\rm 2D}$. (This is obvious since both $p_{\rm 2D}$ and the mean curvature $-W^{2}$ are functions only of the meridional flow, and $p_{\rm swirl}$ is not.) Necessarily then the swirl pressure always has zero mean curvature.

The core cause for the inequality in the pressure curvatures, $|P|<|Q|$, is immediately evident. The near symmetry of the meridional pressure maximum implies that $P_{\rm 2D}\simeq Q_{\rm 2D}\simeq-W^{2}<0$, while for the saddle swirl pressure $P_{\rm swirl}>0>Q_{\rm swirl}$. Hence

[TABLE]

Stated simply – the pressure maximum from the meridional flow is flattened by the swirl pressure the axial direction, but it is steepened by the swirl pressure in the radial direction. This is seen in the visualisations of Fig. 3 and shown quantitatively along the cylinder wall in Fig. 4(a). To exploit fully this insight, more detailed information is required on the meridional and swirl pressure fields.

V.2 Meridional pressure

The pressure field $p_{\rm 2D}$ exists within the fluid to counter divergences that would be otherwise generated just by meridional flow. As the radial gradient of this field is zero at the cylinder wall, (15a), it does not contribute to fluid confinement. It is determined only by the instantaneous state of the meridional flow. In a region around the critical ring the meridional velocity field is a saddle that is approximately anti-symmetric under interchange of the axial and radial directions. See the streamlines in Fig. 3(a) where the Stokes streamfunction locally satisfies $\psi(r,z)\simeq\psi(1-z,1-r)$ for $z\geq 0$, and similarly for $z\leq 0$. Although the flow cannot globally respect such a symmetry, very near to the critical ring it does, approximately. Such a saddle flow is to be anticipated Kiselev_Sverak_2014 ; CKY_2015 ; Kiselev_Tan_2018 and the associated approximate rotational symmetry of $p_{\rm 2D}$ near a local maximum is not particularly surprising. The source term ${{S}}_{\rm 2D}$ is quadratic in velocity so an approximately anti-symmetric streamfunction implies an approximately symmetric pressure field.

However, it is the principal pressure curvatures on the critical ring that matter for singularity formation. So while the near symmetry of $p_{\rm 2D}$ seen in Fig. 3(b) suggests that $P_{\rm 2D}\simeq Q_{\rm 2D}$, it is necessary to examine these curvatures quantitatively, in particular to understand in what way they are not exactly equal. Figure 4(b) shows second derivatives of $p_{\rm 2D}$ along slices at the midplane, $z=0$, and at the cylinder wall, $r=1$. The general agreement between the two curves is a manifestation of the near symmetry of $p_{\rm 2D}$. However, the curves behave differently approaching the critical ring. Necessarily $\partial^{2}_{z}p_{\rm 2D}$ is even about $z=0$, since $p_{\rm 2D}$ is. There is no such constraint on $\partial^{2}_{r}p_{\rm 2D}$ at $r=1$. The important observation is that $Q_{\rm 2D}<P_{\rm 2D}<0$, and hence that $|P_{\rm 2D}|<|Q_{\rm 2D}|$. This means that the pressure curvatures associated with just the meridional flow are unequal with the ordering that promotes, rather than acts against, singularity formation. While this ordering does not seem a priori obvious, it appears from Fig. 4(b) to be a natural consequence of the conditions at the wall and symmetry plane. We will return to the importance of this ordering in Sec. VIVI.3.

V.3 Swirl pressure

The swirl pressure $p_{\rm swirl}$ not only maintains incompressibility of the flow, it also confines the fluid within the cylinder wall. Fully decoupling these two effects is not achievable for flow in a cylinder, but we can mostly separate them via the decomposition $p_{\rm swirl}=p_{a}+p_{b}+p_{c}$, where

[TABLE]

where $\langle\cdot\rangle$ denotes axial mean and tilde denotes axial fluctuations. These pressure components are plotted in Figs. 5 and 6(a). We also decompose the pressure curvatures, $P_{\rm swirl}=P_{a}+P_{b}+P_{c}$ and $Q_{\rm swirl}=Q_{a}+Q_{b}+Q_{c}$, with the obvious meanings. (See Appendix A for details of this decomposition as well as relationships that hold for the component curvatures.)

The most significant fact from the decomposition is best seen in Fig. 6(a). Near the critical ring, the axial variation of $p_{\rm swirl}$ is given almost exclusively by the component $p_{a}$. The positive axial curvature of $p_{\rm swirl}$ comes about from the $p_{a}$ component: $P_{\rm swirl}\simeq P_{a}>0$. (The radial curvatures satisfy $Q_{\rm swirl}\simeq Q_{a}<0$; see Table 1. We return to this shortly.) The stress field associated with the pressure $p_{a}$ exists throughout the fluid solely to provide force at the wall necessary to confine the flow within the cylinder – the accelerations it generates within the fluid have no effect on the divergence of the flow field. The pressure field $p_{a}$ is the essence of the teacup effect near the critical ring – axial variation of the swirl at the cylinder wall necessitates a pressure whose radial gradient at the wall confines the fluid and whose axial gradient then forces axial flow toward the critical ring.

The pressure field $p_{b}$ exists within the fluid to accelerate the flow and suppress divergences that would otherwise arise from spatial variations of $u_{\theta}$. This component is very weak near the critical ring: the range of values for $p_{b}$ is small in Fig. 5 and the curve corresponding to $p_{b}$ in Fig. 6(a) is nearly flat. Its curvatures have signs $P_{b}<0<Q_{b}$ meaning that it acts against singularity formation. However, these curvatures are an order of magnitude smaller than those of $p_{a}$, so the effect is very weak. (See Table 1.) Further from the critical ring, $p_{b}$ makes a more substantial contribution to the momentum balance, but this is not important to singularity formation.

The pressure component $p_{c}$ is easy to interpret physically. It is the axially-independent pressure field that would be generated in the pure swirl flow $\sqrt{\langle u_{\theta}^{2}\rangle}(r)\hat{e}_{\theta}$, whose speed at each $r$ is the axial r.m.s. of $u_{\theta}$. The radially-inward force $-\nabla p_{c}(r)$ is such as to curve each circular streamline of this r.m.s. swirl flow, both maintaining incompressibility and confining the fluid at the cylinder wall. While this pressure component is significant in the radial momentum balance, it contributes minimally, if at all, to the singularity. By definition $p_{c}$ does not vary with $z$, so it does not enter the axial momentum balance and $P_{c}=0$. While $Q_{c}$ is not zero, it is the smallest of all component pressure curvatures (see Table 1). (I suspect that $Q_{c}$ does not diverge at the singularity and hence plays no role in the blowup. See Appendix A.)

The essential aspect is this: the pressure component $p_{a}$ is the mechanism by which the swirl on the wall couples to the pressure field. It is the sole pressure component driving singularity formation. It is visually evident in Fig. 6(a) that $p_{a}$ is responsible for the positive curvature of $p_{\rm swirl}$ along the axial direction, and hence the favourable axial pressure gradient. It is less evident comparing Fig. 5(a) to Fig. 3(c) that $p_{a}$ dictates the negative radial curvature of $p_{\rm swirl}$. This is because $\left.{\nabla p_{\rm swirl}}\right|_{c}=0$, while $\left.{\nabla p_{a}}\right|_{c}\neq 0$. To account for this, in Fig. 5(d) we show $p_{a}+A(r-1)$, where $A=\left.{\partial_{r}p_{c}}\right|_{c}=-\left.{\partial_{r}p_{a}}\right|_{c}$. The field $A(r-1)$ is a linear approximation to $p_{c}$ at $r=1$, and from the contours in Fig. 5(c), $p_{c}$ is nearly linear in the region shown. Comparing Fig. 5(d) to Fig. 3(c) it becomes clear that near the critical ring $p_{\rm swirl}\simeq p_{a}+A(r-1)$. Since the term $A(r-1)$ is linear, the curvatures of $p_{a}+A(r-1)$ are dictated solely by those of $p_{a}$.

(Briefly, the issue here is directly related to the impossibility of separating the interior problem from the boundary condition for the axial mean in a cylinder (18c). This is only problematic in that $p_{a}$ does not “look like” $p_{\rm swirl}$ because $p_{a}$ lacks the component of the pressure field responsible for confining the axial mean flow and hence its gradient is non-zero on the critical ring. Conceptually though, even if one could separate the interior from the boundary for the axial mean, it would not necessarily be desirable to add the resulting field to $p_{a}$, other than for visual comparisons, because $p_{a}$ is the fundamental field that alone is responsible for the opposite-signed curvatures driving singularity formation. The axial mean is not important. This decomposition of the swirl pressure is a rich problem that I will not discuss further other than to note that in the Boussinesq system (Appendix B and Sec. VI) confinement and incompressibility can be fully separated.)

VI One-dimensional model and closure

I will now use facts learned from the pressure decomposition to gain insight into mechanics of the blowup scenario. The preceding analysis describes in detail the situation at one time instant, but does not address the persistence of this mechanism as the flow evolves. To do this I will examine a model based on a primitive-variables formulation of the Euler equations on the cylinder wall, with closure coming from our knowledge of how pressure is determined from velocity.

VI.1 Background

There is a rich literature on one-dimensional modelling of singularities in inviscid flow. See Choi_etal_2017 for a recent summary. For the cylinder flow, LH propose the model LH_PNAS ; LH_MMS

[TABLE]

with the identifications $\omega(z)\sim\left.{\omega_{\theta}}\right|_{r=1}$, $\theta(z)\sim\left.{u_{\theta}^{2}}\right|_{r=1}$, and $u(z)\sim\left.{u_{z}}\right|_{r=1}$. (We abuse notation, by conflicting with usage elsewhere in the paper and by not strictly distinguishing between model quantities and their full-flow counterparts.) Eqs. (19) are closed by determining $u$ from $\omega$ via the Hilbert transform (41)

[TABLE]

The model and closure are natural from a vorticity-formulation viewpoint. The model captures very well features of the teacup flow LH_MMS and exhibits a finite-time singularity Choi_etal_2017 .

Details arise in interpretation of the LH model and the model presented below. These are mostly relegated to Appendix B. The essential point is that away from the cylinder axis $r=0$, the axisymmetric Euler equations with swirl have the same structure as the inviscid 2D Boussinesq equations posed on a half-plane (what I shall refer to simply as the Boussinesq system; see Appendix B). In particular, because the two systems are equivalent on the wall where the singularity occurs, it is convenient to invoke the structure of the simpler Boussinesq system when considering model closures. In the Boussinesq system, closure models for evolution on the boundary come via the Hilbert transform. In discussing the model below, I will continue to use the language of the axisymmetric Euler equations with swirl, but will invoke the equivalence to the Boussinesq system as needed to close the model.

Of the three variables that appear in the LH model (19), two of them, $\omega$ and $u$, are related via the Hilbert transform. One can ask – what about the Hilbert transform of the third variable $\theta$? From (4) and (18a) we have that

[TABLE]

We have used linearity of $H$ and $H(\langle{{b}}\rangle)=H({\rm const})=0$. The final equality is exact, with the understanding that we are invoking the equivalence to the Boussinesq system. (See Appendix B.) Hence the Hilbert transform of $\theta$ is, uniquely, the axial gradient of the pressure field $p_{a}$ on the boundary. This is the unique physical meaning of $H(\theta)$. Hence, even if one did not set out to study the Euler equations in a primitive-variable formation, one is lead to a decomposition of the pressure field just in seeking to understand the meaning of $H(\theta)$. It is important that the variable $\theta$ in the LH model is equivalent to the axial gradient of $p_{a}$. For any model to capture the correct singularity mechanism, it must capture $p_{a}$. The LH model does. This helps to explain why the model can so successfully capture the singularity using only variables on the cylinder wall.

VI.2 Primitive-variable model

The preceding suggests a different approach to closure – working in a primitive-variable formulation and obtaining pressure by Hilbert transform. In the notation of this section, the Euler equations for the axial and swirl flow on the wall are (exactly)

[TABLE]

Recall that $\partial_{z}p_{c}=0$.

These equations can be closed with a single modelling assumption that can be justified from the simulation data. First recall that near the critical ring the contours of the meridional pressure $p_{\rm 2D}$ are nearly circular arcs centred on the critical ring (Fig. 3(b)). If we assume that $p_{\rm 2D}$ is exactly rotationally symmetric about the critical ring in the meridional plane, (hence that its contours are exactly circular arcs), then the meridional pressure is expressible just from source term ${{S}}_{\rm 2D}$ evaluated on the wall. From this the associated axial pressure gradient is

[TABLE]

See Appendix C for details. The superscript $s$ distinguishes this model symmetric meridional pressure from the true meridional pressure gradient. The single modelling approximation we make is to replace the actual adverse pressure gradients $-\partial_{z}p_{\rm 2D}-\partial_{z}p_{b}$ in (22) by the symmetric pressure gradient $-\partial_{z}p_{\rm 2D}^{s}$:

[TABLE]

I address the validity of this approximation below.

The favourable pressure gradient $-\partial_{z}p_{a}$ is given by the Hilbert transform

[TABLE]

This is an exact statement with the appropriate interpretation in terms of the Boussinesq system and requires no other assumptions.

Thus, we arrive at the model

[TABLE]

where $-\partial_{z}p_{\rm 2D}^{s}$ is given by expression (23) and depends only on the axial flow $u$; $-\partial_{z}p_{a}$ is given by expression (25) and depends only on $\theta$, the square of the swirl.

The results from simulations of this model are shown in Fig. 7. The initial condition is

[TABLE]

where the factor $1/2$ is included in $\theta$ because with it the solution almost immediately exhibits scaling behaviour.

Consider first just the case labelled “full model” – meaning the model as written in (26). The dynamics is illustrated in Fig. 7(b) with snapshots of the axial velocity $u$, and the swirl velocity $\pm\sqrt{\theta}$. The slopes of these curves at $z=0$ are $W$ and $\Omega$, the velocity gradients on the critical ring. As the system evolves in time, these gradients steepen from the incoming axial flow.

To establish that the gradients blow up in finite time, we plot $(\Omega/\Omega_{0})^{-1/\gamma}$ and $W^{-1}$ as a function of time in Fig. 7(d). Recall the form of the divergence given in Eqs. (10) and note that the initial condition has $W_{0}=W(0)=-1$. The linearity of these data, together with the common extrapolated zero crossing at $T\simeq 2.250$, is strong evidence that $W$ and $\Omega$ blow up in finite time.

The plot of $(\Omega/\Omega_{0})^{-1/\gamma}$ requires a value for $\gamma$. This can be estimated from the simulation data in two ways. First, the slope of $W^{-1}$ versus $t$ in Fig. 7(d) gives an estimate of $1/\gamma$. A least-squares fit of data over the range $1\leq t\leq 2$ gives $\gamma\simeq 2.31$. The best-fit line is plotted. A second estimate of $\gamma$ is the value that minimises the residual error of a least squares fit of $(\Omega/\Omega_{0})^{-1/\gamma}$ versus $t$. Using the fitting range $1\leq t\leq 2$, the minimum residual error is obtained for $\gamma\simeq 2.31$, the same value to three digits of accuracy. The plot of $(\Omega/\Omega_{0})^{-1/\gamma}$ in Fig. 7(d) uses $\gamma=2.31$ and the corresponding best-fit line is shown. Note this exponent is not very different from the value $\simeq 2.46$ obtained by LH for the full Euler simulation.

Before discussing the implications of the model singularity for full Euler flow, I want to return to the two special cases introduced in Sec. IV. The model contains both. Dropping all the pressure terms, the model (26) reduces to Burger’s equation together with the advection of the passive scalar $\theta$. This system has a finite-time singularity with divergences given in (12). This singularity is shown in Fig. 7(a) and (d). The data have been obtained from a computer simulation that differs from the full model simulation only in the non-evaluation of the pressure terms within the computer code.

The second special case corresponds to the situation with equal axial and radial pressure curvatures on the critical ring: $P=Q$. For the model, $-\partial_{z}p_{\rm 2D}^{s}$ is obtained under the assumption that the pressure field $p_{\rm 2D}$ is exactly rotationally symmetric in the meridional plane about the critical ring (that the contours in Fig. 3(b) are exactly circular arcs). This assumption implies that the the axial and radial curvatures are equal. Hence dropping only $-\partial_{z}p_{a}$ from the model (26), but keeping $-\partial_{z}p_{\rm 2D}^{s}$ gives this special case. The dynamics are shown Fig. 7(c) and (d). Here there is strong deceleration of the axial flow, as seen by the ordering of the $u$ snapshots in Fig. 7(c) compared with the other two cases. The gradient of $u$ at $z=0$ is constant: $W(t)=W_{0}$, as in (14). This system does not blow up. Again the data have been obtained from a computer simulation that differs from the full model simulation only in the non-evaluation of the $p_{a}$ term within the computer code.

The model captures the interplay between inertia and pressure on the cylinder wall. When no stresses are included, the equations blow up in a Burger’s singularity and when only the stress associated with the two dimensional meridional saddle is included, there is no blowup. In both these cases swirl is advected as a passive scalar, leading to vorticity blowing up in the Burger’s case and exponentially growing vorticity in the non-blowup case. Between these two is the case of interest, where the stress from the confinement of swirl on the wall is included. Crucially, here swirl is not a passive scalar – as it is advected toward the critical ring by the axial velocity, it generates increasingly large positive pressure curvature such that the total pressure gradient is insufficient to decelerate in incoming flow. Inertia overwhelms pressure gradients and blowup occurs. The model clearly shows how the teacup effect from the wall swirl is able to drive a finite-time singularity.

VI.3 Connection to Euler

It remains to relate the model to the full Euler flow. The approximation made in the model closure (24) encapsulates the difference between the two, and so I begin with this.

Recall the exact Euler equations on the cylinder wall (22), where the axial pressure gradient $-\partial_{z}p$ separates into contributions $-\partial_{z}p_{\rm 2D}$, $-\partial_{z}p_{b}$ and $-\partial_{z}p_{a}$. These pressure fields are plotted in Fig. 6(b) for the Euler solution at the standard time considered in this paper. The two components with negative curvatures have been combined into $p_{\rm 2D}+p_{b}$. The corresponding adverse pressure gradient $-\partial_{z}(p_{\rm 2D}+p_{b})$ produces outward force decelerating the axial flow approaching the critical ring. Also plotted is $p_{\rm 2D}^{s}$, the pressure obtained from (23) using the actual Euler flow on the cylinder wall. One sees that $p_{\rm 2D}^{s}\leq p_{\rm 2D}+p_{b}\leq 0$, with equality only on the critical ring. This implies that in the vicinity of the critical ring $|\partial_{z}p_{\rm 2D}^{s}|\geq|\partial_{z}p_{\rm 2D}+\partial_{z}p_{b}|$, meaning that the adverse pressure gradient $-\partial_{z}p_{\rm 2D}^{s}$ based on the symmetry assumption provides more deceleration to the incoming flow than the actual adverse pressure gradient $-\partial_{z}(p_{\rm 2D}+p_{b})$.

This justifies the closure approximation (24), where the actual pressure fields acting against singularity formation are replaced by a field that acts more strongly against singularity formation. In other words, the closure approximation suggests that the model should be less liable to blow up than the full Euler equations. This is important if we want to draw inferences about singularities in the Euler equations from singularity formation in the model. We want to know that we have not, at least not in an obvious way, introduced a singularity mechanism through the model closure.

Another way to view the connection between the model and the Euler equations is via the velocity gradient dynamics on the critical ring. Taking the $z$-derivative of model equations (26) and evaluating at $z=0$ gives

[TABLE]

By construction, the curvature of the symmetric pressure $p_{\rm 2D}^{s}$ exactly balances the inertial nonlinearity on the critical ring, leaving only the pressure curvature $P_{a}=-\left.{\partial_{z}^{2}p_{a}}\right|_{z=0}$ driving the velocity gradient $W$.

For actual Euler flow we have instead the inequality

[TABLE]

This follows from (6) under the condition that $Q_{\rm 2D}<P_{\rm 2D}+2P_{b}$. This brings us back to the key observation seen in Fig. 4(b), namely that the meridional pressure curvatures are not exactly equal, $P_{\rm 2D}\neq Q_{\rm 2D}$. From the data in Table 1, $P_{\rm 2D}$ and $Q_{\rm 2D}$ are sufficiently different that the inequality $Q_{\rm 2D}<P_{\rm 2D}+2P_{b}$ holds ($2P_{b}$ is an order of magnitude smaller than the difference between $P_{\rm 2D}$ and $Q_{\rm 2D}$). In fact the inequality follows from the previous observation that $p_{\rm 2D}^{s}\leq p_{\rm 2D}+p_{b}\leq 0$ with equality only on the critical ring.

The difference between the equality in (28) and the inequality in (29) quantifies the previous point that the closure approximation appears to be safe, in that it does not (obviously) enhance singularity formation over that of Euler flow. (The singularities occur with $W\to-\infty$.) For simplicity of discussion, throughout this section I have not strictly distinguished between model and full-Euler quantities. Here it is essential to be clear. Equations (28) and (29) use the same symbols, but apply to different (but closely related) systems – (28) holds for solutions of the model equations (26), while (29) holds for solutions of the full Euler equations. More specifically, (28) holds exactly by construction; (29) holds by numerical observation of the Euler solution and is presumed to hold up to the singularity time.

Although (28) and (29) do not establish a rigorous relationship between the model and Euler flow, they nevertheless reduce the mechanism for singularity formation, in both cases, to its most basic form. Within the Boussinesq analogy, $P_{a}$ is given by the same function of wall swirl in both cases. Using (25) we have $P_{a}\equiv\left.{\partial_{z}^{2}p_{a}}\right|_{z=0}=-H(\partial_{z}\theta)(0)$. Depending on the case, either $\theta$ comes from the solution of the model (26) or else from the swirl on the wall from Euler flow. Referring to (8), for either system to blowup, the pressure curvature due to swirl on the cylinder wall must diverge as $W^{2}$. Specifically, we can relate $W^{2}/\gamma$ to $P_{a}$ in (28) and (29) to give

[TABLE]

Either system will blow up if the left-hand side remains bounded above zero by any finite amount.

In principle, (30) provides a selection mechanism for the exponent $\gamma$. It selects $\gamma$ sharply in the case of the model and bounds $\gamma$ in the case of Euler flow. It is a global condition relating the swirl everywhere on the wall to the velocity gradient on the critical ring. For the model, one can verify numerically that the axial velocity and swirl evolve together such that $-H(\partial_{z}\theta)(0)/W^{2}$ gives the value of $\gamma$. However, this is a triviality given the evidence of a singularity already presented in Fig. 7. Other than numerical simulations, I have been unable to find any convincing arguments or insights into how a particular value of $\gamma$ is selected. I leave this for future work.

VII Conclusion

The potential Euler singularity discovered by Luo and Hou LH_PNAS ; LH_MMS has significantly advanced our mathematical understanding of finite-time singularities and it provides a concrete, easily reproducible case to explore computationally. Here I have sought to understand this singularity from a mechanics point of view and from this gain physical insights into why this particular flow configuration permits velocity gradients to blow up in finite time.

The analysis focuses on the interplay between inertia and pressure. A direct connection is established between the singularity mechanism and flow confinement. The pressure field at the heart of the teacup effect is present solely to confine the rotating fluid within the cylinder; it is determined only by the swirl on the cylinder wall and it plays no role in maintaining incompressibility of the flow. This field is responsible for unequal axial and radial pressure curvatures on the critical ring. This inequality of pressure curvatures is precisely the condition needed for fluid inertia to overwhelm the adverse pressure gradient on the cylinder wall and for velocity gradients to blow up.

To understand how this scenario plays out, a new model has been proposed based on a primitive-variable formulation of the Euler equations. The model describes axial and swirl velocities on the cylinder wall, with closure coming from the dependence of pressure on these velocities. For the swirl the pressure is known exactly. For the axial velocity an approximation is made that has a physical meaning and is supported by Euler simulations. This approximation appears to be distinctly different from those used in related models LH_PNAS ; LH_MMS ; CKY_2015 .

The model captures the interplay between inertia and pressure gradients on the cylinder wall and moreover is embedded in a broader class of problems. In one limit, there are no stresses acting and hence no deceleration of axial flow. This leads to an easily understood Burger’s singularity, accompanied by vorticity blowup from the transport of swirl as a passive scalar. At the other limit, there is only the stress associated with the acceleration of flow around a two dimensional saddle point. This leads to substantial deceleration of the axial flow and no blowup. Between these limits is the case that includes both the stress due to the saddle point and the stress generated from confinement of the swirl on the wall (the teacup effect). Swirl is then not a passive scalar – as it is advected toward the critical ring by the axial velocity, it generates an increasingly large positive pressure curvature (and associated favourable pressure gradient) such that the total pressure gradient is insufficient to decelerate incoming flow. Velocity gradients blow up in a singularity.

There is an important connection between this mechanism and other recent popular models for singularity formation LH_PNAS ; LH_MMS ; CKY_2015 . These models contain two variables, vorticity and square swirl on the cylinder wall. The Hilbert (or similar) transform of the vorticity is used to obtain velocity. The Hilbert transform of the square swirl is, uniquely, the axial gradient of the confining pressure at the core of the mechanism described here.

There are many future directions suggested by this work. Pressure could possibly provide physical insight into the role of the boundary in the rapid growth of vorticity gradients shown by Kiselev and Šverák Kiselev_Sverak_2014 . The model closure proposed here could be connected to the hyperbolic system studied by Kiselev and Tan Kiselev_Tan_2018 . Along these same lines, to impart greater equivalence between axial flow near the cylinder wall and radial flow near $z=0$, one could simulate a cylindrical configuration with a no-penetration condition at $z=0$. The Euler simulations presented here are only for the specific case (initial condition and cylinder aspect ratio) used by Luo and Hou and this leaves open the question of how singularity formation in wall-bounded swirling flows depends on these. One presumes that the scaling exponent $\gamma$ is independent of such factors, as long as blowup occurs, but this too is not presently known since the selection mechanism for the exponent $\gamma$ remains open. It would be highly desirable to investigate these issues and to consider other geometries such as swirling flow within a sphere and to understand the role of pressure in other configurations, such as anti-parallel vortices Bustamante_Kerr . It seems likely that the blowup observed numerically in the model is self-similar, but this is unknown at present, and currently there is no proof of blowup in the model equations. It should be possible to develop precise theorems along the lines of Chae and collaborators chae2008incompressible ; chae2008blow ; constantin2008singular ; chae2010lagrangian ; Chae_etal_2012 to address the specific pressure fields described here. This could possibly lead to a new line of attack on proof of a singularity in the Euler equations. Finally, and most fundamentally, flow confinement is key to the mechanism described here and hence the question remains open as to whether an Euler solution can exhibit blowup in a configuration without a pressure field originating from flow confinement.

Acknowledgements.

This work was partially supported by a grant from the Simons Foundation (Grant number 662985, NG). I am grateful to Guo Luo for pointing out an error in an earlier manuscript and for providing data with which the present simulations could be validated. I became interested in this problem during the IPAM program on the Mathematics of Turbulence and I thank IPAM for their support.

Appendix A: Pressure decomposition

Here we provide details of the pressure decomposition and summarise the relationships that exist between pressure curvatures on the critical ring. In component form, the Euler equations for axisymmetric flow with swirl are

[TABLE]

where $\hat{u}=(u_{r},u_{z})$ and $\hat{\nabla}=(\partial_{r},\partial_{z})$.

Taking the divergence of the nonlinear terms gives the source term ${S}$ on the right-hand-side of the pressure Poisson equation

[TABLE]

The first and third terms are independent of the swirl velocity $u_{\theta}$, while the middle term depends only on $u_{\theta}$. This leads us to define

[TABLE]

Thus the pressure Poisson equation, with boundary condition, is

[TABLE]

This allows for the pressure to be decomposed as $p=p_{\rm 2D}+p_{\rm swirl}$, as given in (15).

Then ${{S}}_{\rm swirl}$ and ${b}$ can be further decomposed into axial mean and fluctuating terms

[TABLE]

where $\langle\rangle$ denotes axial mean,

[TABLE]

This allows for the swirl pressure to be decomposed as $p_{\rm swirl}=p_{a}+p_{b}+p_{c}$, as given in (18).

For the velocity gradient dynamics we require the pressure Hessian $\nabla(\nabla p)$. The pressure field satisfies $\partial_{\theta}p=0$ everywhere. Since $p$ is even in $z$, it also satisfies $\left.{\partial_{z}p}\right|_{z=0}=0$. Hence at $z=0$ the pressure Hessian is

[TABLE]

with the ordering of components $r,\theta,z$. On the critical ring, $\left.{\partial_{r}p}\right|_{c}=0$ since $\left.{b}\right|_{c}=0$, and the pressure Hessian is

[TABLE]

The Laplacian of $p$ is the trace of the Hessian, so on the critical ring $\nabla^{2}p=Q+P$.

From the decomposition, the curvatures for the component fields obey

[TABLE]

There are relationships that hold for the component pressure curvatures on the critical ring. From $\left.{\nabla^{2}p_{\rm swirl}}\right|_{c}=\left.{{{S}}_{\rm swirl}}\right|_{c}=\left.{\frac{1}{r}\partial_{r}u_{\theta}^{2}}\right|_{c}=0$, we have immediately $P_{\rm swirl}+Q_{\rm swirl}=0$ leading to (16).

Less trivial relationships hold for the decomposition of $p_{\rm swirl}$ into $p_{a}+p_{b}+p_{c}$. The reason is that $-\left.{\partial_{r}p_{a}}\right|_{c}=\left.{\partial_{r}p_{c}}\right|_{c}=\langle{{b}}\rangle\neq 0$. Hence these terms appear in the pressure Hessian for $p_{a}$ and $p_{c}$. From (18a), $\left.{\nabla^{2}p_{a}}\right|_{c}=0$, giving

[TABLE]

From (18b) and (18c), $\left.{\nabla^{2}(p_{b}+p_{c})}\right|_{c}=0$, giving

[TABLE]

where we have used that $P_{c}=0$. Note that while the second derivatives of pressure blowup at the singularity, the first derivatives do not. This is because $-\left.{\partial_{r}p_{a}}\right|_{c}=\left.{\partial_{r}p_{c}}\right|_{c}=\langle{{b}}\rangle=\langle u_{\theta}^{2}\left.{\rangle}\right|_{r=1}$, and $u_{\theta}$ does not blowup. Hence, close to the singularity

[TABLE]

where approximately zero means here that the sums are not diverging even though the individual terms are. The simulations suggest that $Q_{c}$ does not blowup at the singularity and can be dropped from (37). This is reasonable since $Q_{c}+\left.{\partial_{r}p_{c}}\right|_{c}=\partial_{r}\langle u_{\theta}^{2}\left.{\rangle}\right|_{r=1}$, and so for $Q_{c}$ to blowup, the gradient of the axial mean must blowup. There is possibly an easy demonstration that this cannot occur.

Appendix B. Connection to 2D Boussinesq system and Hilbert transform

There is a well-known relationship between axisymmetric flow with swirl and two-dimensional thermal convection in the inviscid Boussinesq approximation. See in particular Majda_Bertozzi ; Choi_etal_2017 . The Euler equations for axisymmetric flow with swirl (31) can be recast as,

[TABLE]

where (38a) is a vector equation for the meridional flow $\hat{u}=(u_{r},u_{z})$ obtained by combining (31a) and (31c). The centripetal acceleration term has been written in terms of $(ru_{\theta})^{2}$ and moved to the right-hand side. Equation (38b) is just a reformulation of (31b) into a form that expresses conservation of $ru_{\theta}$ as it is advected as a passive scalar by the meridional flow. Letting $\theta=(ru_{\theta})^{2}$, the equations take the simple from

[TABLE]

In this form, $\theta$ can be viewed as providing a radial driving to the meridional flow. (It should be emphasised, however, that the $\theta$-term in 39a comes from inertia seen in cylindrical coordinates. This term is not associated with stresses acting within the fluid.)

For the inviscid 2D Boussinesq system, consider two-dimensional flow $u(x,y)=u_{x}(x,y)\hat{e}_{x}+u_{y}(x,y)\hat{e}_{y}$ in the region $y\geq 0$. In the Boussinesq approximation, one allows for density variations within the fluid due to thermal expansion from temperature variations. Gravity acts on the density field, here pointing in the $-\hat{e}_{y}$ direction, and the governing equations are

[TABLE]

where $\rho$ represents the density variation relative to some background density. Eq. (40a) describes momentum balance, while Eq. (40b) describes the advection of the density field as a passive scalar. Just as viscosity is zero, thermal diffusivity is zero in this system (both molecular effects are omitted).

The correspondence between the two systems is illustrated in Fig. 8(a). Quoting from (Majda_Bertozzi, , p. 187), “we see that the 2D Boussinesq equations are formally identical to the equations for 3D axisymmetric, swirling flows provided that we evaluate all external variable coefficients” $\dots$ “at $r=1$. Thus away from the axis of symmetry $r=0$ for swirling flows, we expect the qualitative behaviour of the solutions for the two systems of equations to be identical.”

The advantage of the cylindrical system is it is straightforward to simulate numerically. The advantage of the Boussinesq system is that it is mathematically simpler. We may consider either system on a periodic or on an infinite domain in the axial, $z$, or horizontal, $x$, direction. The infinite case is the simplest to consider conceptually and it what I will mean by the Boussinesq system.

This brings us to the Hilbert transformation. Consider a harmonic function $\phi$ in the upper half plane. In our case $\phi$ will be the pressure component $p_{a}$ associated with the boundary swirl, or boundary density for the Boussinesq system. The Hilbert transform of the normal derivative of $\phi$ along $y=0$ is the tangential derivative of $\phi$ along $y=0$. This gives the fundamental relationship between the swirl and axial pressure gradient of $p_{a}$ on the cylinder wall. Hence, the Hilbert transform appears in Sec. VI. A minus sign arises in (21) and (25) because $y$ analogous to $1-r$, so $\partial_{y}=-\partial_{r}$. Concretely, the Hilbert transform of a function $f(x)$ is defined as (Majda_Bertozzi, , p. 173)

[TABLE]

where $PV$ denotes principle value.

Appendix C. Meridional pressure with symmetry assumption

Assume that in a meridional plane a pressure field $p_{\rm 2D}^{s}$ is exactly rotationally symmetric about the critical ring (that the contours in Fig. 3(b) are exactly circular arcs). This assumption necessarily requires invoking the Boussinesq analogy because such a symmetry is impossible within a cylinder. As elsewhere, I nevertheless use here the language of the axisymmetric Euler equations with swirl. Let $(\zeta,\phi)$ be polar coordinates centred on the critical ring as shown in Fig. 8. Then $p_{\rm 2D}^{s}$ is a function only of $\zeta$. We assume $p_{\rm 2D}^{s}$ is determined by a pressure Poisson equation $\nabla^{2}p_{\rm 2D}^{s}={{S}}_{\rm 2D}^{s}$, where the source ${{S}}_{\rm 2D}^{s}$ must also be rotationally symmetric and hence only a function of $\zeta$. Considering the ray $\phi=0$ and identifying $\zeta$ with the positive $z$ axis, we set ${{S}}_{\rm 2D}^{s}={{S}}_{\rm 2D}(r=1,z)$, where ${{S}}_{\rm 2D}(r=1,z)$ is the source term for Euler flow (32) evaluated on the cylinder wall. Straightforward calculation gives ${{S}}_{\rm 2D}(r=1,z)=-2(\partial_{z}u)^{2}$, where $u=u_{z}(r=1,z)$, from which

[TABLE]

Integrating this once gives (23).

Appendix D. Numerical simulations

The Euler equations have been simulated in the vorticity-streamfunction formulation as given by Eqs. (2) in LH_PNAS . The essential difference between the simulations here and those of LH LH_PNAS ; LH_MMS is that here a fixed computation grid is used. A Fourier pseudospectral representation is used in $z$ with dealiasing given by Hou and Li Hou_Li_2007 . A Chebychev grid is used in $r$ with no dealiasing. Fourth-order Runge-Kutta time stepping is used with an adaptive time step such that the CFL number is less than 0.2. Exploiting the separation in the Fourier representation, the Poisson problem for the streamfunction is solved directly. Solving similar Poisson problems, pressure fields are computed in a post-processing step.

For all results reported the computation grid has 769 radial points for $r\in[0,1]$ and 2048 axial points for $z\in[0,L/4)$. At time $t=0.0031$ simulations produce a vorticity maximum $\|\omega\|_{\infty}=1.54276898\times 10^{5}$, agreeing to about 8 digits of precision with the value $\|\omega\|_{\infty}=1.54276901\times 10^{5}$ from simulations by Luo and Hou (private communication).

The simulations of the model equations (26) are mostly straightforward. The $z$ coordinate is mapped to $x\in(-1,1)$ via $\lambda\pi z/2=\tan(\pi x/2)$, where the parameter $\lambda=8$ is used to increase the resolution near $z=0$. The integrals (23) and (25) are computed by quadrature (taking into account the symmetry of the solution). Derivatives are computed spectrally with 2048 equally spaced grid points in $x$. Fourth-order Runge-Kutta time stepping is used with a time step such that the CFL number is fixed at 0.2.