1. Introduction
Swirling flows are employed in a wide range of engineering applications such as hydraulic turbines (Zhang et al. Reference Zhang, Mao, Wu, Chen, Wu and Liu2009), jet engines (Musa et al. Reference Musa, Weixuan, Xiong, Lunkun and Wenhe2018), chemical reactors (Yang et al. Reference Yang, Holemans, Lagrain, Sels and Vanierschot2023), mixing devices (Keeton et al. Reference Keeton, Nomura, Sánchez and Williams2023) and other applications. It is important to understand the dynamics of unstable swirling flows and the complicated evolution to vortex breakdown. Reviews of the vortex breakdown phenomena by Leibovich (Reference Leibovich1984), Escudier (Reference Escudier1988), Rusak & Wang (Reference Rusak and Wang1996) and Rusak (Reference Rusak2000) show various types of vortex breakdown such as bubble, spiral type and the temporal transition between them. Cary, Darmofal & Powell (Reference Cary, Darmofal and Powell1997), Tromp & Beran (Reference Tromp and Beran1997), Billant, Chomaz & Huerre (Reference Billant, Chomaz and Huerre1998), Dennis, Seraudie & Poole (Reference Dennis, Seraudie and Poole2014), Meng, Liu & Luo (Reference Meng, Liu and Luo2018) and Vanierschot et al. (Reference Vanierschot, Müller, Sieber, Percin, van Oudheusden and Oberleithner2020) presented numerical simulations and experimental results of vortex breakdown states.
The study of the instability of swirling flow has a long history. Previous stability theories of vortex flow (Kelvin Reference Kelvin1880; Rayleigh Reference Rayleigh1917; Synge Reference Synge1938; Howard & Gupta Reference Howard and Gupta1962; Leibovich & Stewartson Reference Leibovich and Stewartson1983) focus on the tendency of a small perturbation imposed on a base swirling flow in an infinitely long pipe. The perturbation is in the form of Fourier modes both in the azimuthal and axial directions which can be analysed using the normal mode method. Rayleigh (Reference Rayleigh1917) stated a linear stability criterion that an axisymmetric base rotating flow in an infinitely long pipe is neutrally stable with respect to axisymmetric perturbations if the base flow circulation function
$K$
satisfies
$({1}/{r^3}) ({\textrm{d}K^2}/{\textrm{d}r})\gt 0$
. This criterion is also valid for a rotating flow with a uniform axial velocity. Synge (Reference Synge1938) pointed out that Rayleigh’s criterion, which was initially developed for inviscid flow, can also apply to viscous flow. Later Howard & Gupta (Reference Howard and Gupta1962) studied a more general base flow with axial velocity profile
$W(r)$
and showed that it was stable to small axisymmetric perturbations if
$({1}/{r^3})({\textrm{d}K^2}/{\textrm{d}r})-( {1}/{4})({\partial W}/{\partial r})^2\gt 0$
. However, these criteria are valid only when a swirling flow is in an infinitely long pipe or a finite-length pipe with periodic boundary conditions at the inlet and outlet. Experimental results, for example by Harvey (Reference Harvey1962) showed instability and bubble breakdown phenomena under approximately axisymmetric conditions, contrary to classical stability theory according to which the base flow is stable to axisymmetric perturbations. Although there are many phenomena in swirling flows that are inconsistent with the classical stability theory, the limitations of the classical stability theory have often been overlooked.
Wang & Rusak (Reference Wang and Rusak1996a ) developed an inviscid flow stability theory, which eliminated the crucial limitation of the classical stability theory. They found that inlet and outlet boundary conditions dramatically altered the stability of the base swirling flow in a finite-length pipe. By considering solid body rotation as the base flow and solving the linearized perturbation equations under a set of non-periodic inlet and outlet conditions, Wang & Rusak (Reference Wang and Rusak1996a ) obtained an axisymmetric non-Fourier mode of perturbation that became unstable when the swirl ratio of the base flow was above a critical level. Wang & Rusak (Reference Wang and Rusak1996b ) pointed out that this kind of instability was directly related to the physical mechanism leading to the axisymmetric breakdown state. The intrinsic assumption of an infinitely long pipe or periodic boundary condition is the reason for Rayleigh’s criterion and classical works of a similar type to fail to predict the instability that has been observed in practice.
Recently, Wang et al. (Reference Wang, Rusak, Gong and Liu2016) and Gong (Reference Gong2017) extended the inviscid stability analysis of the axisymmetric (
$m = 0$
) modes to viscous three-dimensional instability modes for a general swirling flow in a finite-length pipe. Their work showed that a series of both axisymmetric and spiral instability modes appears in sequence as the swirl ratio increases above certain critical levels. The stability results of Wang et al. (Reference Wang, Rusak, Gong and Liu2016) predict a stability boundary in the Reynolds number versus swirl ratio diagram. The diagram indicates that with the imposed inlet and outlet boundary conditions, the spiral type of perturbation becomes unstable at swirl levels below the critical swirl ratio for the onset of the axisymmetric type of instability. Feng et al. (Reference Feng, Liu, Rusak and Wang2017a
,
Reference Feng, Liu, Rusak and Wangb
, Reference Feng, Liu, Rusak and Wang2018) studied three-dimensional viscous flow dynamics of perturbations on a solid-body rotation via direct numerical simulations (DNS). Their simulation results of the neutral stability line and the global dynamics of the perturbed solid-body rotation are in agreement with the stability analysis of Wang et al. (Reference Wang, Rusak, Gong and Liu2016).
The solid-body rotation vortex profile does not have a finite-size vortex core. While it is of both theoretical and practical importance, most vortices in real flows consist of a small vortex core surrounded by a nearly inviscid potential vortex. The Lamb–Oseen vortex is an analytical model for such practical vortices. Rayleigh criterion yields the conclusion that vortices with the Lamb–Oseen vortex velocity profile are stable to axisymmetric perturbations, which is contrary to experimental observations that the Lamb–Oseen vortex evolves to axisymmetric vortex breakdown when the swirl exceeds a critical swirl ratio (Billant et al. Reference Billant, Chomaz and Huerre1998). Wang & Rusak (1997a) and Rusak, Wang & Whiting (Reference Rusak, Wang and Whiting1998) studied the Lamb–Oseen vortex with uniform axial velocity. They analysed the axisymmetric Euler equations and established the evolution path from the columnar base flow to the onset of vortex breakdown. Wang (Reference Wang2008) discovered the important effect of energy transfer occurring in the vortex core and at the pipe outlet boundary. Wang & Rusak (1997b) also studied the stability limit and the effect of slight viscosity when the swirl ratio is at a near-critical level. No DNS based on the full Navier–Stokes equations has been performed yet for such flows.
Weakly non-parallel flow convective instability and absolute instability analyses were used to explain the vortex breakdown phenomena in swirling flow (see Ruith et al. Reference Ruith, Chen, Meiburg and Maxworthy2003; Gallaire et al. Reference Gallaire, Ruith, Meiburg, Chomaz and Huerre2006). They focused on the vortex breakdown in a semi-infinite domain and found that the bubble breakdown state was unstable with a spiral instability mode and evolved to a spiral breakdown state. Moreover, they claimed that the onset of the spiral vortex breakdown was initiated by a nonlinear global mode due to the absolute instability of the bubble vortex breakdown. However, the axisymmetric bubble breakdown state in Ruith et al. (Reference Ruith, Chen, Meiburg and Maxworthy2003) did not satisfy the weakly non-parallel assumption. Later, Meliga & Gallaire (Reference Meliga and Gallaire2011) relaxed the locally parallel hypothesis and imposed a realistic, general perturbation mode similar to the analysis of Wang & Rusak (Reference Wang and Rusak1996a ). Meliga & Gallaire (Reference Meliga and Gallaire2011) verified the accuracy of their stability analysis by comparing with DNS results in Ruith et al. (Reference Ruith, Chen, Meiburg and Maxworthy2003). Neither the classical stability theory nor the convective instability–absolute instability analysis appear to adequately explain instability modes in a finite domain with strong axial inhomogeneity under the non-periodic inlet–outlet boundary conditions.
In this paper, we study the three-dimensional, incompressible and viscous flow dynamics of perturbations on a base Lamb–Oseen vortex flow with a uniform axial velocity entering a rotating, finite-length, straight circular pipe. The present work focuses on identifying the dynamics of three-dimensional instability of the base swirling flow. We identify different linear perturbation modes in the perturbed Lamb–Oseen vortex and their full linear–nonlinear evolution paths to two types of final stable states: a steady non-columnar but symmetric flow state and a spiral breakdown state. Linear growth rates of disturbances of both axisymmetric and spiral modes and the corresponding mode shapes are extracted from data obtained by DNS. Effects of swirl ratio and Reynolds number on the instability modes are investigated. The complete stability boundaries of the perturbation modes are determined in the Reynolds number and swirl ratio
$(\textit{Re},\omega )$
plane. We also reveal the detailed kinetic energy transfer mechanisms responsible for the instabilities by an energy method based on the Reynolds–Orr equation.
The mathematical model and numerical method are discussed in § 2. Section 3 presents computational results and analysis of the linear and nonlinear evolution processes. Two typical linear instability modes that appear in different regimes of the
$(\textit{Re},\omega )$
plane are discovered. Their linear growth rates, nonlinear development and then long-time behaviour until the flow reaches their final states are presented. The stability boundaries corresponding to these two types of instability are mapped out in the
$(\textit{Re},\omega )$
plane. Section 4 presents analysis of the energy transfer mechanisms of the two instability modes. Finally, § 5 draws conclusions from the present study.
2. Physical model and numerical method
2.1. Governing equations
The incompressible three-dimensional Navier–Stokes equations are solved in the cylindrical coordinate system
$(r,\theta ,x)$
. Let the three components of
$(u_r,u_\theta ,u_x)$
represent the radial, azimuthal and axial velocity components, respectively. A new set of velocity variables
\begin{equation} (q_r,q_\theta , q_x)=(ru_r, ru_\theta ,u_x) \end{equation}
is introduced to avoid the singularities at the axis of symmetry in numerical computations. Notice the quantities
$q_i$
do not all have the same physical dimensions as velocity. In the following presentations, dimensional variables are marked with asterisks. The flow is assumed to have constant density
$\rho ^*$
and viscosity
$\mu ^*$
. Flow quantities, space and time coordinates are non-dimensionalized using the pipe radius
$R^*$
as the length scale and the axial velocity
$U^*$
as the velocity scale. Time is scaled by
$R^*/U^*$
and the pressure is scaled by
$\rho ^* U^{*2}$
. The Reynolds number
$\textit{Re}={\rho ^*}U^*R^*/\mu ^*$
is based on the pipe radius. The dimensionless Navier–Stokes equations in terms of the new variables are
\begin{align}&\qquad\qquad\qquad\qquad\qquad \frac {\partial q_r}{\partial r}+\frac {1}{r}\frac {\partial q_{\theta }}{\partial \theta }+r\frac {\partial q_x}{\partial x}=0, \\[-12pt] \nonumber \end{align}
\begin{align}& \frac {\textrm{D}q_r}{\textrm{D}t}=-r\frac {\partial p}{\partial r}+\frac {1}{\textit{Re}}\left [ \frac {\partial }{\partial r}\left ( r\frac {\partial }{\partial r}\left ( \frac {q_r}{r} \right ) \right ) \right . -\frac {q_r}{r^2}+\left . \frac {1}{r^2}\frac {\partial ^2q_r}{\partial \theta ^2}+\frac {\partial ^2q_r}{\partial x^2}-\frac {2}{r^2}\frac {\partial q_{\theta }}{\partial \theta } \right ]\!, \\[-12pt] \nonumber \end{align}
\begin{align}&\, \frac {\textrm{D}q_{\theta }}{\textrm{D}t}=-\frac {\partial p}{\partial \theta }+\frac {1}{\textit{Re}}\left [ \frac {\partial }{\partial r}\left ( r\frac {\partial }{\partial r}\left ( \frac {q_{\theta }}{r} \right ) \right ) -\frac {q_{\theta }}{r^2}+ \frac {1}{r^2}\frac {\partial ^2q_{\theta }}{\partial \theta ^2}+\frac {\partial ^2q_{\theta }}{\partial x^2}+\frac {2}{r^2}\frac {\partial q_{r}}{\partial \theta } \right ]\!, \\[-12pt] \nonumber \end{align}
\begin{align}&\qquad\qquad\,\, \frac {\textrm{D}q_x}{\textrm{D}t}=-\frac {\partial p}{\partial x}+\frac {1}{\textit{Re}}\left [ \frac {1}{r}\frac {\partial }{\partial r}\left ( r\frac {\partial q_x}{\partial r} \right ) +\frac {1}{r^2}\frac {\partial ^2q_x}{\partial \theta ^2}+\frac {\partial ^2q_x}{\partial x^2} \right ]\!. \end{align}
For enhanced computational accuracy, the substantial derivatives in the above equations are written in conservative form as
\begin{align}& \frac {\textrm{D}q_r}{\textrm{D}t}=\frac {\partial q_r}{\partial t}+\frac {\partial }{\partial r}\left ( \frac {q_{r}^{2}}{r} \right ) +\frac {\partial }{\partial \theta }\left ( \frac {q_rq_{\theta }}{r^2} \right ) +\frac {\partial q_rq_x}{\partial x}-\frac {q_{\theta }^{2}}{r^2}, \\[-12pt] \nonumber \end{align}
\begin{align}& \frac {\textrm{D}q_{\theta }}{\textrm{D}t}=\frac {\partial q_{\theta }}{\partial t}+\frac {\partial }{\partial r}\left ( \frac {q_{\theta }q_{r}^{}}{r} \right ) +\frac {1}{r^2} \frac {\partial q_{\theta }^{2}}{\partial \theta } +\frac {\partial q_{\theta }q_x}{\partial x}+\frac {q_{\theta }^{}q_r}{r^2}, \\[-12pt] \nonumber \end{align}
\begin{align}&\qquad\,\, \frac {\textrm{D}q_x}{\textrm{D}t}=\frac {\partial q_x}{\partial t}+\frac {1}{r}\frac {\partial q_rq_x}{\partial r}+\frac {1}{r^2}\frac {\partial q_{\theta }q_x}{\partial \theta }+\frac {\partial q_{x}^{2}}{\partial x}. \end{align}
2.2. The base flow
The dimensionless velocity components of the Lamb–Oseen vortex with uniform axial velocity and pressure distribution are
\begin{equation} { u_r=0, \, u_\theta =\frac {\omega }{4r}\big[1-\exp\big(-br^2\big)\big], \, u_x=1,\, p=\int _0^r \frac {1}{r}u_\theta ^2\ \textrm{d}r, } \end{equation}
where
$\omega$
and
$b$
are two dimensionless parameters defined as
\begin{equation} { \omega =\frac {\varGamma _0^*}{2{\pi }U^*R^*}, \,b=\frac {R^{*2}}{4\nu ^* t^*}. } \end{equation}
Here
$\varGamma _0^*$
and
$\nu ^*$
are the circulation of the vortex core and the kinematic viscosity, respectively. Notice
$\varGamma _0^*/2{\pi }R^*$
is the swirl velocity of the vortex at the outer edge
$r=1$
. Therefore, the
$\omega$
defined above is the ratio of the swirl velocity to the axial velocity, which measures the swirl strength of the Lamb–Oseen vortex.
The Lamb–Oseen vortex defined by (2.9) has a uniform axial velocity and zero radial velocity. Unlike the solid-body rotation vortex, it consists of an inner vortex core surrounded by a potential flow with a fixed total circulation
$\varGamma ^{*}_{0}$
. It is an idealization that captures the critical features of a real-world vortex, such as a vortex in a combustion chamber or the trailing vortex from an aeroplane wing tip, all of which have a tight finite vortex core size. The Lamb–Oseen vortex is one of those few exact analytical solutions to the full Navier–Stokes equations. However, it is unsteady. Its vortex core expands with time due to diffusion (Wu, Ma & Zhou Reference Wu, Ma and Zhou2015). If the vortex-core radius is defined as the location where the azimuthal velocity is maximum, its non-dimensional radius is then
\begin{equation} { r_0=2.2418 \sqrt {\frac {\nu ^* t^*}{R^{*2}}}=\sqrt {\frac {1.2564}{b}}. } \end{equation}
The parameter
$b$
can be used as a control parameter that specifies the size of the vortex core to match that of a real-world vortex to be studied. Garg & Leibovich (Reference Garg and Leibovich1979) confirms that experimental data of time-averaged swirl velocity profiles are well fitted by the analytical Lamb–Oseen vortex profile ‘frozen’ at chosen fixed dimensionless time
$b$
. In the present study,
$b$
is set as
$b=4$
, representing a medium core-size
$r_0=0.5604$
.
Notice that the Lamb–Oseen vortex frozen in time does not exactly satisfy the steady Navier–Stokes equations because of the slow nevertheless unsteady diffusion of the vortex core. The ratio of the diffusion time scale,
$t_{\nu }^*=R^{*2}/\nu ^*$
, to the convection time scale,
$t_{U}^*=R^*/U^*$
, is
\begin{equation} { \frac {t_{U}^*}{t_{\nu }^*}=\frac {R^*/U^*}{R^{*2}/{\nu }^*}=\frac {\nu ^*}{U^*R^*}=\frac {1}{\textit{Re}}. } \end{equation}
We assert that the vortex core expands slowly compared with the axial flow velocity for large Reynolds numbers. Nevertheless, after a sufficient amount of time, which may be short for low Reynolds numbers, the vortex core will diffuse to infinite size and the magnitude of vorticity dissipates to zero, making it useless as a base flow to represent the state of a real-world stationary vortex with fixed core size and strength.
To enable a meaningful stability study, a small constant body-force term
$\boldsymbol{f}=\omega b^2( r^2e^{-br^2})/\textit{Re}$
is added to the right-hand side of (2.4) in our computations. The constant body-force term is designed to make the columnar Lamb–Oseen vortex a steady base flow that captures the major features of a real-world vortex. Its use is to represent the effect of small local radial flow velocities in a real-world vortex that counters viscous diffusion to keep the vortex tight. It makes the idealized columnar Lamb–Oseen vortex an exact steady solution of the Navier–Stokes equations. In the terminology of dynamical systems, this makes the steady base flow a fixed point. A detailed discussion on the relevance and effectiveness of the source term is provided in Appendix A. There, we will provide evidence to justify the use of this source term for stability studies.
2.3. Boundary conditions
The inlet boundary conditions at
$x=0$
are set to be the base Lamb–Oseen vortex with uniform axial velocity expressed in terms of
$(q_r,q_\theta , q_x)$
as
\begin{equation} q_r=0, q_{\theta }=\frac {\omega }{4}\big [ 1-\exp\big ({-}br^2 \big) \big], q_x=1. \end{equation}
At the outlet
$x = L$
the following convective boundary conditions are specified:
\begin{equation} \left \{ \begin{aligned} \frac {\partial q_r\left ( r,\theta ,L,t \right )}{\partial t}+C\frac {\partial q_r\left ( r,\theta ,L,t \right )}{\partial x}=0, \\ \frac {\partial q_{\theta }\left ( r,\theta ,L,t \right )}{\partial t}+C\frac {\partial q_{\theta }\left ( r,\theta ,L,t \right )}{\partial x}=0, \\ \frac {\partial q_x\left ( r,\theta ,L,t \right )}{\partial t}+C\frac {\partial q_x\left ( r,\theta ,L,t \right )}{\partial x}=0, \end{aligned} \right . \end{equation}
where
$C$
is a constant value, representing the advection speed of the mean flow. The magnitude of
$C$
is not critical to the solution (Feng et al. Reference Feng, Liu, Rusak and Wang2017a
). In our simulations
$C$
is set to 1. This convective boundary condition, also known as non-reflective boundary condition, can effectively avoid non-physical feedback at the outlet and is successfully used in simulations of axisymmetric coaxial jets (Salvetti, Orlandi & Verzicco Reference Salvetti, Orlandi and Verzicco1996), rotating circular flow (Ortega-Casanova & Fernandez-Feria Reference Ortega-Casanova and Fernandez-Feria2008) and three-dimensional pipe flow (Ruith et al. Reference Ruith, Chen, Meiburg and Maxworthy2003). Slip conditions are specified at
$r=1$
which means that a Lamb–Oseen vortex confined in an inviscid round pipe is simulated.
Classical vortex stability theory assumes homogeneous or periodic boundary conditions along the vortex. However, any vortex found in reality is confined in a finite domain. Therefore, it is important to apply proper boundary conditions that will closely represent the real physical conditions and minimize any non-physical reflections in numerical computations at the same time. In an experimental set-up and also most real flow situations, the flow conditions upstream are specified, while at the exit the flow is allowed to go downstream without immediate obstruction or disturbance. The choice of the Dirichlet boundary conditions at the inlet and convective outlet boundary conditions at the exit described above is designed to best represent the physical situation. It is this set of non-periodic boundary conditions that gives rise to the possibility of instability of the Lamb–Oseen vortex in a finite-length pipe.
2.4. Stability analysis by DNS
Let (
$u_{r,0}$
,
$u_{\theta ,0}$
,
$u_{x,0}$
,
$p_0$
) denote the base Lamb–Oseen vortex velocity and pressure. Disturbances of velocity and pressure in the DNS stability analysis are defined as the difference between the computed instantaneous velocity and pressure at time
$t$
and (
$u_{r,0}$
,
$u_{\theta ,0}$
,
$u_{x,0}$
,
$p_0$
),
\begin{equation} \left \{ \begin{aligned} u_r^{\prime }\quad&=u_r(t) -u_{r,0},\\u_{\theta }^{\prime }\quad&=u_{\theta }(t) -u_{\theta ,0},\\ u_x^{\prime }\quad&=u_x(t) -u_{x,0},\\ p^{\prime }\quad&=p(t) -p_{0}. \end{aligned} \right . \end{equation}
The stability analysis examines the evolution of the above perturbation velocity components starting with an arbitrary initial field under specified well-posed boundary conditions.
In a theoretical linear stability analysis one may assume the following form of perturbation solutions as Wang et al. (Reference Wang, Rusak, Gong and Liu2016) did:
\begin{align} &u_{x}^{\prime }(x,r,\theta ,t)=\textrm{Re}(\mathscr{U}_m(x,r)e^{\textit{im}\theta +\sigma t}), &&u_{r}^{\prime }(x,r,\theta ,t)=\textrm{Re}(\mathscr{V}_m(x,r)e^{\textit{im}\theta +\sigma t}), \notag\\ &u_{\theta }^{\prime }(x,r,\theta ,t)=\textrm{Re}(\mathscr{W}_m(x,r)e^{\textit{im}\theta +\sigma t}), &&p^{\prime }(x,r,\theta ,t)=\textrm{Re}(\mathscr{P}_m(x,r)e^{\textit{im}\theta +\sigma t}), \end{align}
where
$\sigma = \sigma _r +i \sigma _i$
is the linear growth rate of the perturbations. Its real part
$\sigma _r$
represents the growth rate of the velocity amplitude and its imaginary part gives the angular frequency of the possible oscillatory behaviour of the perturbation modes. The above form of solution represents one single azimuthal component of the full perturbation solution, where
$m=0,1,2,\ldots $
is the azimuthal wave number and
$\mathscr{U}_m, \mathscr{V}_m, \mathscr{W}_m$
and
$\mathscr{P}_m$
are the corresponding mode shapes for the
$m$
th mode. The axisymmetric mode is given by the
$m=0$
mode. Any
$m \gt 0$
mode yields a non-axisymmetric solution. In a linear stability analysis, the growth rate
$\sigma$
for each perturbation mode may be determined by solving an eigenvalue problem independent of other modes since the individual modes do not interact within the linear assumption for small amplitude perturbations.
In the present numerical study, we add to the base stationary flow a random initial perturbation velocity field and perform direct numerical computation to obtain the full velocity field. Perturbations at every time instant are computed by using (2.15). The resulting perturbation field consists of all possible ‘modes’ and frequencies. Separation of spatial modes and constant growth rates may not always be possible. However, for the present problem, we can always decompose our full DNS perturbation solutions into an infinite number of Fourier modes in the azimuthal direction due to the retained azimuthal homogeneity. For example, at any given time, the axial perturbation velocity function
$u_x^{\prime }(x,r,\theta ,t)$
can be decomposed as
\begin{align} u_{x}^{\prime }(x,r,\theta ,t) = \sum _{-\infty }^{+\infty }{g_m(x,r,t) e^{\textit{im}\theta }} \end{align}
where
\begin{align} g_m(x,r,t)=\frac {1}{2\pi }\int _0^{2\pi }{u_{x}^\prime \left (x,r, \theta , t \right ) e^{-im\theta }\textrm{d}\theta }. \end{align}
We might call
$g_m(x,r,t)$
the
$m$
th azimuthal ‘mode’ of
$u_x^{\prime }$
, recognizing that its spatial structure may change with time. This is a basic method used in this article for extracting the linear stability modes from DNS results.
In the linear growth range, it will be demonstrated that
$g_m(x,r,t)$
can be expressed in the following space–time separated form:
\begin{equation} g_m(x,r,t) = g_m^{\textit{norm}}(x,r) e^{\sigma _{m} t} \end{equation}
where
$g_m^{\textit{norm}}(x,r)$
is the true spatial mode shape for the
$m$
th mode, consistent with the assumptions by Wang et al. (Reference Wang, Rusak, Gong and Liu2016) in (2.16). For real functions of
$u_x^\prime (x,r,\theta , t)$
, we need only examine the
$m= 0, 1, 2, \ldots $
modes. The corresponding complex growth rate for each mode is then
\begin{equation} \sigma _m = \sigma _{m,r} +i \sigma _{m,i}. \end{equation}
The real part
$\sigma _{m,r}$
stands for the amplitude decay or growth of the
$m$
th azimuthal perturbation mode. Any non-zero imaginary part
$\sigma _{m,i}$
signifies a travelling wave in the perturbations.
In the DNS stability analysis, one can easily calculate the time history of the kinetic energy of the
$m$
th azimuthal mode without assuming the above spatial and time separation by using
\begin{equation} E_{m} \left ( t \right ) =\frac {1}{2}\int _V{\left [g_m(x,r,t)\right ]^2\textrm{d}V}, \end{equation}
where
$V$
is the domain of the flow field. A local instantaneous growth rate for this energy can be defined as
$({\textrm{d}E_m}/{\textrm{d}t})/E_{m}$
. For consistency with the growth rate definition for the velocity magnitude, we here define a local instantaneous growth rate by using
$\sqrt {E_m}$
instead of
$E_m$
. Thus, we define
\begin{equation} \sigma _{m,r}= \frac {\textrm{d}(\ln {\sqrt {E_m}})}{\textrm{d}t}. \end{equation}
A negative
$\sigma _{m,r}$
indicates decay of the perturbation mode. All modes must have negative
$\sigma _{m,r}$
for stability. A positive
$\sigma _m$
value for any
$m$
proves instability.
The growth rate calculated by the above equation will be a constant if the evolution of the perturbation velocity of mode
$m$
is in the linear range and the value of this constant should be equal to that defined in (2.19) and (2.20). Therefore, it is convenient to extract the linear growth rate for each azimuthal mode from the DNS data by using (2.22).
By Parseval’s equality, we can also obtain the total kinetic energy
\begin{equation} E \left ( t \right ) =\frac {1}{2}\int _V{u_{x}^{\prime }}^2\textrm{d}V = \sum _{-\infty }^{+\infty } E_{m}. \end{equation}
A local instantaneous growth rate based on this total energy may also be calculated:
\begin{equation} \sigma _r=\frac {\textrm{d}}{\textrm{d}t} \ln \sqrt {E} = \frac {1}{2E} \frac {\textrm{d}E}{\textrm{d}t}. \end{equation}
The total perturbation energy will be dominated by the mode that has the highest positive
$\sigma _{m,r}$
as time increases. Very often, one single mode becomes unstable. In that case, the overall growth rate based on the total energy from (2.24) is practically the same as the
$\sigma _{m,r}$
of the unstable mode. (Similar analysis can be applied to the total kinetic energy of the perturbation including all three velocity components:
$ E ( t ) =({1}/{2})\int _V{ ( { u^{\prime }_{\theta }}^2+{u^{\prime }_r}^2+{u^{\prime }_x}^2 ) \textrm{d}V}.$
).
To extract the imaginary part of the growth rate
$\sigma _{m,i}$
in (2.20), however, one may need to consider phase shift of the mode shapes with time or perform Fourier analysis of the DNS data in time to determine the frequency of the perturbations. We will demonstrate the use of such methods in §§ 3.1.1 and 3.2.1. It suffices here to note that
\begin{align} \sigma _{m,i} = 2\pi f_{m} = 2\pi /T_{m} \end{align}
where
$f_{m}$
and
$T_{m}$
are the frequency and period of mode
$m$
should it behave in travelling wave nature.
2.5. Numerical scheme, grid convergence and parallel computation
The continuity and Navier–Stokes equations (2.2) to (2.5) are solved by the fractional step method developed by Kim & Moin (Reference Kim and Moin1985). The convective and viscous terms are advanced temporally by the explicit third-order Runge–Kutta scheme and the implicit Crank–Nicolson scheme, respectively. Spatial derivatives are discretized by second-order central differences on staggered grids. Velocities are located at the centre of the cell face and pressure is at the centre of the cell. In each time step, there are two computational processes. First, an intermediate non-solenoidal velocity field is introduced and the Navier–Stokes equations are advanced using the pressure at the previous time step. This intermediate velocity field is used in the calculation of an intermediate pressure field. In the second stage, the intermediate fields have been obtained thus the solenoidal velocity field and the pressure field can be calculated by solving the Poisson equation. The time step is determined according to the Courant–Friedrichs–Lewy (
$\textit{CFL}$
) condition of the third-order Runge–Kutta method. The theoretical stability criterion is
$\textit{CFL}\lt \sqrt {3}$
. We choose
$ \textit{CFL}_{max}=1.5$
to limit the time step. The reader is referred to Verzicco & Orlandi (Reference Verzicco and Orlandi1996) and Ruith et al. (Reference Ruith, Chen, Meiburg and Maxworthy2003) for more details on the specific method to solve the above equations.
Under the initial and boundary conditions described above, the Lamb–Oseen vortex flow is a stationary solution to the full three-dimensional unsteady Navier–Stokes equations. Care must be taken to ensure this is true also in the discrete form of the numerical computation. In the present study, we make sure that the base Lamb–Oseen vortex velocity field with the specified boundary conditions described in § 2.2 is divergence free and the time derivatives of velocity components evaluated by computing the right-hand sides of (2.3)–(2.5) are all within machine zero.
Although the Navier–Stokes equations are cast in the cylindrical coordinate system in (2.2)–(2.5) and the initial Lamb–Oseen vortex is axisymmetric, no axisymmetry is assumed in the computation.
A small random initial perturbation velocity field is then added to the base flow to study its time evolution for stability analysis. In principle, the added velocity perturbations must be divergence free to ensure it is an allowed solution to the Navier–Stokes equations. In practice, the base flow can be perturbed by a random three-dimensional disturbance to one or all velocity components. In the present study, only the axial velocity component is perturbed by a three-dimensional random function. The computational algorithm corrects the velocity field to be divergence free within the first few iterations by solving the Poisson equation for pressure. The initial random perturbations in the axial velocity induce perturbations in all other directions.
We present here the case of dimensionless pipe length
$L=6$
. The choice of pipe length and its impact are investigated in § 5. Meshes are uniform in the axial, azimuthal and radial directions. The accuracy of the flow calculations is validated by a grid refinement study. The computational domain (
$0 \leqslant x\leqslant L= 6$
,
$0 \leqslant r \leqslant R=1$
,
$0\leqslant \theta \leqslant 2\pi$
) is resolved on three successively refined grids with dimensions
$128\times 64\times 64$
,
$256\times 128\times 128$
and
$512\times 256\times 256$
. At each time step, solution of the Poisson equation for pressure ensures that the divergence of the velocity field is within
$10^{-13}$
. Table 1 shows the computed linear growth rate of the perturbation kinetic energy on those three grids. The values computed on the
$256\times 128\times 128$
and
$512\times 256\times 256$
grids show little difference (only
$1.1\times 10^{-3}$
relative error) and the converged growth rate also agrees with the theoretical value by Gong (Reference Gong2017) with less than
$1.4\times 10^{-2}$
relative error. For consideration of computational cost and accuracy, the following results are obtained on the medium
$256\times 128\times 128$
grid.
Table 1. The linear growth rates computed by different grids and the theoretical value at
$\textit{Re}=5000$
and
$\omega =3.8$
.
Because the three-dimensional DNS requires much central processing unit (CPU) time and computer memory, the numerical simulation is accelerated by parallelizing the code with a message-passing interface. The cost of time is reduced greatly by parallelized simulations. For example, the time needed to compute the flow field is approximately 18 days with the serial version of the code but only approximately 10 hr with the parallel code on 112 CPU cores for a typical three-dimensional case of
$\textit{Re}=5000$
and
$\omega =3.8$
on a
$256\times 128\times 128$
grid. Three meshes are used to test the performance of our parallel programme, see figure 1. It illustrates that the speedup and efficiency (the slope of speedup lines) increase with mesh size. Although the speedup achieved by our code is lower than the ideal linear speedup represented by the black line, the speedup scales almost linearly with the number of cores on sufficiently fine grids.
Figure 1. Parallel speedup factor versus number of used computer cores.
3. Computational results and analysis
We present in this section computational results and detailed analysis of two typical cases of complete linear and nonlinear evolution of the perturbations to the Lamb–Oseen vortex in a finite-length pipe using the methodology described in § 2.
Figure 2. Time history of growth rate and perturbation energy for the case
$\omega = 3.85$
and
$\textit{Re} = 5000$
.
Linear perturbation modes become unstable when the swirl ratio is beyond a certain critical swirl level at a given Reynolds number or when the Reynolds number is beyond a certain critical Reynolds number at a given swirl ratio. During the linear phase, different modes behave independently. Each mode may be easily identified with its own mode shape and growth rate. The mode that has the highest growth rate determines the linear stability boundary and dominates the development of instability in the linear stage. With further growth of perturbations, unstable modes lead to large-amplitude nonlinear modes and may evolve to different final equilibrium states. The identification of the different unstable modes provides a fundamental connection to the nonlinear growth of the instability modes and the eventual passage to a final equilibrium state. Two representative cases at
$\textit{Re}=5000, \omega =3.85$
and
$\textit{Re}=700, \omega =4.3$
are identified here. The first case is dominated by an axisymmetric instability mode (
$m=0$
) and it leads to a final steady axisymmetric but non-columnar flow. The second case is dominated by a spiral mode (
$m=1$
) and it leads to a final unsteady periodic spiral vortex breakdown state.
For each case, we present the velocity trajectory at a given point to depict the global nonlinear evolution towards the corresponding final equilibrium state. We examine the structure of the linear mode shapes and the final equilibrium states.
Finally, stability boundaries of different modes in the
$(\textit{Re}, \omega )$
plane are determined through the DNS simulations. These boundaries identify the most unstable mode and demarcate regions of stability of the Lamb–Oseen vortex in the
$(\textit{Re}, \omega )$
plane for both the axisymmetric (
$m=0$
) and the spiral (
$m=1$
) modes.
3.1. Case 1: flow dominated by instability of the axisymmetric
$m=0$
mode
We start by examining the flow at
$\textit{Re}=5000$
. An initial low-amplitude random velocity perturbation field is added to the base Lamb–Oseen vortex flow. As will be shown later in § 3.1.3, the initial perturbations of the
$m=0$
component increase when the swirl ratio
$\omega$
is above a critical value
$\omega = 3.762$
, indicating that the base Lamb–Oseen vortex flow becomes unstable at this Reynolds number and swirl ratio condition. For the flow with
$3.762 \lt \omega \lt 3.91$
, the
$m=0$
perturbation mode grows quickly with time and we examine such a case with
$\textit{Re}=5000, \omega =3.85$
in this range.
3.1.1. Linear growth of the perturbed flow
Figure 2 shows the time history of the square root of the perturbation energy (
$E^{0.5}$
) and its instantaneous growth rate based on our DNS simulations. The growth rate is found to be negative at the beginning and reaches zero at
$t = 30$
. The initial random disturbances imposed on the base flow contain energies of a wide range of modes. All stable modes of the perturbation decay quickly except the most unstable mode, whose growth overtakes the decay of other initial perturbation components. At around
$t = 400$
the growth rate
$\sigma$
approaches a near-constant value of
$0.0133$
and stays at this value for 500 non-dimensional time units. This period is the linear growth stage of the perturbation. The computed growth rate
$\sigma =0.0133$
agrees closely with the value (
$\sigma =0.0131$
) predicted by Gong (Reference Gong2017) through linear analysis. As the disturbance grows larger, the flow enters into a nonlinear growth stage, for which any linear analysis including that by Gong (Reference Gong2017) does not apply. The present full three-dimensional Navier–Stokes solution, however, has no such limitation. Figure 2 shows that the linear stage ends at
$t=900$
with
$\sigma$
decreasing quickly to near zero at
$t = 1250$
. From then on, the total energy of the perturbation appears to saturate with minor fluctuations until it reaches a steady level beyond
$t=2400$
. Indeed, the flow has reached a new stable steady state of non-columnar accelerated flow shown in figure 3 at
$t=6000$
. The flow in the
$\theta -r$
plane is axisymmetric, as shown by the velocity and pressure field in figure 3(b). Low-pressure regions correspond to accelerated axial velocities.
We see from this case the process of the initial unstable axisymmetric columnar Lamb–Oseen vortex flow evolves to another steady axisymmetric state due to instability of an axisymmetric perturbation mode. Rusak et al. (Reference Rusak, Wang, Xu and Taylor2012) predicted a similar phenomenon of the flow starting from a columnar base flow to a non-columnar accelerated flow type by using an inviscid weakly nonlinear analysis method. Their solution has the same qualitative features as the present solution despite the use of only an inviscid flow model and with slightly different boundary conditions. Recently, Qiao et al. (Reference Qiao, Shi, Wang and Liu2025) revealed that the final state of this flow is one solution belonging to a branch of steady, axisymmetric non-unique solutions of the Navier–Stokes equations for the Lamb–Oseen vortex base flow in a finite-length pipe.
Figure 3. Projected perturbation velocity vectors and contours of perturbation pressure for the case
$\textit{Re}= 5000$
and
$\omega =3.85$
at
$t = 6000$
: (a) in the
$x{-}z$
plane, (b) in the
$\theta -r$
plane at
$x=3$
.
In order to identify and examine the particular modes of instability, we decompose the perturbation velocity field by using the method described in § 2.4. Figure 4 presents the time history of the energies of mode
$m = 0$
and
$1$
of our DNS results for this case. The perturbation energy
$E$
and growth rate
$\sigma$
for the
$m=0$
and
$1$
modes are denoted by the subscripts
$0$
and
$1$
, respectively, while those for the total perturbation are denoted by the subscript DNS. The energy of the perturbation mode
$m = 1$
decays rapidly with a negative linear growth rate
$\sigma _{1,r}=-0.0039$
until it reaches and then remains at a machine-zero level of
$10^{-13}$
. We thus assert that mode
$m=1$
is stable to three-dimensional disturbances and does not have any impact on the subsequent flow evolution. In addition, we find that all higher azimuthal modes decay faster than mode
$m=1$
. For example, modes
$m=2$
and
$m=3$
decrease to a magnitude of
$10^{-15}$
at
$t=240$
and remain at this value thereafter. Consequently, it is not surprising to see in figure 4 that the growth history for mode
$m = 0$
overlaps fully with that for the total perturbation field except in the very short early stage (
$t\lt 30$
) when some of the initial stable perturbation modes have not died out yet. This analysis shows that the axisymmetric mode
$m = 0$
is the only one that dominates the flow dynamics during the linear growth stage. When the amplitude of this mode grows to a sufficiently high level the evolution goes into the nonlinear stage.
Figure 4. Time history of growth rate and energy of perturbations of mode
$m = 0$
and mode
$m = 1$
for the case
$\textit{Re}=5000$
and
$\omega =3.85$
.
Figure 5. Normalized profiles of perturbation axial velocity for mode 0 at
$t=560$
,
$640$
and
$720$
: (a) along axial direction and (b) along radial direction.
Next, we examine the spatial structure of the unstable
$m=0$
mode. Figure 5 plots the axial perturbation velocity profiles of
$g_0^{\textit{norm}}$
at
$t=560$
,
$640$
and
$720$
in the linear growth stage. For each time instant, the axial profile is normalized by its value at the outlet while the radial profile by its value at the central axis. The spatial profiles of the perturbation velocity at the different time instants apparently collapse into a single line, proving that the dependence of the solution on
$t$
is separable from that on space in the linear growth stage in the form of (2.19), and the velocity profiles in figure 5 are essentially those of the
$g^{\textit{norm}}_0$
in (2.19). These results by solving the full nonlinear Navier–Stokes equations support the use of the method of separation of space and time variables and eigenvalue analysis by Gong (Reference Gong2017) in the linear range of the evolution of the perturbations.
Furthermore, figure 5(b) shows the radial profiles of
$g^{\textit{norm}}_0$
at different axial
$x$
locations almost collapse into one single line, while figure 5(a) shows only slightly dissimilar axial profiles for different radial locations. These findings based on DNS data support the assumed separation-of-variables solution
$\varPsi (X,r,T)=\phi _B(r)A(X,T)$
in the long-wave asymptotic analysis by Rusak, Granata & Wang (Reference Rusak, Granata and Wang2015). Here
$X$
and
$T$
are rescaled axial coordinate and time,
$\varPsi$
is the stream function and
$\phi _B$
and
$A$
are functions to be determined. Figure 5(b) shows that the spatial structure of the unstable
$m=0$
mode may be further separated approximately as standing waves in
$r$
and
$x$
, respectively.
Although the instability mode for this case is the axisymmetric
$m=0$
mode, it is insightful to examine the decaying
$m=1$
mode. The growth rate calculated by the kinetic energy of this mode yields a constant
$\sigma _{1,r} = -0.0039$
in the linear range as shown in figure 4. Figure 6 shows the azimuthal profiles of the extracted axial perturbation velocity
$u_1$
of mode
$m=1$
multiplied by
$e^{-\sigma _{1,r}t}$
at five consecutive time instants. These profiles are clearly footprints of one single azimuthal travelling wave moving in the positive
$\theta$
direction. The phase difference is
$\Delta \theta =1.2665$
in an interval
$\Delta t=1$
, which determines the phase speed as
$\sigma _{1,i}=-1.2665$
. Its sign is determined by the direction of the travelling wave. The perturbation flow field is seen to rotate anticlockwise around the pipe centreline as shown in figure 6. Therefore, we determine
$\sigma _{m=1}=-0.0039-1.2665i$
and reveal that the
$m=1$
mode is an unsteady spiral mode despite its decaying nature.
Figure 6. The phase difference and direction of travelling wave
$u_1$
in the azimuthal direction at
$x=3$
with time interval
$\Delta t=1$
.
Figure 7. Normalized profiles of perturbation axial velocity for mode
$m=1$
at
$t=480{-}483$
: (a) along axial direction and (b) along radial direction.
Figure 8. The perturbation velocity trajectory for
$\textit{Re}= 5000$
and
$\omega = 3.85$
: (a)
$t$
= 0–30, the initial transient stage; (b)
$t$
= 0–900, the linear growth stage; (c)
$t$
= 30–1250, the complete linear and nonlinear stage; (d)
$t$
= 1250–4000, final stage spiralling towards a stable fixed point; (e)
$t$
= 30–4000, the long-time behaviour of evolution.
Figure 7 shows the profiles of the axial perturbation velocity
$g_1$
at
$t=480{-}483$
in the linear growth stage. For each time instant, the axial profile of
$g_1$
is normalized by its value at the outlet while the radial profile by its value at the central axis. Once again, we see the normalized spatial shapes are fixed in time showing that the variable
$t$
is separable at selected linear growth stage following the form of (2.19):
$g_m(x,r,t) = g_m^{\textit{norm}}(x,r) e^{\sigma _{m} t}$
. In addition, it is observed that the axial profiles of
$g^{\textit{norm}}_1$
are almost fully similar for different
$r$
locations, but the radial profiles show significant variations for different
$x$
locations. Therefore, it is not proper in a theoretical analysis to assume solutions of the separation-of-variables type in space for the
$m= 1$
mode.
Although not shown here, the radial and azimuthal velocity components
$u_r^{\prime }$
and
$u_\theta ^{\prime }$
and the pressure perturbation field
$p^\prime$
all show the same characteristics.
The DNS simulations yield only full velocity fields. We have demonstrated how one can separate the perturbation field into its composing modes and identify their growth/decay rates, spatial structure, wave speeds (if any) by using the present method of analysis.
3.1.2. The global dynamics and the final state of the flow
At the beginning of the flow evolution, the perturbation is small and the growth of any unstable modes is controlled by a linear process. When the perturbation reaches a finite size, a conventional linear stability analysis becomes invalid while the present method based on DNS remains valid. The long-time behaviour and global nonlinear dynamics of the swirling flow for this case are discussed in this section.
Figure 8 is the phase portrait for the perturbation velocity at the spatial point
$x=3$
,
$r=0.25$
and
$\theta =0$
during several important time intervals for
$0\lt t\lt 4000$
. It is helpful to refer back to figure 4 in the discussion of the different stages of the evolution of perturbation velocity. In the initial stage from
$t = 0$
to
$t = 30$
, various perturbation modes except
$m = 0$
decay in time and the perturbation velocity of the flow field changes rapidly and irregularly as shown in figure 8(a). The decay of the stable modes, which comprise most energy of the initial random perturbations, dominates the evolution of the flow so that the velocity trajectory is drawn towards a point around (0, 0,
$10^{-5}$
). Figure 8(b) plots the trajectory of the perturbation velocity from
$t=0$
to
$t=900$
. The velocity change during
$30\lt t\lt 900$
(red line) is much larger than that in the time period
$0 \lt t \lt 30$
(black line, which is almost invisible in the figure). During this time period, the most unstable
$m = 0$
mode grows with a constant rate
$\sigma = 0.0133$
and dominates the flow dynamics. During the linear growth stage, the velocity trajectory moves along a straight line with constant ratios between the three perturbation velocities. Notice that the scales in all three velocity coordinate directions are made equal in the plots. The axial and azimuthal velocities grow at larger rates compared with the radial velocity, which indicates a tendency towards axial symmetry. The time period
$900 \lt t \lt 1250$
marks the end of the linear stage with a decreasing growth rate from 0.0133 to 0. The trajectory of the perturbation velocity moves fast along a slightly curved line due to nonlinearity as shown in figure 8(c). The axial and azimuthal velocities still grow significantly more than the radial component. In the time period
$1250 \lt t \lt 4000$
shown in figure 8(d), the trajectory moves along a spiral orbit that is eventually drawn into a fixed point far away from the initial base flow, signifying that the flow has reached another stable steady state (red square in figure 8
d). Finally, figure 8(e) plots the trajectory of the perturbation velocity evolving from near-zero initial perturbations to the base Lamb–Oseen flow to the final fixed point, an attractor of the dynamic system finite-distance away from the initial base flow.
Figure 9. Structure of final state of flow at
$t=4000$
for the case
$\textit{Re}=5000$
and
$\omega =3.85$
: (a) isosurface of vorticity
$|\varOmega |=0.7$
and contours of axial vorticity in the azimuthal plane; (b) contours of azimuthal vorticity in the
$x{-}z$
plane; (c) contours of axial vorticity in the
$x{-}z$
plane.
Figure 9 shows the isosurface of vorticity
$\left |{\varOmega }\right |=0.7$
as well as contours of azimuthal and axial vorticity components in the
$x{-}z$
plane of the perturbation field at
$t=4000$
. The light yellow isosurface represents the surface where vorticity magnitude
$|\varOmega | =0.7$
. Contours of the axial vorticity component in the
$y=0$
azimuthal plane are also displayed. Figures 9
(b) and 9(c) show that both the azimuthal and the axial vorticity components are axisymmetric. It is observed that the final state is a steady axisymmetric accelerated flow swirling around the axis in the same direction as the base flow.
It is to be noted that the structure of the final state of the flow shown in figure 9 bears close resemblance to that of the disturbance flow field of the linearly unstable
$m=0$
mode shown in figure 3. Despite the fact that the growth of the initial disturbance goes from the linear stage and through the nonlinear stage, the final flow state appears to have evolved from the dominant linearly unstable mode shape. In addition, computations demonstrate that the final flow field is stable to three-dimensional disturbances. Rusak et al. (Reference Rusak, Granata and Wang2015) studied the inviscid and axisymmetric Lamb–Oseen vortex by a long-wave asymptotic analysis. They also find a branch of accelerated flow solutions that bifurcates from the first critical swirl ratio. However, their velocity profile of inviscid accelerated flow (see their figure 2) differs from those in the present work shown in figure 3, indicating significant impact of viscosity on the flow dynamics.
3.1.3. Effect of swirl ratio
The above discussions detail the flow evolution of disturbances to the base Lamb–Oseen vortex for the condition
$\textit{Re}=5000$
and
$\omega =3.85$
. Computations have been performed for
$\omega$
ranging from 3.75 to 3.95 at the same fixed
$\textit{Re}=5000$
. The curve of the growth rate of the perturbation energy in the linear range as a function of
$\omega$
is shown in figure 10. At swirl ratios below a critical number the Lamb–Oseen vortex base flow is asymptotically stable to three-dimensional perturbations. This critical swirl ratio is computed as
$\omega =3.762$
. When
$\omega$
is above this critical swirl ratio, the flow becomes unstable with
$m = 0$
as the unstable perturbation mode. The flow evolves into a final steady non-columnar flow. In cases for
$3.762 \lt \omega \lt 3.91$
all flow features are qualitatively the same as the representative case of
$\textit{Re}=5000$
and
$\omega =3.85$
where the
$m=0$
axisymmetric mode is the dominant linearly unstable mode.
Figure 10 shows that the growth rate of disturbance increases with increasing swirl ratio until it reaches the maximum at
$\omega =3.86$
before decreasing with increasing swirl ratio. The results by Gong (Reference Gong2017) at the same swirl ratio range are also presented in figure 10. There is a surprisingly good agreement between the DNS results and predictions based on the linear stability analysis.
Figure 10. Linear growth rate versus
$\omega$
at fixed
$\textit{Re}= 5000$
: circle, obtained from DNS; black line, predicted by linear theory (Gong Reference Gong2017).
Figure 11. Time history of energy and growth rate of total perturbation velocity and those of mode
$m=0$
and mode
$m=1$
for the case
$\textit{Re}=700$
and
$\omega =4.3$
.
3.2. Case 2: flow dominated by instability of the spiral
$m=1$
mode
The Lamb–Oseen vortex may respond differently to initial random perturbations at different Reynolds numbers and swirl ratios. Results for the case at
$\textit{Re}=700$
and
${\omega }=4.3$
are presented and analysed in this section as another typical case in the
$(\textit{Re}, \omega )$
plane where the axisymmetric
$m=0$
mode is stable while the
$m=1$
spiral mode becomes unstable. Unlike in Case 1, where the final state of the flow is a steady axisymmetric but non-columnar flow, instabilities at this condition lead to a spiral vortex breakdown state. In addition, it is found that the growth of the linearly unstable
$m=1$
mode induces growth in the axisymmetric
$m=0$
component as a forced response due to a weak nonlinear interaction mechanism.
3.2.1. Linear growth of the perturbed flow
Figure 11 plots the instantaneous growth rate and perturbation energy versus time based on DNS simulation data for this case. Like in Case 1, the total energy of the initial perturbations experiences a short decay period during
$0 \lt t \lt 8$
. By the end of the period, the initial energy contained by the stable modes has decayed, while the energy of the dominant unstable mode is the only one that grows at that time. Then the system enters its linear growth stage, during which the growth rate stabilizes at a constant value of
$\sigma _r=0.0059$
from
$t =180$
to
$t = 2000$
. With further increase of time, the linear growth process is saturated with
$\sigma _r$
decreasing to 0 at
$t = 2681$
. From that time on, the flow seems to have reached a ‘steady’ state for approximately 400 non-dimensional time units before the total energy starts to oscillate again and eventually enters into a periodic state beyond
$t=4000$
.
We examine the energy composition in the linear growth period (
$180\lt t\lt 2681$
). Contrary to the findings of Case 1, the
$m=0$
axisymmetric mode is initially stable while the
$m=1$
mode becomes the dominant unstable mode for this case.
The
$m=1$
mode grows with a constant positive growth rate of
$\sigma _{1,r} = 0.0059$
over the entire linear growth period. Figure 12 shows the azimuthal profiles of the extracted axial perturbation velocity
$u_1$
of mode
$m=1$
multiplied by
$e^{-\sigma _{1,r}t}$
at five consecutive time instants in the linear growth period. These profiles are clearly the same profile displaced over an equal phase difference between the evenly spaced time instants, proving that the
$m=1$
mode is a single azimuthal travelling wave moving in the positive
$\theta$
direction. Based on the phase difference (
$\Delta \theta =1.2566$
) and the time interval
$\Delta t=1$
, the phase speed can be determined to be
$\sigma _{1,i}=-1.2566$
. Therefore, we have
$\sigma _{m=1}=0.0059-1.2566i$
.
Figure 12. The phase difference and direction of travelling wave
$u_1$
in the azimuthal direction at
$x=3$
with time interval
$\Delta t=1$
.
Figure 13 shows the axial perturbation velocity profiles of
$g_1$
at
$t=960{-}963$
. For each time instant, the axial profile of
$g_1$
is normalized by its value at the outlet and the radial profile by its value at the central axis. Once again, we see that these normalized profiles collapse into a single line for the four time instants, asserting the validity of the separated time and space modal representation of the perturbations in the form of (2.19). The modal shapes of
$g^{\textit{norm}}_1$
plotted in figure 13 are invariable during the complete linear growth period
$180\lt t\lt 2681$
.
Figure 13. Normalized profiles of
$u_x^\prime$
of mode
$m=1$
at
$t=960{-}963$
: (a) along the axial direction and (b) along the radial direction.
However, figure 13 shows that the radial profiles of
$g^{\textit{norm}}_1$
are very dissimilar for different
$x$
locations and their axial profiles are only slightly similar for different radial locations. Therefore, it is not possible to represent the spatial mode in terms of a separated variable form of
$f(x)g(y)$
in an attempt at analytical study for this case.
The axisymmetric
$m=0$
mode, on the other hand, decays very fast with a negative growth rate of
$\sigma _{0,r} = -0.014$
from
$t=180$
to
$t=520$
. Despite its fast decay, its spatial mode shape can be clearly recovered from the computational data as shown in figure 14. Perfect collapse of the normalized perturbation velocity profiles at the three time instants
$t=240, 245,$
and
$ 250$
and the negative growth rate of its energy affirms it being a stable
$m=0$
eigenmode of the system.
At
$t=520$
, however, the energy of the
$m=0$
component makes a sudden change from decay to growth at a high positive rate of 0.0118 for 1300 non-dimensional time units. The growth rate is exactly two times that of the unstable
$m=1$
mode. Figure 15 shows the normalized axial perturbation velocity profiles of the
$m=0$
component at
$t=1200$
,
$1210$
and
$1220$
in this second linear growth stage. This seemingly unstable
$m=0$
‘mode’, however, is not a true eigenmode of the system. Its spatial profiles shown in figure 15 are distinct from those of the true
$m=0$
mode shown in figure 14, which is stable.
Figure 14. Normalized profiles of
$u^{\prime }_x$
of mode
$m=0$
in the first linear growth period
$180\lt t\lt 520$
: (a) along the axial direction; (b) along the radial direction.
Figure 15. Normalized profiles of
$u^{\prime }_x$
of
$m=0$
component in the second linear growth period
$520\lt t\lt 2000$
: (a) along the axial direction; (b) along the radial direction.
The appearance of this growing
$m=0$
component in the perturbation velocity cannot be explained in a strict linear analysis. It is a response of the system to a forcing by the unstable
$m=1$
mode through a weakly nonlinear effect after the
$m=1$
mode has grown sufficiently large in amplitude.
Following the decomposition in the form of (2.16) for the perturbation modes in the linear range, the
$m=1$
component of the axial perturbation velocity is represented as
\begin{align} &\boldsymbol{u}_{m=1}^{\prime }=\epsilon [ \mathscr{U}_1(x,r)e^{i(\theta +\sigma _{1,i}) t} + \bar {\mathscr{U}}_1(x,r)e^{-i(\theta +\sigma _{1,i}) t} ]e^{\sigma _{1,r} t}, \end{align}
where the overbar stands for complex conjugate;
$\epsilon$
is a small parameter and
$\mathscr{U}(x,r)$
is the function of
$\mathscr{O}(1)$
within the linear small perturbation assumption. Substituting (3.1) into the Navier–Stokes equations and collecting terms of
$\mathscr{O}(\epsilon )$
will yield the linear equations for each component. Within a strict linear analysis, each component is independent of others. However, if we keep terms of
$\mathscr{O}(\epsilon ^2)$
, we will obtain non-homogeneous equations with forcing terms collected from neighbouring modes due to the nonlinear convection terms in the Navier–Stokes equations. In this case, the unstable
$m=1$
mode gives rise to the following forcing term in the equation for the
$m=0$
perturbation velocity component
${\boldsymbol{u}}_{m=0}$
,
\begin{equation} (\partial _t-\mathscr{L}_0){\boldsymbol{u}}_{m=0}^\prime = -\epsilon (\mathscr{U}_1\boldsymbol{\nabla} _{-1}\mathscr{\bar U}+\mathscr{\bar U}_1\boldsymbol{\nabla} _1 \mathscr{U}_1)e^{2\sigma _{1,r}t}, \end{equation}
where
$\mathscr{L}_0$
represents the linear operator of the Navier–Stokes equations,
$\boldsymbol{\nabla} _{\pm 1}=(\partial _r,r^{-1}(\pm i),\partial _x)$
. The forcing term grows exponentially at exactly twice the growth rate of the
$m=1$
mode. The solution will then consist of a homogeneous part, which is the intrinsic stable or unstable
$m=0$
mode depending on the sign of the eigenvalues, and a non-homogeneous part, which is a response due to the forcing term on the right-hand side of (3.2). If we were able to invert the linear operator on the left-hand side of (3.2), we can then write the forced response solution as
\begin{equation} {\boldsymbol{u}}_{m=0}^\prime = -\epsilon (2\sigma _{1,r} I - \mathscr{L}_0)^{-1} (\mathscr{U}_1\boldsymbol{\nabla} _{-1}\mathscr{\bar U}+\mathscr{\bar U}_1\boldsymbol{\nabla} _1 \mathscr{U}_1)e^{2\sigma _{1,r}t} . \end{equation}
The spatial solution in front of the exponential term yields the ‘mode’ shape while the exponential term determines the growth rate of the forced response in the
$m=0$
component by the
$m=1$
mode.
In the present case, the intrinsic
$m=0$
mode is stable and decays to machine zero at
$t=520$
. At the same time, the forced response by the unstable
$m=1$
mode has gained notable amplitude to overpower the small
$\epsilon$
and show the predicted
$2\sigma _{1,r}$
growth rate in figure 11 and the fixed ‘mode’ shape shown in figure 15. To further prove that this forced response in the
$m=0$
component is not an eigenmode of the linearized Navier–Stokes operator, we performed additional computations. We continue the solution from an arbitrary time instant between
$t=1000$
and 2000 when the growing
$m=0$
component has reached significant amplitude (see figure 11) and artificially remove the
$m=1$
mode in the solution at each time step. Invariably we find the total energy of the perturbations decays quickly to zero, proving that this growing
$m=0$
component is indeed dependent on the unstable
$m=1$
mode as a forced response. A similar forced response phenomenon was also found by Sipp & Lebedev (Reference Sipp and Lebedev2007).
Figure 16. Planar projected perturbation velocity vectors and contours of perturbation pressure for the case
$\textit{Re}= 700$
and
$\omega =4.3$
at
$t = 2681$
: (a)
$x{-}z$
plane; (b)
$\theta -r$
plane at
$x=3$
.
Despite its fast growth, the energy of the
$m=0$
component remains almost one order of magnitude lower than that of the primary unstable
$m=1$
mode even at its peak at
$t=2681$
. Therefore, its impact on the growth rate of the
$m=1$
is non-observable. It should be noted that higher azimuthal modes such as
$m=2$
and
$m=3$
are found to decay faster than mode
$m=0$
and their energy decreases to the level of
$10^{-16}$
at
$t=23$
. The growth of the
$m=1$
mode perturbations also forces growth in the higher harmonic
$m=2,3$
components due to the above weakly nonlinear effect. Their evolution process is similar to that of mode
$m=0$
. However, their energy levels stay insignificant compared with even that of the
$m=0$
component. Therefore, we notice that the total perturbation energy is almost all contained in the
$m=1$
mode as shown in figure 11 throughout
$180\lt t\lt 2681$
.
The axial profiles at different
$r$
and the radial profiles at different
$x$
of
$g^{\textit{norm}}_0$
shown in figures 14 and 15 are very dissimilar for both the intrinsic
$m=0$
mode and the forced-response mode, much more so compared with those for Case 1 (see figure 5), which makes the assumption of a separated variable form solution in
$x$
and
$r$
even less plausible. Notice also that the axial profiles of
$g_{0}^{\textit{norm}}$
for the stable
$m=0$
mode in this case consist of three wavelengths while those for the unstable
$m=0$
mode in Case 1 appear to be
$3/2$
waves. This may be due to the increase of the swirl ratio (
$\omega =4.3$
versus
$\omega =3.85$
in Case 1). Because of the above reasons, the long-wave asymptotic analysis by Rusak et al. (Reference Rusak, Granata and Wang2015) assuming separated variable solutions becomes questionable for the complex flow structures appearing at high swirl ratios.
As time approaches
$t = 2681$
, figure 11 shows the growth rates for
$m=0$
and
$m=1$
reduce to 0 at the same time. The flow appears to have reached a ‘steady’ state for approximately 400 non-dimensional time units. It must be noted that the perturbations are really in an unsteady rotational mode because it is dominated by the spiral
$m=1$
mode with an imaginary growth rate although the energy of both the total and the
$m=1$
mode show a constant level.
Figure 16 presents the planar projected velocity vectors together with contours of pressure on the
$x{-}z$
plane and the
$\theta -r$
plane at
$x=3$
for the ‘steady’ state at
$t = 2681$
. The low-pressure regions correspond to accelerated axial and azimuthal velocities. This flow field has directly evolved from the unstable
$m=1$
mode and has a similar spatial structure in its mode shape, demonstrating that the spiral mode of perturbation dominates the process. This state is unstable to three-dimensional disturbances. It undergoes a second nonlinear oscillation process and finally reaches a stable unsteady state after time
$t \gt 4700$
: a temporally periodic spiral vortex breakdown.
3.2.2. The global dynamics and the final state of the flow
As in the study of Case 1, we present in figure 17 the phase portraits for the perturbation velocity at the point
$(x=3$
,
$r=0.5$
,
$\theta =0)$
in different time periods to examine the long-time behaviour and global nonlinear dynamics. Figure 17(a) presents the velocity trajectory in the initial stage when the movement is complicated and irregular in a short transient period
$(0\lt t\lt 62)$
during which various stable perturbation modes decay. The trajectory is eventually drawn towards the velocity point
$(0,0,5\times 10^{-6})$
from
$t = 62$
to
$t = 180$
(black line). Starting from
$t=180$
, however, the unstable
$m=1$
mode takes over and grows exponentially. Consequently, the trajectory of the perturbation velocity during
$180 \lt t \lt 2000$
(blue line) traverses a much larger range so that the trajectory over the time period
$0 \lt t \lt 180$
(black line) is almost invisible in figure 17(b). During that linear growth stage, the spiral mode
$m = 1$
grows with a constant rate
$\sigma _r= 0.0059$
and dominates the flow dynamics. The perturbation velocity is seen to move along a spiral orbit expanding with time. The radial and azimuthal velocities grow at greater rates compared with the axial velocity, which represents a tendency towards asymmetry.
Figure 17(c) zooms further into the time period
$2000 \lt t \lt 2681$
, during which the trajectory continues to move along an outward spiral orbit but with the growth rate decreasing from 0.0059 to 0, reaching a limit cycle at
$t=2681$
. Compared with the time period
$180 \lt t \lt 2000$
(dark blue line), the orbit to the final limit cycle (light blue line) is of much larger size, which corresponds to the nonlinear evolution process. The radial and azimuthal perturbations remain large. This limit cycle, however, is unstable to three-dimensional disturbances. A growing nonlinear oscillation around the final limit cycle of figure 17(c) begins at time
$t=2681$
and continues to
$t=4700$
as shown in figure 17(d). After
$t =4700$
, the flow reaches a stable final state with a temporally periodic behaviour. The velocity trajectory takes the unstable limit cycle at
$t=2681$
as a centre ring and forms a complex but regular orbit shown in figure 17(e), which is the attractor of the flow dynamics.
Figure 18(a) shows the time-periodic axial perturbation velocity for mode
$m=0$
and mode
$m=1$
at
$x=2, r=0.2$
, and
$\theta =0$
versus time. The
$m=1$
velocity experiences a fast oscillation at a high frequency and the
$m=0$
velocity experiences a slow oscillation at a low frequency. Figure 18(b) presents the corresponding fast Fourier transform (FFT) of
$u_x^\prime$
for the two modes. The distinct fundamental frequencies are clearly identified. The frequency of the slow variation of the
$m=0$
signal is
$f_E=0.0073$
, which is also the frequency of the final total energy oscillation. Since the
$m=0$
velocity represents the azimuthal average, the low frequency reflects periodic behaviour in the axial direction, see Supplementary movie 2, available at https://doi.org/10.1017/jfm.2025.10595.
The frequency of the rapid variation of the
$m=1$
signal is
$f_u=0.2089$
. This high frequency reflects the fast spiral motion in the azimuthal direction. This spiral frequency is clearly related to the angular frequency
$|\sigma _{1,i}| =1.2566$
of the
$m=1$
mode identified in § 3.2.1 during its linear growing phase between
$t=180$
and 2681. It means that the fast spiral motion of the final dynamic equilibrium state is inherited from the linearly unstable
$m=1$
spiral mode. Moreover, the low frequency originates from the
$m=0$
mode, which also modulates and influences the dynamics of the
$m=1$
field.
Simultaneously, the spiral mode is subject to spatial modulation by the slow oscillation of the
$m=0$
mode in the axial direction, resulting in a visible envelope in the signal in figure 18(a). The red line in figure 18(b) also shows the modulated frequency centred at
$f_u$
with sidebands around
$f=f_u \pm nf_E,\ n=1,2,\ldots$
, clear evidence that the envelope frequency originates from the
$m=0$
field. The presence of multiple harmonically spaced sidebands in the FFT analysis confirms the frequency-modulated behaviour of the signal.
Figure 17. The velocity trajectory for
$\textit{Re}= 700$
and
$\omega = 4.3$
: (a)
$t$
= 0–180, the initial transient stage; (b)
$t =$
0–2000, the linear growth stage; (c)
$t =$
180–2681, the complete linear and nonlinear stage; (d)
$t =$
2681–4700, second instability growth stage; (e) unstable limit cycle at
$t=2681$
and stable final periodic behaviour.
Figure 18. (a) The time-periodic perturbation axial velocity at
$x=2$
,
$r=0.2$
and
$\theta =0$
versus time. (b) The FFT of
$u_x^\prime$
. The black line represents the extracted mode
$m=0$
and the red line represents the extracted mode
$m=1$
.
The energy of the perturbation velocity experiences periodic growth and decay in a period
$T_E=1/f_E=136.4$
due to the final oscillatory flow, see figure 11. The total energy and those of the four extracted azimuthal wavenumbers
$m=0- 3$
of the final flow state are presented in the form of polar plots in figure 19. The polar angle represents time in a period and the radius represents the perturbation energy
$E$
. Time
$t_0=5085.72$
is chosen as
$0^\circ$
so that eccentric circles in figure 19 show the fluctuating energy of various modes in a period
$t_0+T_E$
. The
$m=1$
spiral mode is the dominant mode. The energy amplitudes of
$m=2$
and
$m=3$
are almost invisible compared with those of components
$m=0$
and
$m=1$
. All azimuthal modes appear as a frequency lock-in phenomenon with
$f_E=0.0073$
, showing a good periodic stability of the dynamical attractor.
Figure 19. Perturbation energy of extracted azimuthal wavenumbers
$m=0\sim 3$
together with total energy in the period
$T_E=136.4$
.
To visualize the final state of the vortex flow, streaklines are obtained by releasing particles at equidistant azimuthal positions on a circle of radius
$r=0.1$
at the inlet. Based on the period
$T_u=4.7870$
for the spiral motion of the vortex flow, two sets of streaklines released consecutively with a separation time of
$\Delta t=T_u/2= 2.3935$
are shown in figure 20, which clearly illustrate two whirling vortex structures in opposing positions (
$180^\circ$
phase difference) for the two time instants. Within each period
$T_u$
of the fast spiral motion of the vortex, the tube formed by the streaklines maintains almost the same size and shape.
Streaklines obtained at two instants
$t=5030.19$
and
$t=5096.51$
are plotted in figure 21, corresponding to the maximum and minimum perturbation energy levels, respectively, within the low-frequency energy oscillation cycle
$T_E$
. It is found that the cross-section of the streaktube with the maximum perturbation energy is bigger than that of the streaktube with the minimum perturbation energy at the outlet and also, has a larger wave amplitude. The spiral vortex bulges as the total perturbation energy peaks. The spatial shape of streaktubes changes at the low frequency
$f_E$
during the long period
$T_E$
. It reflects the slow periodic behaviour of the extracted
$m=1$
spiral mode modulated by the oscillation of the
$m=0$
mode.
Continuous streaklines based on the total velocity field (Supplementary movie 1) and those consisting of only the
$m=0$
(Supplementary movie 2) and
$m=1$
(Supplementary movie 3) modes of perturbations are also provided. The
$m=1$
mode is the dominant mode and its streaktube shows a fast spiral motion around the pipe at the frequency of
$f_u$
. The streaktube of the
$m=0$
mode is axisymmetric and exhibits small and slow bulging motions along the axial direction corresponding to the oscillation of the total perturbation energy at frequency
$f_E$
. Together they form the total vortex flow field since the
$m=0$
and
$m=1$
modes constitute most of the total energy of the perturbation velocity field.
Figure 20. Streaklines with
$180^{\circ}$
phase shift in a small period
$T_u=4.7870$
at time
$t=5094.12$
and
$t=5096.51$
in three-dimensional space.
Figure 21. Streaklines at
$t=5030.19$
and
$t=5096.51$
in the global period
$T_E=136.4$
at the
$x{-}z$
plane.
Figure 22. (a) Isosurface of vorticity
$|\varOmega | =1.61$
. (b) Contours of axial vorticity at
$x{-}z$
plane. (c) Contours of azimuthal vorticity at
$x{-}z$
plane.
Figure 22 presents a different set of visualization of the final spiral breakdown state at
$t = 5200$
. The yellow surfaces in figure 22(a) are isosurfaces of the vorticity magnitude
$|\varOmega | =1.61$
, which are overlapped with the contours of the axial vorticity component
$\varOmega _x$
in the azimuthal plane. Figures 22(b) and 22(c) plot in the
$x{-}z$
plane the contours of the axial vorticity component
$\varOmega _x$
and those of the azimuthal vorticity component
$\varOmega _{\theta }$
, respectively, both of which reveal asymmetric alternating signs of the vorticity components relative to the axis. These figures help illustrate the flow structure of the ultimate spiral vortex breakdown state.
3.3. Stability boundary in the
$(\textit{Re}, \omega )$
plane
In order to determine the stable and unstable regions in the parameter
$(\textit{Re},\omega )$
plane, one may search for ranges of unstable critical swirl ratios for a given Reynolds number
$\textit{Re}$
. Alternatively, we may start with a given swirl ratio
$\omega$
and search for the critical Reynolds number and also determine the mode of instability. Computations have been done in such a way to map out the stability boundaries in the
$(\textit{Re},\omega )$
plane for the Lamb–Oseen vortex. The results are shown in figure 23.
The curves in the figure delineate the neutral stability boundaries for the marked
$m=0$
or
$m=1$
modes. The solid
$m=0$
boundary and the dash–dotted
$m=1$
boundary mark the boundary of critical Reynolds number, below which the base flow is found to be stable to all small disturbances. In other words, they mark the most unstable mode for each
$\omega$
.
Figure 23. Stability boundaries in the (
$\textit{Re}$
,
$\omega$
) plane for the Lamb–Oseen vortex in a pipe of length
$L = 6$
: the solid and dashed lines are neutral boundaries for the
$m=0$
axisymmetric mode with
$\sigma =0$
and
$\sigma _r=0,\sigma _i \ne 0$
, respectively. The dash–dotted line is neutral boundary for the
$m=1$
spiral mode.
For example, at
$\omega =3.85$
, the flow is found to be stable to any small initial random perturbations for Reynolds number
$\textit{Re}\lt 1484$
. At
$\textit{Re}=1484$
, the growth rate of mode
$m=0$
becomes zero, while all
$m \geqslant 1$
modes stay negative. Once the Reynolds number crosses the minimum threshold, the
$m=0$
mode becomes unstable, while all higher-mode disturbances are found to be stable. Case 1 (
$\omega =3.85$
,
$\textit{Re}=5000$
) examined in § 3.1 is one such example in that region. This is found to be the general behaviour of the flow for small swirl ratios.
The critical Reynolds number decreases with increasing swirl ratio initially for small swirl ratios, but then the trend reverses starting at approximately
$\omega =3.875$
. For conditions with swirl ratios higher than
$\omega =3.91$
, the most unstable mode switches to
$m=1$
spiral mode shown by the dash–dotted line. The critical Reynolds number starts to decrease again and monotonically so with increasing swirl ratio. More importantly, in the vicinity of this intersection point, there exists a region at higher Reynolds numbers where both the axisymmetric mode with real growth rate and spiral modes are unstable. Above the dash–dotted line, two regions are identified. In the region between the dash–dotted and dashed lines, the
$m=0$
mode is stable while the
$m=1$
mode is unstable, although the unstable
$m=1$
may induce growth in the
$m=0$
component due to nonlinear forcing. Case 2 (
$\omega =4.3$
,
$\textit{Re}=700$
) examined in § 3.2 is such a typical case. In the high Reynolds number region above the dashed line, both the
$m=0$
and
$m=1$
modes are found to be unstable. In addition, the unstable
$m=0$
mode there contains an imaginary part.
One may conclude from this study that higher swirl ratios favour spiral types of instability and eventual spiral vortex breakdown.
4. Mechanism for instability
The appearance of an instability mode is tightly related to the energy transfer between the perturbations and the base flow inside the flow domain and on its boundaries. The Reynolds–Orr equation provides a useful tool for such analyses (Wang Reference Wang2008; Gong Reference Gong2017).
4.1. The Reynolds–Orr equation
The Reynolds–Orr equation is obtained by multiplying the transport equations for the perturbation velocity with the perturbation velocity and integrating them over the flow domain. Note that although a constant source term is added in (2.4), it does not affect the production of perturbation energy, as the source term is balanced in the equations for the base flow. In the cylindrical coordinate system, it is written as
\begin{align} \frac {\textrm{d}E(t)}{\textrm{d}t}=&-\int _{V_C}{ \big( u_r^\prime ,u_{\theta }^{\prime },u_{x}^{\prime } \big) B \big( u_{r}^{\prime },u_{\theta }^{\prime },u_{x}^{\prime } \big) ^T\textrm{d}V}+\int _{S_1}{u_x^\prime p^\prime \textrm{d}S_1}-\int _{S_2}{u_x^\prime p^\prime \textrm{d}S_2}\nonumber\\&+\frac {1}{2}\int _{S_1}{u_x \big( {u_r^\prime }^2+{u_\theta ^\prime }^2+{u_x^\prime }^2 \big) \textrm{d}S_1}-\frac {1}{2}\int _{S_2}{u_x \big( {u_r^\prime }^2+{u_\theta ^\prime }^2+{u_x^\prime }^2 \big) \textrm{d}S_2}\nonumber\\&-\frac {1}{\textit{Re}} \!\!\!\int _{V_C}{ \left( \! | \boldsymbol{\nabla }u_{r}^\prime |^2+| \boldsymbol{\nabla }u_\theta ^\prime |^2+| \boldsymbol{\nabla }u_{x}^\prime |^2 +\frac {2u_r^\prime }{r^2}\frac {\partial u_\theta ^\prime }{\partial \theta }+\frac {{u_{r}^\prime }^2}{r^2}-\frac {2{u_{\theta }^\prime }^2}{r^2}\frac {\partial u_{r}^\prime }{\partial \theta }+\frac {{u_{\theta }^\prime }^2}{r^2} \! \right) \textrm{d}V}\nonumber\\&+\frac {1}{\textit{Re}}\int _{\partial {V_C}}{\Bigg( u_r^\prime \frac {\partial u_r^\prime }{\partial n}+u_{\theta }^\prime \frac {\partial u_\theta ^\prime }{\partial n}+u_x^\prime \frac {\partial u_x^\prime }{\partial n}\Bigg) \textrm{d}S}. \end{align}
Here,
$V_C$
is the computational domain and
$\partial {V_C}$
is the pipe wall together with inlet and outlet surfaces. Here
$S_1$
and
$S_2$
represent the pipe inlet and outlet sections, respectively. Here
$\partial /\partial {n}$
is the directional derivative in the outer normal direction of
$\partial {V_C}$
. The symmetric strain-rate tensor
$B$
of the base Lamb–Oseen vortex is
\begin{equation} B=\frac {\omega \big ( e^{-br^2}+br^2e^{-br^2}-1 \big )}{r^2}\left [ \begin{matrix} 0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\\ \end{matrix} \right ]\!. \end{equation}
The two integrals on
$S_1$
vanish because there is no velocity perturbation at the inlet. Substituting
$B$
into (4.1), we can express the growth rate
$\sigma _r$
of the perturbation velocity defined in (2.24) by the following equation:
\begin{equation} \sigma _r=\frac {1}{2E\left ( t \right )}\frac {\textrm{d}E(t)}{\textrm{d}t} =\sigma _1 + \sigma _2+ \sigma _3+ \sigma _4, \end{equation}
where
\begin{align} \sigma _1=&\frac {1}{E\left ( t \right )}\int _{Vc}{\left [ \frac {\omega }{r^2}\big ( 1-e^{-br^2} \big ) -\omega be^{-br^2} \right ] u_{r}^\prime u_{\theta }^\prime \textrm{d}V}, \end{align}
\begin{align} \sigma _2=&-\frac {1}{2E\left ( t \right )}\int _{S_2}{u_{x}^\prime p^\prime \textrm{d}S_2}, \end{align}
\begin{align} \sigma _3=&-\frac {1}{4E(t)}\int _{S_2}{u_x\left ( {u_r^\prime }^2+{u_{\theta }^\prime }^2+{u_{x}^\prime }^2 \right ) \textrm{d}S_2}, \end{align}
\begin{align} \sigma _4=&-\frac {1}{2E\left ( t \right ) \textit{Re}}\int _{V_C}{\bigg ( | \boldsymbol{\nabla }u_{r}^\prime |^2+| \boldsymbol{\nabla }u_{\theta }^\prime |^2+| \boldsymbol{\nabla }u_{x}^\prime |^2 } \nonumber \\ &+\frac {2u_{r}^\prime }{r^2}\frac {\partial u_{\theta }^\prime }{\partial \theta }+\frac {{u_{r}^\prime }^2}{r^2}-\frac {2{u_{\theta }^\prime }^2}{r^2}\frac {\partial u_{r}^\prime }{\partial \theta }+\frac {{u_{\theta }^\prime }^2}{r^2} \bigg ) \textrm{d}V \\ &+\frac {1}{2E\left ( t \right ) \textit{Re}}\int _{\partial Vc}{\left ( u_{r}^\prime \frac {\partial u_{r}^\prime }{\partial n}+u_{\theta }^\prime \frac {\partial u^\prime _{\theta }}{\partial n}+u_{x}^\prime \frac {\partial u_{x}^\prime }{\partial n} \right ) \textrm{d}S \nonumber }. \end{align}
The sum of the above four terms gives the total growth rate of the perturbation energy:
$\sigma _1$
represents the growth rate of the integrated perturbation energy produced by the vortex core;
$\sigma _2$
is the growth rate due to pressure work at the pipe outlet;
$\sigma _3$
is the growth rate due to the convection of perturbation energy out of the pipe, which is always negative unless there is back flow at the exit;
$\sigma _4$
stands for the growth rate due to viscous effect in the pipe domain and on the pipe surfaces.
Note that in the classical stability theory, terms related to the outlet will vanish because of the periodic inlet–outlet boundary condition. We will see that the energy transfers at the outlet boundary due to pressure and convection play a critical role in determining the stability of a vortex flow in a pipe of finite length because of the non-periodic flow condition at the inlet and outlet boundaries. In the following subsections we apply the above analysis to the two typical instability cases studied in the previous section.
4.2. Energy balance and sources of instability
We make use of the Reynolds–Orr equation to categorize the production of the perturbation energy and thus identify the sources of energy production responsible for the observed instabilities of both Cases 1 and 2.
Figure 24. Growth rates
$\sigma _1$
,
$\sigma _2$
,
$\sigma _3$
,
$\sigma _4$
,
$\sigma _{sum}$
and
$\sigma _E$
versus
$t$
for case
$\textit{Re}= 5000$
and
$\omega = 3.85$
.
Figure 24 plots the time history of
$\sigma _1$
,
$\sigma _2$
,
$\sigma _3$
,
$\sigma _4$
and their sum
$\sigma _{sum}$
for Case 1 (
$\textit{Re}= 5000$
and
$\omega = 3.85$
). Also shown for comparison is the growth rate
$\sigma _E$
calculated directly from the DNS perturbation velocity (black circles). The sum of the four contributions in the Reynolds–Orr equation fully matches that of
$\sigma _E$
, proving perfect energy conservation in accounting for the sources of energy production in (4.1). The line of
$\sigma _1$
represents kinetic energy exchange between the base flow and the perturbation. During the linear growth stage,
$\sigma _1$
is negative, indicating that the vortex core is absorbing energy and therefore has the effect of stabilizing the flow, which is consistent with Rayleigh’s criterion. However, the energy production from pressure work
$\sigma _2$
and energy loss from flux
$\sigma _3$
are both significant and must be accounted for because of the non-periodic inlet and exit boundary conditions in the case of a finite-length pipe. The computational results show that the pressure work in the linear range is always positive and significantly higher than the energy loss due to convection at the exit of the pipe (
$\sigma _2\gt -\sigma _3$
), resulting in a net positive gain
$\sigma _2 +\sigma _3$
that overpowers the energy drain
$\sigma _1$
by the vortex core and the viscous dissipation
$\sigma _4$
. The net energy balance
$\sigma _{sum}$
is thus positive, causing the observed instability. The viscous effect tends to stabilize the flow because
$\sigma _4$
is always negative, which also explains the fact of an increased growth rate with increasing Reynolds number (lower viscous dissipation).
Figure 25. Linear growth rates
$\sigma _1$
,
$\sigma _2$
,
$\sigma _3$
,
$\sigma _4$
,
$\sigma _{sum}$
and
$\sigma _E$
versus
$\omega$
at
$\textit{Re}= 5000$
(the unstable mode
$m=0$
dominates the flow dynamics).
Figure 25 presents the same terms of energy balance but their values are taken from the linear growth stages of flow conditions with swirl ratios in the range
$3.762\lt \omega \lt 3.91$
and the fixed Reynolds number
$\textit{Re}= 5000$
, where the flow is determined to be unstable to the axisymmetric
$m=0$
mode of perturbations in § 3.3. The energy production due to pressure
$\sigma _2$
and the loss due to convection
$\sigma _3$
at the exit of the pipe remain as the dominant terms. The energy production rate due to viscous effect
$\sigma _4$
remains negative and almost constant for all swirl ratios because the viscous dissipative effect depends primarily on the Reynolds number. Although the pressure work decreases with increasing swirl ratio, it still overpowers all other negative energy sources, resulting in a net positive growth rate
$\sigma _{sum}$
for the total perturbation energy for all cases in figure 25. It is interesting to notice in figure 25 that the energy production rate
$\sigma _1$
due to the vortex core becomes less negative and reaches zero as the swirl ratio
$\omega$
increases towards
$3.91$
which is the critical swirl number.
Figure 26 plots the time history of
$\sigma _1$
,
$\sigma _2$
,
$\sigma _3$
,
$\sigma _4$
and their sum
$\sigma _{sum}$
for Case 2 (
$\textit{Re}=700$
and
$\omega =4.3$
), which belongs to the region of spiral mode instability with
$\omega \gt 3.91$
in figure 23. Similarly to the cases where the
$m=0$
mode of instability dominates, the energy production rates
$\sigma _2$
and
$\sigma _3$
are still of the highest magnitudes. The energy production rate by pressure work
$\sigma _2$
is now slightly less than the energy loss due to convection at the pipe exit (
$\sigma _2\lt -\sigma _3$
). However, the energy production
$\sigma _1$
driven by the Lamb–Oseen vortex core is now positive and of the same order of magnitude as
$\sigma _2$
and
$\sigma _3$
. Together with the positive work by the vortex core, the pressure work at the outlet overcomes the energy loss due to convection at the outlet and viscous dissipation, leading to the spiral mode instability. The vortex core appears to promote growth of non-axisymmetric spiral perturbations at high swirl ratios.
Wang et al. (Reference Wang, Rusak, Gong and Liu2016) revealed similar energy transfer mechanisms for the vortex flow of the solid-body rotation type. We note the following general observations regarding the stability mechanisms that apply to a swirling flow with non-zero axial flow velocity within a tube of finite-length.
The vortex core in a swirling flow serves as a waveguide to sustain perturbation propagation. In an infinitely long (or sufficiently long) waveguide, these waves appear as helical neutral travelling waves as represented by the spatial Fourier modes in the axial direction in a normal mode analysis for which the Rayleigh criterion applies. However, when the waveguide length is comparable to the perturbation’s wavelength, the imposed asymmetric boundary conditions at inlet and outlet will introduce additional energy transfer at the inlet and outlet due to pressure work and flow convection, which can lead to instability. In addition, swirling flows are capable of sustaining upstream propagating helical waves. Imposition of non-periodic boundary conditions strongly affects the wave shape of the perturbation velocity field. Typical axisymmetric and spiral mode shapes as shown in figures 3 and 16 can be considered as the deformed helicity waves (see also figures 5, 7, 13, 14 and 15). They effectively affect the generation of the perturbation energy by the base vortex core.
Axial inhomogeneity can exist in practical flow situations such as the flow in a combustion chamber or a swirling jet where the flow rotation is diminished in the axial direction by viscous forces.
Figure 26. Growth rates
$\sigma _1$
,
$\sigma _2$
,
$\sigma _3$
,
$\sigma _4$
,
$\sigma _{sum}$
and
$\sigma _E$
as a function of
$t$
for case
$\textit{Re}= 700$
and
$\omega = 4.3$
.
5. Effect of pipe length
Wang & Rusak (Reference Wang and Rusak1996a
), Wang et al. (Reference Wang, Rusak, Gong and Liu2016) and Gong (Reference Gong2017) demonstrated that the key to the existence of instability of a vortex in a pipe is the finite pipe length. Their analyses focused on examples with
$L = 6$
(length to radius since
$L$
is non-dimensionalized by the pipe radius) for both the solid-body rotation and the Lamb–Oseen vortex profiles. Feng et al. (Reference Feng, Liu, Rusak and Wang2018) studied the solid-body rotation in an
$L=2$
pipe by DNS. A shorter length was used due to limitations of the computer code and computational resources at the time. We have since improved the efficiency of the computer code, including incorporating parallel computation as described in § 2.5. Together with increased computational power, the new code enabled us to perform DNS simulations of the Lamb–Oseen vortex profile with sufficient grid resolution in the large domain associated with
$L=6$
. The effect of pipe length is investigated in this section by performing additional computations for
$L=5$
and
$L=7$
conditions.
We examine the same two typical instability cases as those discussed in §§ 3.1 and 3.2. For the case where the axisymmetric instability mode dominates, computations are repeated by varying the swirl ratio
$\omega$
at the fixed Reynolds number
$\textit{Re}=5000$
. The linear growth rate of the perturbation energy as a function of
$\omega$
for three different pipe lengths (
$L=5-7$
) is shown in figure 27 (see also figure 10). As the pipe length increases, the range of swirl ratios where flow is unstable becomes narrower and the linear growth rate decreases where instability occurs. This indicates that as the pipe length continues to increase (approaching an infinitely long pipe or a pipe with periodic inlet–outlet boundary conditions), the stability behaviour approaches the result of classic stability analysis as predicted by the Rayleigh criterion: the Lamb–Oseen vortex is centrifugally stable. More computations are needed in the future for full proof of the asymptotic limit. But the trend is clear here.
Figure 28 plots the axial perturbation velocity profiles of
$g_0^{\textit{norm}}$
along the central axis in the linear growth stage for the three different pipe lengths obtained at
$\omega$
values that correspond to the maximum growth rate for each pipe length. The profiles are normalized by their corresponding values at the outlet. The axial distance is normalized by the pipe length
$L$
. These normalized mode shapes of the perturbation velocity appear to be little affected by the pipe length.
Figure 27. Linear growth rate versus
$\omega$
for three different pipe lengths at
$\textit{Re}=5000$
.
We next examine the same case as in § 3.2, where the flow is dominated by an
$m=1$
spiral instability at
$\textit{Re}=700$
and
$\omega =4.3$
. The linear growth rates of the perturbation mode
$m=1$
for different
$L$
are listed in table 2. The trend is the same as that of the axisymmetric mode: as the pipe length increases, the growth rate decreases, indicating a more stable flow. Figure 29 plots the axial perturbation velocity profiles of
$g_1^{\textit{norm}}$
at
$r=0.5$
and
$\theta =0$
in the linear growth stage for the three pipe lengths. These mode shapes of
$m=1$
instability mode show more noticeable variations with pipe length compared with the
$m=0$
instability mode shown in figure 28. The mode shape for the shorter
$L=5$
pipe shows significantly more oscillatory behaviour, influenced by the non-periodic inlet–outlet boundary conditions.
Table 2. Linear growth rate at
$\textit{Re}=700$
and
$\omega =4.3$
for three different pipe lengths.
Figure 28. Normalized profiles of
$u^{\prime }_x$
at the central axis of the linear perturbation mode
$m=0$
for three different pipe lengths.
Figure 29. Normalized profiles of
$u^{\prime }_x$
of the linear perturbation mode
$m=1$
for three different pipe lengths.
6. Conclusions
Direct numerical simulations are performed on the evolution of perturbations to a Lamb–Oseen vortex in a pipe of finite length to study the stability of the base flow and its linear and nonlinear dynamics. Contrary to conclusions based on the Rayleigh criterion (Rayleigh Reference Rayleigh1917), the Lamb–Oseen vortex in a pipe of finite length is found to be unstable for certain combinations of Reynolds number and swirl ratio. A methodology based on azimuthal Fourier decomposition is proposed to analyse the evolution of the perturbation velocity field obtained by the DNS. Axisymmetric (
$m=0$
) and higher azimuthal (
$m=1, 2, \ldots $
) types of instability modes are identified in the evolution of the perturbations and final states of the perturbed Lamb–Oseen vortex flow. Total energy and the energies contained in the different modes of the perturbation velocity field and their instantaneous growth rates are examined during the full linear and nonlinear evolution process and their connections are made to the final state of the flow.
The linear stability boundaries are determined in the
$(\textit{Re},\omega )$
plane. For a given swirl ratio
$\omega$
, the flow is found to become unstable to small initial random perturbations when the Reynolds number
$\textit{Re}$
is above a critical value. However, different modes are found to be responsible for the initial linear instability and to influence its subsequent nonlinear evolution and the final equilibrium state of the perturbed vortex flow. For small swirl ratios, it is the axisymmetric
$m=0$
mode that dominates the initial instability and leads to a final axisymmetric but non-columnar flow state. For large swirl ratios, it is the non-axisymmetric
$m=1$
spiral mode that dominates the initial instability and leads to a final unsteady spiral vortex breakdown state. These two scenarios are exemplified by two cases and are analysed in detail as Case 1 and Case 2, respectively.
For Case 1 (
$\textit{Re}=5000$
and
$\omega =3.85$
), the axisymmetric mode
$m=0$
grows at a constant rate
$\sigma _{m=0}=0.0133$
while all other modes decay in the linear growth stage. The extracted non-dominant mode
$m=1$
has a negative and complex growth rate
$\sigma _{m=1}=-0.0039-1.2665i$
. The existence of the non-zero imaginary part of growth rate indicates it is a spiral travelling wave in the azimuthal direction. This spiral mode and all other higher modes of perturbations decay fast until they remain at the machine zero level. The spatial shapes of
$m=0$
and
$m=1$
modes,
$g_0(x,r,t)$
and
$g_1(x,r,t)$
, are both found to be time-invariant functions (
$g_0^{\textit{norm}}(x,r)$
and
$g_1^{\textit{norm}}(x,r)$
) in the linear range when normalized properly by the assumed linear growth behaviour described by (2.19). For the axisymmetric mode, the mode shape
$g_0^{\textit{norm}}(x,r)$
may also be well represented in separated variables form in
$x$
and
$r$
because the axial profiles of
$g_0^{\textit{norm}}(x,r)$
at different
$r$
locations and the radial profiles of
$g_0^{\textit{norm}}(x,r)$
at different
$x$
locations are, respectively, similar. This is, however, not the case for the
$m=1$
spiral mode. The radial profiles of
$g_1^{\textit{norm}}(x,r)$
are quite dissimilar at different
$x$
locations although the profiles along the
$x$
direction show a more similar wave-like shape. After the saturated nonlinear process, the perturbation evolves to an attractor of the global flow dynamics: a non-columnar accelerated flow, which is stable to any further three-dimensional perturbations. During the whole evolution process including the nonlinear stage, the axisymmetric mode dominates the perturbation velocity field.
For Case 2 (
$\textit{Re}=700$
and
$\omega =4.3$
), a representative case with high swirl ratios in the
$(\textit{Re},\omega )$
plane, the asymmetric
$m=1$
spiral mode is found to be the dominant instability mode. The
$m=0$
and
$m=2,3, \ldots $
modes are all found to be initially stable. The unstable spiral mode has a linear growth rate
$\sigma _{m=1}=0.0059-1.2566i$
and the stable axisymmetric mode has a negative linear growth rate
$\sigma _{m=0}=-0.014$
. With the negative growth rate, the initial symmetric
$m=0$
mode appears to decay fast before its energy suddenly starts to increase at a positive rate
$\sigma _{m=0}=0.0118$
. This sudden change of behaviour is identified as a mean-flow response forced by the unstable
$m=1$
mode through a weakly nonlinear effect after the perturbations of the
$m=1$
mode have grown sufficiently large in amplitude. The normalized profiles
$g_1^{\textit{norm}}$
of the
$m=1$
perturbation mode are found to be time invariant in the complete linear growth stage. The normalized profiles
$g_0^{\textit{norm}}$
of the
$m=0$
perturbation mode are found to be time invariant in the two separate linear growth stages, the first of which is in the decay range of the initial
$m=0$
mode, the second is in the growing stage forced by the unstable
$m=1$
mode. Their spatial structures are very different in the two stages, revealing that the axisymmetric
$m=0$
components in the two stages are really two different axisymmetric mechanisms. The first has a three axial wavelengths shape while the second appears to be a 1/4 axial wave shape.
The growth of the
$m=1$
mode perturbations also forces growth in the higher harmonics
$m =2,3, \ldots $
due to the same weakly nonlinear effect. However, their energy levels are insignificant. Despite its fast growth rate, the energy of the forced
$m=0$
response stays almost one order of magnitude below that of the unstable asymmetric
$m=1$
spiral mode throughout the linear and nonlinear stages and until the flow is drawn into its final attractor of the dynamic system: a spiral vortex breakdown state, in which the vortex is seen to move in a fast spiral motion at frequency
$f_u=0.2089$
while pulsating in size and range of the spiral motion giving rise to a low-frequency oscillation of the total perturbation energy at
$f_E=0.0073$
. The fast spiral frequency seems to have been inherited from that of the
$m=1$
linear perturbation mode and the slow energy frequency is from the secondary instability mode of the mean flow,
$m=0$
.
The energy method based on the Reynolds–Orr equation is used to identify the sources of energy production responsible for the observed instabilities based on the DNS simulation data. Four terms of perturbation energy production are identified: the perturbation energy produced by the vortex core; the pressure work at the pipe outlet; the energy that is carried out of the pipe by the flow; the energy dissipation due to viscosity. The viscous effect always dissipates the perturbation energy and depends primarily on the Reynolds number. In the classical stability theory for the vortex flow in an infinitely long pipe, the pressure work and energy convection terms cancel at the inlet and outlet boundaries due to the assumed periodic or homogeneous boundary conditions, thus the energy production by the vortex core is the only potential source of instability to which the Rayleigh criterion can be applied for an inviscid flow. In the case of a finite pipe, however, the pressure work and energy convection at the pipe exit are found to be of the highest magnitudes among the four energy production sources and thus play a critical role in determining the stability of the flow.
For the conditions represented by Case 1, where the swirl ratio is small and the axisymmetric
$m=0$
instability dominates, the energy production by the vortex core is negative and small, consistent with the Rayleigh criterion. However, the energy production from the pressure work at the outlet overpowers all other negative energy sources resulting in a net positive growth rate for the perturbations.
For the conditions represented by Case 2, where the swirl ratio is high and the
$m=1$
instability is found to be dominant, the energy production by the vortex core becomes positive and comparable in magnitude to those of the pressure work and energy loss due to flow convection at the pipe exit.
Imposition of the non-periodic boundary conditions not only introduces the pressure work and energy convection terms but also strongly affects the waveshape of the perturbation velocity field, which in turn can change the sign of the energy production by the vortex core. In the high swirl ratio case, the energy produced by the vortex core together with the positive work by the pressure at the outlet overcomes the energy loss due to convection and viscous dissipation, leading to a spiral mode instability.
The effect of pipe length is investigated by performing computation of the Lamb–Oseen vortex in different pipe lengths
$L=5-7$
. For cases where the axisymmetric instability mode dominates, the linear growth rate of the perturbation energy as a function of
$\omega$
decreases as the pipe length increases, and the unstable range of
$\omega$
becomes narrower. For cases where the spiral instability mode dominates, the linear growth rate of the perturbation energy also decreases as the pipe length increases. For a shorter pipe
$L=5$
, the mode shapes show more oscillatory variation due to the impact of the non-periodic inlet–outlet boundary conditions. The results suggest that as the pipe length increases and approaches an infinitely long pipe (or a pipe with periodic inlet–outlet boundary conditions), the Lamb–Oseen vortex becomes centrifugally stable, consistent with the classical Rayleigh criterion.
The Lamb–Oseen vortex is a model for practical vortical flows. Practical situations introduce boundary conditions over finite domains. The present work demonstrates the use of DNS for the study of stability and global dynamics of vortex flows. The study reveals important findings on the nature of the Lamb–Oseen vortex that provide critical insights into the physical mechanisms of unstable vortex flows in a finite domain. The stability analysis is general and can be applied to other vortices of interest. Moreover, the results help connect unstable base flow to experimentally observable transitions such as breakdown phenomena and offer guidance for controlling flow regimes in practical systems such as combustors where flow instability is utilized to enhance mixing and flame stabilization.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10595.
Acknowledgements
The authors would like to thank the referees, whose critique and insightful questions have significantly improved the quality of this paper.
Funding
The first and second authors received support by the National Natural Science Foundation of China (NSFC) (Y.-P.S., grant number 11988102).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Effect of the source term
In § 2.2, we discussed the difficulty of using the Lamb–Oseen vortex with uniform axial velocity and zero radial velocity as a base flow for stability analysis because of the slow expansion of its vortex core and dissipation of vorticity magnitude due to viscous diffusion.
There are three options to remedy this problem.
-
(i) Directly use the unsteady Lamb–Oseen vortex as the base flow. Starting from a given vortex size and swirl ratio and Reynolds number, add disturbances to the base flow and perform DNS simulations to observe how the flow and disturbances evolve.
-
(ii) Fix the Lamb–Oseen vortex profile as boundary conditions at the pipe inlet, impose relevant wall and pipe outlet boundary conditions. Perform computations first to establish a steady flow and use that as the base flow for stability studies.
-
(iii) Add a small constant body-force term as described in § 2.2 and obtain the results as presented in the main body of the paper.
We examine each of the above options.
A.1. Option 1: Lamb–Oseen vortex without body force
The artificial source term is removed and the inlet velocity profile is allowed to evolve over time in accordance with the analytical Lamb–Oseen solution. We perform computations for the same two Cases studied in §§ 3.1 and 3.2.
Figure 30. Time history of perturbation energy of the Lamb–Oseen vortex without body force under large disturbances: (a)
$\textit{Re}=5000$
and
$\omega =3.85$
; (b)
$\textit{Re}=700$
and
$\omega =4.3$
.
We start at time
$t_0$
, the inlet Lamb–Oseen vortex evolves with
$b=\textit{Re}/[4(t+t_0)]$
, where
$t_0$
is chosen to ensure that the inlet boundary conditions and the base flow at
$t = 0$
yield an initial vortex core size corresponding to
$b = 4$
. A large initial perturbation is introduced to trigger nonlinear evolution. Figures 30(a) and 30(b) present the time history of the perturbation energy for the two cases: Case 1
$\textit{Re}=5000$
and
$\omega =3.85$
and Case 2
$\textit{Re}=700$
and
$\omega =4.3$
, respectively. The initial growth observed in the early stage indicates that unstable modes allow the perturbations to grow and overcome viscous diffusion. The flow fields shown in figures 31 and 32 closely resemble those in figures 3 and 16, demonstrating that the nonlinear modes captured under the presence of the source term inherently exist in the unforced, naturally evolving flow.
Figure 31. Planar projected perturbation velocity vectors and contours of perturbation pressure for the case
$\textit{Re}=5000$
and
$\omega =3.85$
at
$t=10$
(time-dependent inlet velocity profile): (a)
$x{-}z$
plane; (b)
$\theta -r$
plane at
$x=3$
.
Figure 32. Planar projected perturbation velocity vectors and contours of perturbation pressure for the case
$\textit{Re}=700$
and
$\omega =4.3$
at
$t=2$
(time-dependent inlet velocity profile): (a)
$x{-}z$
plane; (b)
$\theta -r$
plane at
$x=3$
.
As time progresses, the base flow evolves continuously under the influence of viscous diffusion, with the vortex core gradually expanding. This results in a time-dependent variation of the stability characteristics of the base flow. Ultimately, the flow always evolves towards a stable state. The subsequent decrease in perturbation energy shown in figure 30 is due to the instantaneous base flow becoming stable over time.
The above results demonstrate that if the instantaneous profile of the Lamb–Oseen vortex is unstable, perturbations can grow temporarily over a short time before the base vortex weakens under the effect of the viscous diffusion, but eventually will become unobservable as discussed in § 2.2, leading to much difficulty and ambiguity in drawing definitive conclusions on the stability of the original vortex. On the contrary, our present approach of adding a small body-force term in the Navier–Stokes equations to render the Lamb–Oseen flow a clearly defined base flow reveals the same physics in clear and definitive terms.
A.2. Option 2: frozen Lamb–Oseen vortex fixed at inlet
In this case, the flow will evolve in the pipe with the given inlet boundary conditions and it does not satisfy the steady Navier–Stokes equations. One has to wait for the computations to reach a steady-state solution before adding perturbations to it for stability studies. There are two issues with this approach. First, the resulting solution may by itself be unstable, making it impossible to obtain a definitive base flow in spite of the computational effort. Second, the resulting steady base flow may not necessarily be representative of real-world vortices with defined vortex strength and size. At best, it represents a particular flow configuration. Therefore, this approach is not pursued here.
A.3. Option 3: use of a small constant body force to keep the vortex steady – the approach of the present work
The Lamb–Oseen vortex defined by (2.9) has a uniform axial velocity and zero radial velocity. As discussed in § 2.2, it is an idealization that captures the critical features of a real-world vortex with a tight finite vortex core size. Unfortunately, it contains a slight unsteady term because of the expansion of its vortex core due to diffusion.
To enable a meaningful stability study, a small constant body-force term
$\boldsymbol{f}=\omega b^2( r^2e^{-br^2})/\textit{Re}$
is added to the right-hand side of (2.4) in our computations. The constant body-force term is to represent the effect of small local radial and axial flow velocities that counter viscous diffusion to keep the vortex tight in a real-world vortex. It makes the idealized columnar Lamb–Oseen vortex an exact steady solution of the Navier–Stokes equations.
In a real-world vortex, vortex transport due to small local radial and axial velocities balances diffusion to keep the vortex core fixed and tight in its stationary state. The Burgers vortex is an analytical model for it. The Burgers vortex has the same vortical structure (azimuthal velocity profile) as the Lamb–Oseen vortex, but contains variable radial and axial velocities. The inward radial velocity and axial vortex stretching keeps the vortex core fixed. In order to demonstrate that our approach of using the small body-force for this counter-diffusion effect preserves the fundamental stability nature of the problem, we perform computation and stability analysis by using the exact Burgers vortex solution as the base flow without the use of a body force.
The Burgers vortex is an exact steady solution to the Navier–Stokes equations. The non-dimensional velocity components and pressure of the Burgers vortex are (Wu et al. Reference Wu, Ma and Zhou2015)
\begin{equation} u_r=-\frac {1}{2x_0}r, \, u_\theta =\frac {\omega }{4r}\big[1-\exp\big(-br^2\big)\big], \, u_x=\frac {x}{x_0},\, p=\int _0^r \frac {1}{r}u_\theta ^2\ \textrm{d}r, \end{equation}
where
$\omega$
and
$b$
are defined as
\begin{equation} { \omega =\frac {\varGamma _0^*}{2{\pi }U^*R^*}, \,b=\frac {\textit{Re}}{4x_0}. } \end{equation}
Here,
$x_0$
denotes the axial location corresponding to the outlet of the finite-length pipe, where the axial velocity satisfies
$u_x = 1$
. All other scaling parameters are the same as those used for the Lamb–Oseen vortex. Figure 33 compares the growth rate and perturbation energy of the Lamb–Oseen and Burgers vortex at
$\textit{Re}=5000$
and
$\omega =3.85$
, demonstrating that the Burgers vortex is more stable. The normalized axial perturbation velocity profiles in figure 34 and the corresponding velocity vectors in figure 35 together reveal the linear perturbation structure and the final accelerated state, similar to those observed in the Lamb–Oseen vortex results (see figures 3 and 5).
Figure 33. Time history of growth rate and perturbation energy of the Lamb–Oseen vortex (black lines) and Burgers vortex (red lines) for the case
$\textit{Re}=5000$
and
$\omega =3.85$
by DNS.
Figure 34. Normalized profiles of perturbation axial velocity for mode 0 along the axial direction.
Figure 35. Planar projected perturbation velocity vectors and contours of perturbation pressure for the case
$\textit{Re}=5000$
and
$\omega =3.85$
by DNS (Burgers vortex as the base flow): (a)
$x{-}z$
plane; (b)
$\theta -r$
plane at
$x=312.5$
.
Figure 36 presents the perturbation energy evolution of the Burgers vortex at
$\textit{Re}=700$
and
$\omega =4.45$
. After the decay of perturbation energy of mode
$m=0$
in
$0 \leqslant t \leqslant 15$
, the flow exhibits a sequence of stages, including initial linear growth, first nonlinear saturation, secondary unstable growth and nonlinear spiral breakdown that closely resemble the behaviour observed in the Lamb–Oseen vortex at
$\textit{Re} = 700$
and
$\omega = 4.3$
. The normalized axial perturbation velocity profiles shown in figure 37, along with the perturbation velocity vectors of ‘steady’ spiral states in figure 38, show similar dominant spiral perturbation modes when compared with the Lamb–Oseen vortex results (see figures 13 and 16). The transient state is unstable to three-dimensional perturbations and finally develops into a periodic spiral breakdown with oscillating perturbation energy.
Figure 36. Time history of energy of total perturbation velocity and modes
$m=0,\ 1$
for the case
$\textit{Re}=700$
and
$\omega =4.45$
of the Burgers vortex.
Figure 37. Normalized profiles of perturbation axial velocity for mode 1 along the axial direction.
Figure 38. Planar projected perturbation velocity vectors and contours of perturbation pressure for the case
$\textit{Re}=700$
and
$\omega =4.45$
(Burgers vortex as the base flow): (a)
$x{-}z$
plane; (b)
$\theta -r$
plane at
$x=40.75$
.
Our simplified treatment of the Lamb–Oseen with an artificial body force successfully reproduces both the linear instability modes and the nonlinear final states observed in the exact, steady Burgers vortex, indicating that the introduced body force plays a role equivalent to the axial stretching in the Burgers vortex. The results show that the use of an artificial source term is a justified and effective modelling strategy for maintaining a stationary base flow in stability analysis. It enables a clear interpretation of the disturbance growth mechanism while reflecting the physical behaviour of a fixed vortex profile or the instantaneous stability property of a slowly diffusing vortex profile.
Our theoretical model and numerical simulations on a typical base flow, the Lamb–Oseen vortex, reveal general instability mechanisms and capture characteristic linear modes and nonlinear final flow states that are relevant to other vortex flows. The flow parameters, such as swirl numbers, Reynolds number and vortex core size are identified as key parameters that govern the onset and evolution of instability of vortices of specific vortex profiles. The knowledge thus acquired can serve as guidance for further study of a wide range of real swirling flows, such as the Burgers vortex and swirling jets (Billant et al. Reference Billant, Chomaz and Huerre1998).
Open access
