Three-Dimensional Unsteady Axisymmetric Viscous Beltrami Vortex Solutions to the Navier–Stokes Equations

: This paper is aimed at eliciting consistency conditions for the existence of unsteady incompressible axisymmetric swirling viscous Beltrami vortices and explicitly constructing solutions that obey the conditions as well as the Navier–Stokes equations. By Beltrami ﬂow, it is meant that vorticity, i.e., the curl of velocity, is proportional to velocity at any local point in space and time. The consistency conditions are derived for the proportionality coefﬁcient, the velocity ﬁeld and external force. The coefﬁcient, whose dimension is of [length − 1 ], is either constant or nonconstant. In the former case, the well-known exact nondivergent three-dimensional unsteady vortex solutions are obtained by solving the evolution equations for the stream function directly. In the latter case, the consistency conditions are given by nonlinear equations of the stream function, one of which corresponds to the Bragg–Hawthorne equation for steady inviscid ﬂow. Solutions of a novel type are found by numerically solving the nonlinear constraint equation at a ﬁxed time. Time dependence is recovered by taking advantage of the linearity of the evolution equation of the stream function. The proportionality coefﬁcient is found to depend on space and time. A phenomenon of partial restoration of the broken scaling invariance is observed at short distances.


Introduction
Three-dimensional swirling vortex solutions to the Navier-Stokes equations have been of researchers' particular interest because of their immediate relevance to phenomena in nature, laboratory and practical engineering.Besides numerical solutions, exact analytical solutions should make it possible to deepen the understanding of the details of the structure of vortices characterized by speed, vorticity, pressure, shearing, stretching, stagnation, etc.
In solving the non-linear Navier-Stokes equations, it is a common practice to assume the functional form of a solution that reduces the original partial differential equations to a more tractable system of ordinary differential equations.The well-known vortex solutions found in this way involve the steady axisymmetric one-cell solution by Burgers [1] and Rott [2], the two-cell solution by Sullivan [3], and the unsteady multi-cell solutions by Bellamy-Knights [4,5].Solutions with a line singularity were presented in [6,7].Solutions which temporally approach the Burgers' solution were studied in [8].
It was not clear if the unsteady solutions in [4,5] were linked to the steady solutions in [1][2][3].Recently, these steady and unsteady solutions were treated in a unified way in [9], wherein the "instability" of these solutions were shown to be attributed to their extreme sensitivity to boundary conditions.For other vortex solutions and a review, see Drazin and Riley [10].
All the solutions given above are similarity solutions that have been obtained by taking advantage of the invariance of the Navier-Stokes equations under the scaling transformation.As a consequence, some components of their velocity fields linearly diverge at spatial infinity.Hence, the solutions are regarded as modeling flow patterns near the symmetry axis or plane.
Observable physical excitations on a uniform background flow must have a finite energy density, except for logarithmic singularity [9].In particular, the divergence at spatial infinity will be problematic when we are interested in investigating the global nature of a vortex.The primary interest in this paper is to find three-dimensional unsteady viscous swirling vortex solutions whose energy density is globally not divergent.(In this paper, by "three-dimensional", we mean that there is no Cartesian coordinate system in which one of the velocity components identically vanishes).
Beltrami flow, whose vorticity is parallel to velocity in the whole space under consideration, is expected to be a promising candidate for the above purpose because of the existence of the proportionality coefficient, a scalar quantity of dimension [length −1 ], thereby raising the possibility of generating vortices of finite extension in space and/or time.Incidentally, we recall that tilting the vorticity is also accompanied by the breaking of the scaling invariance [8].
The marked characteristics of the Beltrami flow, the study of which was originated by Trkal [11] and Gromeka [12], is that the vorticity equations are apparently linear by virtue of the vanishing vector product of vorticity and velocity.In particular, inviscid Beltrami flows obeying the Euler equations have been of researchers' intense concern because of their close relevance to the flows with large Reynolds numbers that are expected to model the flows in meteorology [13,14] or engineering [15][16][17].The so-called ABC flow that shows rich physical content in the Lagrangian dynamics provides a key to understanding chaotic or dynamo phenomena in fluids [18][19][20][21].
Steady axisymmetric inviscid Beltramian solutions in bounded or unbounded space have been found by focusing on the so-called Bragg-Hawthorne equation [17,22,23].How the structure of the fundamental modes admitted by the linear equation varies by depending on the coordinate system chosen for solving the equation has been revealed in [22,24].Non-axisymmetric solutions with vanishing helicity, i.e., inner products of velocity and vorticity, as well as the vortices with a constant proportionality coefficient, c hereafter, are presented in [24].A Beltrami flow whose c is constant is called a Trkalian flow.
Despite the seemingly quite artificial setting for the Beltrami flows, the physical mechanism of generating a Beltramian (i.e., force free magnetic) field was considered as early as about 70 years ago by Chandrasekhar [25] and Chandrasekhar and Woltjer [26] in an astrophysical context with a suggestion of a principle of minimum dissipation and maximum energy.Let us recall that, with the correspondence of magnetic field to vorticity, the force-free motion of charged particles in magnetic field has mathematical similarity to Beltrami flow.The idea in [25,26] parallels the expectation in the Newtonian fluid dynamics that an increase in the Reynolds number tends to align velocity and vorticity [27,28].Later, the importance of the maximization of helicity [21] and the minimization of magnetic energy due to dissipation in the evolution of magnetic field of nontrivial topology has been pointed out in [29,30].It is convincing to explore the role of Beltrami flows in meteorology [31] and plasma physics [32,33] where strong flows participate.There exist profound physical meanings in studying Beltrami flows.
It is known that a steady incompressible viscous Beltrami flow is possible when external force is nonconservative.Steady compressible solutions were presented in [34].Examples of unsteady Beltrami flow with a constant c were given by [11,35].That threedimensional unsteady incompressible viscous Beltrami flows with a constant c decay exponentially in time is also argued by [36,37].These solutions are specified by, in addition to c, a wavevector that describes spatial oscillations of the flow.In this paper, we shall call such a flow (solution) a single mode flow (solution).
The works by [11,[35][36][37] assert that the viscous Trkalian velocity field of a single mode in an unbounded space exponentially decays in time with the characteristic time proportional to 1/νc with ν being the kinematic viscosity.The appearance of this time scale is due to the parameter c of the dimension [length −1 ], thereby violating the scaling invariance of the original Navier-Stokes equations.With a characteristic length scale being defined, it is thus expected that the spatial extension of vortex solution may also be rendered finite.The generalized Beltrami flows [10,38,39] and the extended Beltrami flows [40,41] will share the same features in common.
An axisymmetric genuine meridional Beltrami vortex does not exist (see the next section).Local alignments of velocity and vorticity, meanwhile, are observed in numerical simulations of turbulence [27,[42][43][44].A constraint on the time evolution of a spatially periodic flow from an arbitrary initial condition was given in [45].That the incompressible Beltrami flow decays timewise exponentially if c does not depend on time was shown in [46].These works naturally let us nourish an expectation that the Beltrami flows with the spatially variable proportionality coefficient exist.If that is the case, recalling that the decay time is 1/νc, we may address a question: how is the anticipated "rapid" decay of the flow avoided especially for small Reynolds numbers?To the author's knowledge, any works on this problem have not been reported thus far.
Our plan in this paper is as follows.We derive and solve the evolution equation and the constraint equation for the stream function and the coefficient c on the projected meridional motion, i.e., a projection of the three-dimensional flow onto the local meridian plane in the cylindrical coordinates.The derived constraint equation corresponds to the Bragg-Hawthorne equation that has been frequently used for studying steady inviscid flows.However, our constraint does not depend on viscosity.An obstacle is that the constraint equation, although it reduces to the Bragg-Hawthorne equation in a special case, is generally nonlinear and is difficult to be solved analytically in the whole space and time.However, since the evolution equation is linear and the temporal behaviors of single mode solutions are uniquely determined, we can construct the entire time-dependence by superposing the single mode solutions obtained by mode-decomposition of a numerical solution of the constraint equation.Of course, it must be assured that the superposed flow is Beltramian.
The outline of the paper is as follows.After the implication of the axisymmetric viscous Beltrami vortex is reviewed and summarized, we derive the constraint equation for unsteady three-dimensional axisymmetric viscous swirling Beltrami vortices.It is proved that, if the velocity field is time-dependent and c is variable, c is time-dependent, too.Finally, the well-known exact vortex solution for constant c and a new solution with nonconstant c are constructed.

Dynamical and Constraint Equations for the Axisymmetric Beltrami Vortex
The Navier-Stokes equations for the respective components of the velocity field in cylindrical coordinate (r, θ, z) are expressed as where v = (v r , v θ , v z ) T is the velocity, ν the kinematic viscosity, p the pressure, ρ the density, and F = (F r , F θ , F z ) T the external force.In this paper, we restrict ourselves to axisymmetric flows so that all the θ-derivative terms in (1)~(3) are ignored.
A Beltrami flow is one whose velocity and vorticity are parallel: where c may generally be a function of time and spatial coordinates.Hereafter, we use a symbol ∂ α for partial derivative with respect to the variables α = r, z.Note that, if v θ ≡ 0, then v ≡ 0; that is, a genuinely meridional Beltrami flow does not exist.Note that rv θ plays the role of the stream function for the meridional components v r and v z .Employing the parallelism of velocity and vorticity, the variables v r and v z in the second row of (4) can be eliminated to give a kinematic relation between v θ and c: Furthermore, the incompressible flow must fulfill the continuity condition Since ω given by ( 4) is divergence-free, too, these conditions lead to another one on c: ∂ r c and ∂ z c are determined by ( 5) and (7), respectively, as where Equation ( 8) implies that, if c is nonzero and constant, then Λ = 0, provided that ω r = 0 = ω z .If c depends on either r or z only, either ω r or ω z may vanish.

Consistency of Dynamical and Constraint Equations
The two constraints (8) for the partial derivatives of the (pseudo) scalar c are generally incompatible with each other.Consistency of (8) for arbitrary regular solution requires the partial derivatives ∂ r and ∂ z to commute.From (8), a straightforward calculation yields Let φ denote any function that satisfies ω • ∇φ = 0.Then, the function Λ given by is a solution to (10).In terms of such a function φ , Equation ( 9) is expressed as On the other hand, employing (11) and referring to (4), the conditions (8) read and Equations ( 13) and (14) imply that c will depend on the spatial coordinates via the stream function ψ ≡ rv θ only.Notice that ∂(φ , ψ)/∂(r, z) = 0 for regular c so that φ and ψ are mutually dependent when the spatial variations are concerned.We thus may write with a relation φ = ∂φ(ψ, t)/∂ψ.Namely, the conditions ( 12) with ( 15) are equivalent to (8).Note that ψ can be time-dependent, too.Recalling the definition (4) of ω, (12) is written as This equation elicits the hidden nonlinearity of the dynamics of the axisymmetric Beltrami vortex when φ = 0 and will henceforth play the central role in our arguments.Equation ( 16) governs the spatial variation of the stream function ψ and in this sense corresponds to the Bragg-Hawthorne equation for inviscid flow [47,48]: where ψ is the stream function (not necessarily equal to rv θ ), H ≡ p/ρ + v 2 /2 is called the (pressure or energy) head, and C ≡ rv θ (= ψ) is supposed to be a function of the stream function ψ for the axisymmetric inviscid flow.Note that ψ and C are not related a priori.For recent applications of the Bragg-Hawthorne equation, see, e.g., [17,22,23].The Bragg-Hawthorne equation is an expression of the equation of motion for inviscid flow under steadiness assumption.Therefore, the solution of this equation is a solution of equation of motion.On the other hand, our Equation ( 16) has been derived not from the dynamical equations but from the constraints on the unsteady incompressible axisymmetric Beltrami flow.There is no immediate reason that these two equations should coincide with each other except the same derivative terms that result from taking the velocity multiplied by the coordinate as the dependent field.Nevertheless, two equations become essentially identical when they are linear.In fact, assuming that dH/dψ = 0 and C ∝ ψ, the Bragg-Hawthorne equation writing dC/dψ as c (being constant) reduces to ( 16) with φ = 0.
Because the advection term is the vector product of vorticity and velocity added by the gradient of the specific kinetic energy, the Navier-Stokes Equation (2) for axisymmetric (i.e., ∂ θ = 0) Beltrami flow is linearized as Equivalently, from ( 16) and ( 17), we have with lnc = φ.This is the dynamical equation for ψ that expresses the effect of the dissipation as well as the torque rF θ on the rate of temporal change of ψ.Noting that ψ is the angular momentum component per unit mass about the symmetry axis, the first term with a negative sign in the large braces of ( 18) represents the direct resistance on the fluid element.This resistance term, being proportional to c 2 and an invariable concomitant of alignment of velocity and vorticity, reveals the hallmark of the Beltrami flow.When ψ is steady and has negligible spatial variations, the direct resistance is balanced with the torque.The second term involves the mixed effects of the rate of strain and vorticity.The constraint (16) Is written in terms of the velocity component as When c is constant and φ = 0, it is obvious that the fundamental spatial length scale of the velocity field is given by c −1 .Many analytic Beltrami flow solutions are known for constant c.It remains to examine in what way the conditions found so far for unsteady axisymmetric Beltrami flows are consistent with the rest of the Navier-Stokes equations, i.e., (1) and (3).The details are relegated to Appendix A. One condition is that 1/c(ψ, t) is integrable over ψ.The components F r and F z are related to c through (A2) and (A3) in Appendix A. Then, we have the second condition When c is constant, (20) can be decomposed to H must be spatially constant if external force is absent.Suppose that c is constant.If either the external force is derived from an axially symmetric potential V or F θ = F θ (r) ∝ 1/r, then ( 21) and ( 22) immediately lead to Bernoulli's law Suppose that v θ of a Beltrami flow that obeys the Navier-Stokes Equation ( 2) is known.Notice that, in solving (2), information of the pressure is not needed because of the system's axisymmetry.It has been shown in Appendix A that the Navier-Stokes Equations ( 1) and ( 3) are transformed to (A2) and (A3), and the equations of χ ≡ ψ c −1 dψ.v r and v z derived from χ as described in Appendix A are automatically guaranteed to fulfill the continuity condition as well as the Beltrami relation (4).Equations (A2) and (A3) implicitly involve p/ρ through the head H. Therefore, what we need to obtain v r and v z from χ, besides the boundary conditions which are compatible with the Beltrami relation (4), are H and the external force F related by (20) or ( 21) and (22).Once each component of v is determined as a function of r, z, and t, c is obtained from the Beltrami relation.That this procedure is possible is a peculiarity of the axisymmetric Beltrami flow.F will be present when c is not constant.The details on the relation of nonconstant c and F are elaborated in Appendix B.

Vortex Solutions
In this section, we consider viscous Beltrami vortices with constant c (Trkalian flow) and nonconstant c separately.The Beltrami flows, inviscid or viscous, with constant c were studied in many works, as was cited in Introduction.The case in which c is variable is studied formally in [46] for axially asymmetric flows in a Cartesian coordinate system.In particular, if c has no explicit time dependence and the external force is conservative and the Bernoulli function is constant, then v 2 decays timewise exponentially.The generalized ABC flow solution for compressible flow was also found for c that is dependent on a single spatial variable in [34].

Constant c
We first solve the dynamical equations with the constraints of constant c and F θ = 0. From ( 13) and ( 14), c can be constant when ψ is constant or v θ ∝ 1/r, too.However, this does not satisfy (16).
In case ψ is not constant, constant c implies φ = 0.Then, ( 19) is linear and homogeneous.It is easily solved for a single mode flow by the separation of variables: where J 1 is the Bessel function of the first order, α(t) a function of t to be determined later and δ a phase.k and δ are both real or pure imaginary.If we require the solution to be nondivergent in the whole spatial region, k must be real and |k| < c.Single mode solutions of this type frequently appear, irrespective of the symmetry when the nonlinear terms are somehow removed [22,35].See also [49], wherein the Coriolis force is taken into account.From ( 23) and ( 4), it follows that Referring to (4), it is readily confirmed that the equality ω = cv holds.Equation (18) gives Since c = 0, this equation has a nontrivial solution The time-dependent velocity fields given by ( 23) and ( 24) then fulfill the Navier-Stokes Equations ( 1)~(3).As was expected, solutions exist whose kinetic energy density v 2 /2 is finite in the whole space.
The characteristic decay time is given by 1/νc 2 .Only steady vortex solutions are allowed in the inviscid limit.Separation of the variables r and t in the above solution is a consequence of breaking the scaling invariance of the Navier-Stokes equations.Exact three-dimensional solutions with the decay-time inversely proportional to the kinematic viscosity have been presented by several authors [11,35,36,50].

Nonconstant c
The proportionality coefficient c will be variable when F is nonconservative or H + V is not spatially constant.Unfortunately, we do not know the functional form of c(ψ, t) or φ(ψ, t), whose temporal evolutions will be governed by the dynamics.However, it may be sensible to set their instantaneous form as an initial condition.We consider a simplest case given by φ(ψ, t) where A and B, whose dimensions are [length −1 ] and [length −1 •velocity −1 ], respectively, are spatially constant (but possibly time-dependent).The assumption (26) amounts to keeping the first two terms in the power series expansion of φ in ψ.With this choice, the scales of the spatial extension and velocity of the solution will be given by A −1 and A/B, respectively.The turn-over time of the vortex rotation will be the order of B/A 2 .The representative Reynolds number is given by Re = 1/νB.In passing, the assumption of time-independent c employed in [46] and the form (26) are compatible with each other for a single mode solution whose stream function decays as e −νc 2 t .In this case, the parameter B may be chosen as to be proportional to e νc 2 t to cancel the decay factor e −νc 2 t .Investigating whether this is possible or not for axisymmetric flow is beyond the scope of the present paper.
We begin with exploring the spatial dependence of the velocity fields with the time fixed.Since φ = B, we rewrite Equation ( 16) as Changing the variable by Bψ = − ln|g|, Equation ( 27) is put into the form The analytic solution of ( 28) is not known.Since the term A 2 g −1 ln|g| results from a variation of U(g) = A 2 (ln|g|) 2 /2, we may approximate it around the minimum point g = 1 by a harmonic potential U(g) ≈ A 2 (g − 1) 2 /2 and linearize (28) as where ĝ ≡ 1 − g ≈ Bψ.Equation ( 29) is identical to (19) with φ = 0. Therefore, when | ĝ| is sufficiently smaller than unity, the solution to ( 28) is given by ( 23), with c replaced by A. If v θ decreases as fast or faster than 1/r for r → ∞ , the flow is approximately given by ( 23) and (24).The velocity decays exponentially with time.Equation ( 27) permits a variety of solutions.Shown in Figure 1 are three examples, (i), (ii), and (iii), obtained numerically under the condition that ∂ z ψ = ∂ 2 z ψ = 0 .That is, these solutions have no z-dependence.Each solution for v θ has one peak and decays gradually at long distances.Equations ( 28) or (29) show that 1/|A| defines the system's length scale, so that, grossly speaking, the (half) width of the peak becomes larger with the decrease in A. From Figure 1, the radial distance to the peak and the azimuthal velocity at the peak are grossly read off as 0.5A −1 and A|B| −1 , respectively.Taking these as the rules of thumb, the turn-over time of the vortex at the peak is given grossly by 2π(0.5A−1 )/(A B −1 ) ∼ O( B /A 2 ) , as was already expected.The existence of the solutions (i), (ii), and (iii) implies the existence of an infinite number of continuous sequences of solutions in accordance with the continuous change in the parameters A, B, and the boundary condition ∂ r ψ(r 0 ).It is easily confirmed by substitution that, ignoring the last term on the l.h.s. of ( 27), the solution to (27) for large r and fixed z is approximated by ψ ≈ ψ 0 ≡ −(2/B) ln r with negative B. The tail of the solution is therefore approximated by v θ ≈ −(2/B) ln r/r, which, from (26), implies that c is nearly proportional to 1/r 2 at long distances.Vortex solutions with such a long tail were numerically confirmed to exist up to r = 30 for |A| ≤ 6.84 • • • .We disregard the case of larger |A|, for which, although not shown here, numerical calculations exhibited that the solutions keenly oscillate around z = 0.
As regards the behavior in the finite z region, Equation (27) suggests that nondivergent ψ will be expanded in a power series of e −|A|z .By substituting the expression to (27) and comparing the coefficients of e −n|A|z , the equations of ψ n are, to n = 2, where c 0 ≡ e ψ 0 for the choice (26) for c(ψ).The prime stands for a differentiation with respect to r.The first equation in (31) gives the solutions shown in Figure 1.Specifically, ψ 0 ∼ −2 ln r for large r.
Splitting the term A 2 c 2 0 ψ n as the sum of linear and nonlinear terms by ψ n and discarding all the nonlinear terms, the equations in (31) for n ≥ 1 may be linearized as The n 2 emerges from ∂ 2 z in (27).This means that ψ n /r obeys the Bessel differential equation.Thus, restricting ourselves to regular solutions to (31), ψ n will behave as 1/2 A r) both near r = 0 and at medium distances provided that the nonlinear terms are negligible.At very large distances, the term involving the derivative of ψ 0 on the r.h.s. will have to be retained.Thus, for the leading form of the solutions, we solve the equation Regular solutions exist and are given by r −1 J 1 ((n 2 + 1) 1/2 |A|r).
The expected exponential behavior in z specified by (30) can be observed analytically in the vicinity of the symmetry axis, as is shown below.
Near the z-axis, v θ behaves as v θ ∝ r.With the time being fixed, let us put for r ≈ 0. Substituting (32) to (27), we have where prime stands for a derivative with respect to z. Referring to Figure 1, it will not be unreasonable as a first approximation to ignore the O(r 4 -term at the very vicinity of z-axis. Multiplying ∂ z a to the both sides of ( 33) and integrating once, we have where C is an integration constant, which is generally time-dependent.Solutions to (34) may be expressed in terms of Jacobi's elliptic function or are easily found via numerical integration.Here, we give an exact solution that is constructed by elementary functions for a special combination of the parameters, i.e., 48B 2 C = A 6 .In this case, the polynomial in the square root in ( 34) is factorized to (8B/3)(a + q) 2 (q/2 − a), q = A 2 /4B.Using the integration formula with α = q = A 2 /4B and β = q/2 = A 2 /8B, we have where The numerical factor in front of e Az in (36) has been chosen for a convention a(0) = 0.
(Recall that the origin of the coordinate z is arbitrary.)As was already expected, a(z) has a power series expansion in e −|A|z and converges exponentially to a constant value as z → ∞ near the symmetry axis.
The quantity 4a(z)/A 2 as a function of z was numerically calculated, and the result is shown in Figure 2 for four cases: A = −10, −5, 10, and 5.When A is positive, |a| monotoni- cally increases with z, while a negative A gives rise to an additional zero, z 0 .Numerical calculation gives |A|z 0 = 2.6391 • • • , which is independent of |A| because a(z) is a func- tion of Az only.At this point v z vanishes.In each case depicted in the figure, 4a/A 2 exponentially approaches a constant, either 1 (exact value) or 0.101020 Note also that, when A is negative, there is a z at which da/dz = 0 or v r = 0 from (4).The distance of the stream line from the symmetry axis minimizes at this z and then increases to infinity as the stream line approaches the plane z = 0 or z 0 .The case of negative A may be suited to describing a flow between two parallel planar boundaries.As regards the solutions depicted in Figure 2, the flow directing toward the symmetry axis just below the plane z = z 0 forms a downdraft (i.e., a(z) < 0) near the axis and then changes the direction to far infinity near z = 0.The flow above the plane z = z 0 forms an updraft (a(z) > 0).These characteristics of the flows will be altered if the solutions are extrapolated to negative z.Note also that, when A is negative, there is a z at which 4).The distance of the stream line from the symmetry axis minimizes at this z and then increases to infinity as the stream line approaches the plane 0 z = or 0 z .The case of negative A may be suited to describing a flow between two parallel planar bounda- ries.As regards the solutions depicted in Figure 2, the flow directing toward the symmetry axis just below the plane

Time Dependence
All of the calculations of the velocity field so far are at some fixed time.In order to determine the time-dependence of the solutions, we start from the fact that the time-dependent solution to the linear equation of motion (17) is generally expressed as a superposition of the single-mode solutions as where ( , , ) with the integration volume element d rdrdzdt σ = . The integration regions are such that

Time Dependence
All of the calculations of the velocity field so far are at some fixed time.In order to determine the time-dependence of the solutions, we start from the fact that the time-dependent solution to the linear equation of motion (17) is generally expressed as a superposition of the single-mode solutions as ψ = r v(q, k, t)J 1 (qr)e ikz qdqdk + r G 0 (σ; σ )F θ (σ )dσ , v(q, k, t) = f (q, k)e −ν(q 2 +k 2 )(t+t 0 ) , (37) where σ ≡ (r, z, t) and v(q, k, t) is the Fourier-Hankel transform of v θ = ψ/r with the integration volume element dσ = rdrdzdt.The integration regions are such that 0 ≤ The retarded Green's function G 0 (σ; σ ) satisfies, together with (A6), and is explicitly given in Appendix B. Note that v(q, k, t) has a Gaussian damping factor with respect to q and k.At the same time, it exhibits a temporally exponential decay.Spectrum function f (q, k) is the remaining factor of v.The v(q, k, 0) will be determined from the Fourier-Hankel transformation of v θ obtained (numerically) by solving (27) or (28) for an unknown t 0 and by extracting the Gaussian damping factor from the solution at large q or k.Finally, the time dependence is recovered by substituting the damping factor exp(−ν(q 2 + k 2 )(t + t 0 )) as in (37).
Two hypothetical cases will be instructive.If f (q, k) distributed on a "circle" q 2 + k 2 = c 2 , (37) would be a superposition of (25), allowing the possibility of pure imaginary q and/or k.If f (q, k) depended on q as q µ−1 , the integration over q in (37) would be performed to give the r-dependence v θ ∝ (r/(νt) µ/2+1 ) 1 F 1 (µ/2 + 1; 2; −r 2 /4νt), where 1 F 1 is the confluent hypergeometric function, whose limiting behavior is given by Note that, in the above example, the combination of r and t as the similarity variable r 2 /νt emerges in the velocity field.This phenomenon results from some distribution of the wavenumber and indicates the partial restoration of the scaling invariance that had once been broken by setting c as constant.
The procedure of extracting the Gaussian factor and the spectrum function f from the data for v θ of the solution (i) in Figure 1 is sketched in Appendix C, and the result is shown in Figure 3.In the present case, the q-dependence of v depicted in Figure 3 indicates f (q) → 1 for 20 < q 2 < 40.The spectrum function of the solution (i) is prominent at small q or long wavelengths and exhibits a Gaussian damping at short wavelengths.
The spectrum function shown in Figure 3 is employed to recover the time-dependence of the solution (i) in accordance with (37), with the result shown in Figure 4.In the numerical calculations, the upper bound of each integration was taken to give the tenth zero of J 1 in the integrand.Therefore, calculational errors accumulated near r = 0.It is notable that the resultant solution temporally behaves similarly to the familiar scale invariant solutions [4,5,9].See also Figure 5.17 in [10] for the (two dimensional) scaling invariant Oseen's flow.
Unfortunately, the ψ constructed in this way does generally not satisfy the constraint Equation ( 16) for t > 0 because the constraint equation is entirely independent of the equation of motion.Thus, some temporal changes in the tilt of vorticity from the one in the Beltrami vortex will inhere in the ψ given by (37).
However, we can always choose a continuous and differentiable sequence of solutions to the constraint equation.This is because the differential Equation ( 27) is not singular.Let the sequence start from, e.g., the solution (i) in Figure 1.We specify each solution on the sequence with a labelling parameter t 0 and the spectrum function f (q, k; t 0 ).Shifting the parameter as t 0 → t + t 0 with t 0 to be used for the solution (i), the solution t to the constraint equation is written as The suffix C means that it is the solution to the constraint equation.The t-dependence of f c reflects the t-dependences of the parameters A, B, and ∂ r ψ(r 0 ), too.Unfortunately, the ψ constructed in this way does generally not satisfy the con- straint Equation ( 16) for 0 t > because the constraint equation is entirely independent of the equation of motion.Thus, some temporal changes in the tilt of vorticity from the one in the Beltrami vortex will inhere in the ψ given by (37).Unfortunately, the ψ constructed in this way does generally not satisfy the con- straint Equation ( 16) for 0 t > because the constraint equation is entirely independent of the equation of motion.Thus, some temporal changes in the tilt of vorticity from the one in the Beltrami vortex will inhere in the ψ given by (37).
However, we can always choose a continuous and differentiable sequence of solutions to the constraint equation.This is because the differential Equation ( 27) is not singular.Let the sequence start from, e.g., the solution (i) in Figure 1.We specify each solution On the other hand, the solution to the evolution equation with the initial condition chosen as, e.g., solution (i) in Figure 1, is given by (37): where the suffix D stands for being the solution to the dynamical equation.The above expression for ψ D is obtained by the Green function method for the operator ) similar to the one described in Appendix B. Now, identifying the labelling parameter t in ψ C (r, z; t) with the physical time, we ask a question: Under what condition does the ψ C (r, z; t) coincide with the ψ D (r, z; t)?
This question seems reasonable because of the similarity of Figure 1 to Figure 4.If ψ C (r, z; t) = ψ D (r, z; t), operating ∂ t − ν∇ 2 to (ψ C (r, z; t) − ψ D (r, z; t))/r from the left and using the identity (B3), we have F θ (r, z, t) = ∞ 0 ∂ t f C (q, k; t + t 0 )e −ν(q 2 +k 2 )(t+t 0 ) J 1 (qr)e ikz qdqdk.(38) This is the (necessary) condition for the Beltrami vortex to continuously exist.The procedure of constructing ψ C described above implies that f C (q, k; t + t 0 ) is t-dependent and the r.h.s. of (38) does not generally vanish, so that F θ does not vanish either.The analysis that has led to this consequence is based on the assumption (26) for c(ψ, t).Whether a sequence whose f c does not depend on t exists for a more general form of c(ψ, t) is an open question.In summary, an axially symmetric incompressible viscous swirling Beltrami vortex with variable c is possible when an external force is present.

Conclusions and Discussion
In the axisymmetric viscous swirling Beltrami vortices, the proportionality coefficient c between the velocity and vorticity is constrained by the equation analogous to the Bragg-Hawthorne equation for steady inviscid flow.Our consistency equation, applicable also to unsteady flows, is nonlinear in the stream function except for the case of a constant c.The exact solutions found for a constant c are of the type known up to now for inviscid flow [22,49].
One of the main results for the case of a nonconstant c is that the spatial and temporal dependence of c emerge through the stream function ψ = rv θ , i.e., c = c(ψ, t).The physically meaningful numerical solution sought in the simplified form of c(ψ, t) exhibited a cyclone-like behavior for radial motion and a gradual temporal decay.The velocity field does not diverge at infinite distances, as has been shown analytically or numerically in the present study, in which ln|c| is taken to be linear in ψ.
Breaking the scale invariance due to c of the dimension [length −1 ] generates a characteristic length l = c −1 and a decay time τ = l 2 /ν for a single-mode flow.It should be noted that the decay time does not depend on the wavenumber, k, appearing in (23).For air at normal pressure and room temperature, τ air ≈ 5 × 10 4 (l/m) 2 s.Similarly, for water, we have an estimation τ water ≈ 10 6 (l/m) 2 s.These periods shall be long enough for a macroscopic Beltrami vortex to persist.We also saw that a nonconstant c will result in a longevity of the vortex when the spectrum function obeys a power law.
The second important finding is that an axisymmetric Beltrami vortex with a nonconstant c continues to be Beltramian only when some special kind of external force depending on the consistency equation must be at work.For example, the F θ acts to directly grow the azimuthal component of the flow and induce the meridional components of the motion via the stream function.This is why the vortices undergoing such external forces as the Coriolis force or the Lorentz force are worth studying.The vortex-boundary and vortex-vortex interactions will also play the role.The numerical simulation method is expected to provide powerful tools to elicit information on this issue [27,[42][43][44].
Without the first principle to determine the functional form of c(ψ, t) or C, in the present paper, we have examined only the simplest form of the consistency condition for the viscous vortices.Numerous types of consistency conditions remain to be explored.
Funding: This research received no external funding.

J 2023, 6 ,
FOR PEER REVIEW 11 has a power series expansion in | | e A z − and converges exponentially to a constant value as z → ∞ near the symmetry axis.as a function of z was numerically calculated, and the re- sult is shown in Figure 2 for four cases: 10 A = − , 5 − , 10 , and 5 .When A is positive, | | a monotonically increases with z , while a negative A gives rise to an additional zero,  , which is independent of | | A because ( ) a z is a function of Az only.At this point z v vanishes.In each case depicted in the figure, approaches a constant, either 1 (exact value) or 0.101020 , for z → ∞ .
forms a downdraft (i.e., ( ) 0 a z < ) near the axis and then changes the direction to far infinity near 0 z = .The flow above the plane 0 z z = forms an updraft ( ( ) 0 a z > ).These characteristics of the flows will be altered if the solutions are extrapolated to negative z .

Figure 3 .
Figure 3. Upper panel: The numerical solution (i) in Figure 1 is shown by red circles and red thick curve.Solid black curve is the result of fitting to the numerical solution (i) by equation (A9) in Appendix C. Lower panel: Hankel transform v  (red solid curve), Gaussian factor

Figure 4 .
Figure 4. Temporal variation due to (15) with 0 F θ = of the solution (i) for θ v in Figure 1.Integra-

Figure 3 .
Figure 3. Upper panel: The numerical solution (i) in Figure 1 is shown by red circles and red thick curve.Solid black curve is the result of fitting to the numerical solution (i) by Equation (A9) in Appendix C. Lower panel: Hankel transform v (red solid curve), Gaussian factor exp(−γq 2 ) with = 0.28 and γ = 0.048 (green broken line) and f = v(q, 0)exp(γq 2 ) (blue dotted curve).

Figure 3 .
Figure 3. Upper panel: The numerical solution (i) in Figure 1 is shown by red circles and red thick curve.Solid black curve is the result of fitting to the numerical solution (i) by equation (A9) in Appendix C. Lower panel: Hankel transform v  (red solid curve), Gaussian factor

Figure 4 .
Figure 4. Temporal variation due to (15) with 0 F θ = of the solution (i) for θ v in Figure 1.Integra-

Figure 4 .
Figure 4. Temporal variation due to(15) with F θ = 0 of the solution (i) for v θ in Figure1.Integration was performed up to the tenth zero of J 1 .