2.1 Model equations

An unsteady and axisymmetric flow of an incompressible viscoelastic fluid with swirl is considered in a finite-length straight circular pipe of radius $R$ and length $\bar{x}_{0}$ (see figure 2). The pipe centreline is along the $\bar{x}$ axis, where $\bar{x}$ is the axial distance from the inlet and $0\leqslant \bar{x}\leqslant \bar{x}_{0}$ , $\unicode[STIX]{x1D703}$ is the circumferential angle $(0\leqslant \unicode[STIX]{x1D703}<2\unicode[STIX]{x03C0})$ and $\bar{r}$ is the radial distance from the centreline, where $0\leqslant \bar{r}\leqslant R$ . The radial, circumferential and axial components of the velocity vector $\bar{\boldsymbol{V}}$ are $\bar{u},\bar{v}$ and $\bar{w}$ respectively. The constant density of the fluid is $\bar{\unicode[STIX]{x1D70C}}$ , the pressure gauge from a reference pressure is $\bar{p}$ and the Cauchy stresses are $\bar{\unicode[STIX]{x1D70E}}_{ij}$ , where the indices $i$ and $j$ are $\bar{x}$ , $\unicode[STIX]{x1D703}$ or $\bar{r}$ , and $\bar{\unicode[STIX]{x1D70E}}_{ij}=\bar{\unicode[STIX]{x1D70E}}_{ji}$ (a symmetric stress tensor). The continuity and linear momentum equations that describe the fluid axisymmetric motion are given in a cylindrical coordinate system by

Figure 2.

Schematic diagram of a viscous swirling flow in an open pipe.

The stress tensor $\bar{\unicode[STIX]{x1D70E}}$ of the viscoelastic fluid is described by the generalized Giesekus constitutive model (Bird et al. Reference Bird, Armstrong and Hassager1987a ),

Here, $\text{d}()/\text{d}\bar{t}=\unicode[STIX]{x2202}()/\unicode[STIX]{x2202}\bar{t}+\bar{\boldsymbol{V}}\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D735}}()$ denotes the substantial derivative of a certain property and $\bar{\unicode[STIX]{x1D735}}\bar{\boldsymbol{V}}$ is the second-order tensor of the spatial velocity variations. Moreover, $\bar{\unicode[STIX]{x1D6FE}}_{(1)}$ is the rate-of-strain tensor and $\bar{\unicode[STIX]{x1D70E}}_{(1)ij}$ and $\bar{\unicode[STIX]{x1D6FE}}_{(2)ij}$ denote the contravariant convected derivatives of the extra stress $\bar{\unicode[STIX]{x1D70E}}$ and rate-of-strain $\bar{\unicode[STIX]{x1D6FE}}_{(1)}$ second-order tensors respectively. The viscosity parameter of the fluid is $\unicode[STIX]{x1D702}_{0}=\unicode[STIX]{x1D702}_{s}+\unicode[STIX]{x1D702}_{p}$ , where $\unicode[STIX]{x1D702}_{s}$ is the Newtonian solvent viscosity and $\unicode[STIX]{x1D702}_{p}$ is the added polymer viscosity. The constant $\unicode[STIX]{x1D706}_{1}$ represents the fluid relaxation time, $\unicode[STIX]{x1D706}_{1}\geqslant 0$ , and the constant $\unicode[STIX]{x1D706}_{2}$ represents the retardation time, $\unicode[STIX]{x1D706}_{2}\geqslant 0$ . Furthermore, $\unicode[STIX]{x1D702}_{0}=\unicode[STIX]{x1D702}_{s}(\unicode[STIX]{x1D706}_{1}/\unicode[STIX]{x1D706}_{2})$ . The non-dimensional constant $a$ is the modified mobility parameter. It should be noted that when $\unicode[STIX]{x1D706}_{1}=\unicode[STIX]{x1D706}_{2}=0$ , the fluid is Newtonian. When $\unicode[STIX]{x1D706}_{1}$ is replaced by the Maxwell relaxation time $\unicode[STIX]{x1D706}_{M}=\unicode[STIX]{x1D706}_{1}-\unicode[STIX]{x1D706}_{2}$ and both $\unicode[STIX]{x1D706}_{2}$ and $a$ are set to zero, the fluid is described by the generalized upper convected Maxwell model. When $a=0$ , the fluid is described by the Oldroyd-B model. When $\unicode[STIX]{x1D706}_{1}=\unicode[STIX]{x1D706}_{2}$ , the fluid behaves exactly as Newtonian. When $\unicode[STIX]{x1D706}_{1}=0$ and $\unicode[STIX]{x1D706}_{2}>0$ , the model reduces to the so-called second-order fluid (Bird et al. Reference Bird, Armstrong and Hassager1987a ,Reference Bird, Armstrong and Hassager b ).

In appendix A, we use the Giesekus model to calculate the stresses $\bar{\unicode[STIX]{x1D70E}}_{xx}$ , $\bar{\unicode[STIX]{x1D70E}}_{xy}$ and $\bar{\unicode[STIX]{x1D70E}}_{yy}$ in the case of a two-dimensional unidirectional shear flow between a stationary wall and a parallel infinite plate at distance $h$ moving with velocity $U$ in the $\bar{x}$ direction. The flow is in the $\bar{x}$ direction and the velocity components are $\bar{v}_{x}=\bar{\dot{\unicode[STIX]{x1D6FE}}}\bar{y}$ , $\bar{v}_{y}=\bar{v}_{z}=0$ , where $\bar{\dot{\unicode[STIX]{x1D6FE}}}=U/h$ . We demonstrate that with proper values of $\unicode[STIX]{x1D702}_{0}$ , $\unicode[STIX]{x1D706}_{1}$ , $\unicode[STIX]{x1D706}_{2}$ and $a$ , the Giesekus model may represent the nature of the fluid viscosity $\unicode[STIX]{x1D702}=\bar{\unicode[STIX]{x1D70E}}_{xy}/\bar{\dot{\unicode[STIX]{x1D6FE}}}$ as a function of the shear rate $\bar{\dot{\unicode[STIX]{x1D6FE}}}$ found in the rheological measurements of various Boger fluids shown in figure 12 of Stokes et al. (Reference Stokes, Graham, Lawson and Boger2001).

The governing equations (2.1) and (2.2) are non-dimensionalized by introducing

Here, $U$ is a characteristic speed of the incoming axial velocity to the pipe and all properties are functions of $x$ ( $0\leqslant x\leqslant x_{0}$ , where $x_{0}=\bar{x}_{0}/R$ ) and $r$ ( $0\leqslant r\leqslant 1$ ). Then, the continuity equation and the equations of motion (2.1) become

and the Giesekus constitutive model equation (2.2) becomes

Here, the Reynolds number $Re$ , the Weissenberg number $We$ and the ratio of the fluid characteristic times (the ratio of retardation time to relaxation time) $\unicode[STIX]{x1D706}$ are defined respectively by

The inverse of the Reynolds number, $1/Re$ , gives a measure of the viscous effects with respect to the flow inertia. The Weissenberg number $We$ represents a measure of the fluid elasticity effects. Since $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D702}_{s}/\unicode[STIX]{x1D702}_{0}$ , this parameter provides a measure of shear-thinning effects. The shear thinning is also tuned by the mobility parameter $a$ .

2.2 Boundary conditions

We study the viscous flow in a finite-length rotating straight circular pipe with certain conditions posed on the boundaries; see figure 1. Along the pipe inlet $(x=0)$ , we assume fixed profiles of the axial and circumferential velocity components and of the azimuthal vorticity $\unicode[STIX]{x1D701}=(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x)-(\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}r)$ , given for all $t>0$ by

Here, $\unicode[STIX]{x1D714}$ represents the swirl ratio of the incoming flow. These inlet conditions reflect an incoming flow to the pipe that is created by a vortex generator ahead of the pipe in steady and continuous operation, as found in the experimental apparatuses of Sarpkaya (Reference Sarpkaya1971, Reference Sarpkaya1995), Leibovich (Reference Leibovich1978, Reference Leibovich1984), Malkiel et al. (Reference Malkiel, Cohen, Rusak and Wang1996), Mattner et al. (Reference Mattner, Joubert and Chong2002) and Umeh et al. (Reference Umeh, Rusak, Gutmark, Villalva and Cha2010). We also consider in this paper the special case where $\unicode[STIX]{x1D701}_{0}(r)=-\text{d}w_{0}/\text{d}r$ , which implies a zero axial gradient of the radial velocity at the inlet, i.e. for all $t>0$ , $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x(0,r,t)=0$ . This assumption is for simplicity and does not limit the scope of the analysis, or the results, in any substantial way. It helps to focus the analysis only on the effects of the viscoelastic stresses on the flow dynamics without mixing it with other effects related to non-zero values of $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x(0,r,t)$ . Rusak (Reference Rusak1998) described how general inlet azimuthal vorticity profiles can be included in the study.

In addition, since the stresses are modelled by the differential equations (2.8), the six stresses at the inlet, $\unicode[STIX]{x1D70E}_{ij}(0,r,t)$ , need to be set. We assume that for all $t\geqslant 0$ and $0\leqslant r\leqslant 1$ , the six stresses are given by $\unicode[STIX]{x1D70E}_{ij}(0,r,t)=\unicode[STIX]{x1D70E}_{ij0}(r)$ , where $\unicode[STIX]{x1D70E}_{ij0}(r)$ are the six stresses of a base columnar (parallel) flow state with velocity components $w=w_{0}(r),v=\unicode[STIX]{x1D714}v_{0}(r)$ and $u=0$ . These stresses are computed in § 3.

Along the pipe centreline $(r=0)$ , the symmetry conditions

are used for all $t>0$ .

We assume that the influence of the rotating pipe wall axial velocity boundary layer on the swirling flow in a finite-length pipe is small within the scope of the present analysis, and a no-penetration boundary condition along the rotating wall may be sufficient to reflect the physical situation, i.e. for all $t>0$ , we assume along the wall $(r=1)$

The second condition represents the no slip of the circumferential velocity on the rotating wall. These conditions may reflect the flow behaviour in a pipe with length $x_{0}\ll \sqrt{Re}$ and as long as the axial velocity boundary layer remains sufficiently thin with respect to the pipe radius.

At the pipe outlet $(x=x_{0})$ , we assume for all $t>0$ zero radial velocity and zero axial gradients of the axial and circumferential velocities,

These conditions may reflect the flow development in a pipe where the pipe non-dimensional length satisfies $x_{0}\gg 1$ but is sufficiently shorter than the length for a fully developed viscous flow, $x_{0}\ll \sqrt{Re}$ .

Similar boundary conditions on the flow were considered in the analyses of Wang & Rusak (Reference Wang and Rusak1996, Reference Wang and Rusak1997a ,Reference Wang and Rusak b ) and in the numerical simulations of Beran & Culick (Reference Beran and Culick1992), Beran (Reference Beran1994), Lopez (Reference Lopez1994), Darmofal (Reference Darmofal1996) and Rusak et al. (Reference Rusak, Granata and Wang2015) of swirling flows of a Newtonian fluid in an open pipe. These boundary conditions may also reflect the physical situation as reported in the experiments of Sarpkaya (Reference Sarpkaya1971, Reference Sarpkaya1995), Leibovich (Reference Leibovich1978, Reference Leibovich1984), Malkiel et al. (Reference Malkiel, Cohen, Rusak and Wang1996) and Mattner et al. (Reference Mattner, Joubert and Chong2002).

We emphasize that there is no symmetry between the relatively fixed inlet flow conditions (2.10) and the relatively passive outlet flow conditions (2.13) in the present flow set-up (figure 2). This asymmetry generates an active role of the inlet and outlet flow states in determining the flow structure in the pipe. It forms a dynamic axial inhomogeneity of the convective flow along the pipe which results in the appearance of large axial and radial changes of the flow structure from a parallel (or a near-parallel) flow state, specifically when the swirl is close to the critical level $\unicode[STIX]{x1D714}_{c}$ . The viscous diffusion effects in the bulk induce an additional axial inhomogeneity which increases the size of the flow changes. The viscoelastic flow properties are found in the present study to reduce the viscous stresses and thereby somewhat decrease the viscous diffusion effects.

Equations (2.4)–(2.8) with boundary conditions (2.10)–(2.13) present a well-defined mathematical problem to describe the development of a swirling flow of an incompressible viscoelastic fluid in a pipe for any values of the parameters $\unicode[STIX]{x1D714}$ , $x_{0}$ , $Re$ , $We$ , $\unicode[STIX]{x1D706}$ , $a$ and physically relevant profiles $w_{0}(r)$ and $v_{0}(r)$ of the incoming flow, such as the Q-vortex model (Leibovich Reference Leibovich1984). The general solution of this problem requires extensive numerical computations. In this paper, we concentrate on high-Reynolds-number flows where $Re\gg 1$ , $1\ll x_{0}\ll \sqrt{Re}$ , $0<We\leqslant 2.5$ , and where the swirl ratio of the incoming flow is close to the critical swirl, i.e. $|\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{c}|/\unicode[STIX]{x1D714}_{c}\ll 1$ . In this case, asymptotic analysis can be used to derive approximate solutions of the problem (2.4)–(2.13) which can shed light on the nonlinear flow behaviour, the effect of viscoelasticity parameters on the critical swirl for vortex breakdown, and the flow dynamics and transition to breakdown states in a swirling flow of a viscoelastic fluid.

4.1 The perturbation equations and critical swirl

It is expected that for every inlet swirl level $\unicode[STIX]{x1D714}$ , viscous effects of the order of $O(1/Re)$ create correspondingly small perturbations from the base inviscid-limit columnar flow state where $u=0,w=w_{0}(r),v=\unicode[STIX]{x1D714}v_{0}(r)$ . However, previous studies (Wang & Rusak Reference Wang and Rusak1997b ) have shown that near the critical swirl $\unicode[STIX]{x1D714}_{c}$ the velocity perturbations to a columnar flow solution of order $\unicode[STIX]{x1D716}_{1}$ can be much greater than $O(1/Re)$ . These velocity perturbations are a result of the nonlinear interaction between the near-critical level of swirl of the flow and the small viscous effects, including the viscoelastic effects. Therefore, following Wang & Rusak (Reference Wang and Rusak1997b ), we assume for all $t>0$ and for $0\leqslant x\leqslant x_{0}$ and $0\leqslant r\leqslant 1$ the following asymptotic expansions about the base steady inviscid-limit columnar state:

Here, $K$ denotes the circulation function, with $K_{0}(r)=rv_{0}(r)$ . We assume that $|\unicode[STIX]{x1D716}_{2}|\ll |\unicode[STIX]{x1D716}_{1}|\ll 1$ , $|\unicode[STIX]{x1D716}_{1}|\gg 1/Re$ , and at the second-order terms $\unicode[STIX]{x1D716}_{2}\sim \text{d}\unicode[STIX]{x1D716}_{1}/\text{d}t,\,\sim \unicode[STIX]{x1D716}_{1}^{2}\sim 1/Re$ . (We note that assuming $1/Re\sim \unicode[STIX]{x1D716}_{1}$ results in a singularity and a breakdown of the above asymptotic expansions when $\unicode[STIX]{x1D714}$ approaches the critical level $\unicode[STIX]{x1D714}_{c}$ (Wang & Rusak Reference Wang and Rusak1997b ).) Moreover, it is assumed that $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{c}+\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ , where $|\unicode[STIX]{x0394}\unicode[STIX]{x1D714}|/\unicode[STIX]{x1D714}_{c}\ll 1$ , and that $\unicode[STIX]{x1D716}_{1}$ and $\unicode[STIX]{x1D716}_{2}$ are functions of $t$ , $1/Re$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ . The base-state pressure $p_{0}(r)$ and viscous flow stresses $\unicode[STIX]{x1D70E}_{ij0}(r)$ (where $i,j=r,\unicode[STIX]{x1D703},x$ ) are given by (3.1)–(3.7). It should be noted that when $We\sim O(1)$ , the terms with $\unicode[STIX]{x1D70E}_{ij1}(x,r)$ in (2.5)–(2.7) and in (2.8) may be neglected. These assumptions lead to the richest problem of interaction between inertia, viscosity and elasticity.

The boundary conditions (2.10)–(2.13) result in conditions on the various velocity perturbations,

The continuity equation (2.4) and expansions (4.1) and (4.2) give in the leading orders $O(\unicode[STIX]{x1D716}_{1})$ and $O(\unicode[STIX]{x1D716}_{2})$

respectively. Perturbation stream functions $\unicode[STIX]{x1D713}_{1}(x,r)$ and $\unicode[STIX]{x1D713}_{2}(x,r)$ are defined as

The boundary conditions (4.7)–(4.10) provide conditions for $\unicode[STIX]{x1D713}_{1}$ and $\unicode[STIX]{x1D713}_{2}$ ,

In the leading-order terms of order $O(\unicode[STIX]{x1D716}_{1})$ , the equations of motion (2.5)–(2.7) and expansions (4.1)–(4.6) result in

and in the second-order terms of orders $O(\unicode[STIX]{x1D716}_{2},\text{d}\unicode[STIX]{x1D716}_{1}/\text{d}t,\unicode[STIX]{x1D716}_{1}^{2},\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x0394}\unicode[STIX]{x1D714},1/Re)$ , they result in

Elimination of pressure from (4.19) and (4.21) by cross-differentiation and subtraction leads to the leading-order equation for the azimuthal vorticity,

In addition, (4.13) and (4.20) show that

Integration of (4.26) with respect to $x$ , subject to an inlet condition for $v_{1}$ in (4.7) and for $\unicode[STIX]{x1D713}_{1}$ in (4.15) and (4.16), gives

Let $y=r^{2}/2$ , where $0\leqslant y\leqslant 1/2$ and $\unicode[STIX]{x1D6FA}_{c}=\unicode[STIX]{x1D714}_{c}^{2}$ . Then, substituting (4.27) into (4.25) and using (4.16), followed by some rearrangement and integration with respect to $x$ , results in a partial differential equation for the solution of $\unicode[STIX]{x1D713}_{1}(x,y)$ ,

with boundary conditions

A non-trivial solution of this boundary value problem (4.28) and (4.29) exists only for specific values of $\unicode[STIX]{x1D6FA}_{c}$ , the eigenvalues of the problem. The lowest eigenvalue $\unicode[STIX]{x1D6FA}_{c}=\unicode[STIX]{x1D714}_{c}^{2}$ was defined by Wang & Rusak (Reference Wang and Rusak1996) as the ‘critical swirl for a finite-length pipe’. The corresponding eigenfunction is then given by

where $\unicode[STIX]{x1D6FA}_{c}$ and $\unicode[STIX]{x1D6F7}$ are solved from

For given profiles of $w_{0}$ and $K_{0}$ , the problem (4.31) is solved numerically by using a standard fourth-order Runge–Kutta integration scheme with a scaling condition $\text{d}\unicode[STIX]{x1D6F7}/\text{d}y(0)=1$ (this scaling condition does not affect the results described in § 5 since it scales out of the resulting speeds). The critical swirl level $\unicode[STIX]{x1D714}_{c}$ is a transcritical bifurcation point of the equilibrium solutions of the inviscid-limit flow problem. It should be noted that when the pipe length $x_{0}$ tends to infinity, the critical swirl $\unicode[STIX]{x1D714}_{c}$ and the related eigenfunction approach those of Benjamin (Reference Benjamin1962).

The solution of $\unicode[STIX]{x1D6F7}(y)$ provides the spatial structure of the first-order velocity and pressure perturbations,

The amplitude $\unicode[STIX]{x1D716}_{1}(t)$ of these perturbations is determined in the next section.

4.2 The near-critical flow behaviour

To determine the size of $\unicode[STIX]{x1D716}_{1}(t)$ , we analyse the second-order equations (4.22)–(4.24). Elimination of the pressure perturbation $p_{2}$ by cross-differentiation of (4.22) and (4.24) and subtraction gives

Moreover, from (4.23) we find that

The substitution of (4.34) in (4.33), and the use of (4.13) and (4.14), after some rearrangement and integration with respect to $x$ , yields

with boundary conditions

Equation (4.35) is a non-homogeneous partial differential equation for the solution of $\unicode[STIX]{x1D716}_{2}(t)\unicode[STIX]{x1D713}_{2}(x,y)$ , subject to homogeneous boundary conditions (4.29) and (4.36). Since the related homogeneous problem has the non-trivial solution $\unicode[STIX]{x1D713}_{1c}$ , solvability of (4.35) requires that the forcing term is orthogonal to $\unicode[STIX]{x1D713}_{1c}$ . Multiplying (4.35) by $\unicode[STIX]{x1D713}_{1c}$ , integrating over the domain ( $0\leqslant x\leqslant x_{0},0\leqslant y\leqslant 1/2$ ) and using the homogeneous boundary conditions for $\unicode[STIX]{x1D713}_{2}$ results in the solvability condition,

Here, $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}=2\unicode[STIX]{x1D714}_{c}\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ and

where

It should be noted that to the orders considered here, the stresses according to (3.2)–(3.7) are computed at $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{c}$ , which is a function of the parameters characterizing $v_{0}$ and $w_{0}$ . Moreover, the parameter $P_{3}$ is a function of $We$ , $\unicode[STIX]{x1D706}$ and $a$ . In the case of a Newtonian fluid, where $We=0$ , the base flow stresses are given by $\unicode[STIX]{x1D70E}_{rx0}=\text{d}w_{0}/\text{d}r$ and $\unicode[STIX]{x1D70E}_{r\unicode[STIX]{x1D703}0}=\unicode[STIX]{x1D714}_{c}((\text{d}v_{0}/\text{d}r)-(v_{0}/r))$ . Then, the parameter $P_{3}$ as well as all of the following results reduce to those found in Wang & Rusak (Reference Wang and Rusak1997b ). In addition, it can be seen from (3.2)–(3.7) and (4.39) that, for all values of $We\neq 0$ , when $\unicode[STIX]{x1D706}=1$ (the relaxation and retardation times are equal, $\unicode[STIX]{x1D706}_{1}=\unicode[STIX]{x1D706}_{2}$ ), both base flow stresses $\unicode[STIX]{x1D70E}_{rx0}$ and $\unicode[STIX]{x1D70E}_{r\unicode[STIX]{x1D703}0}$ reduce to their expressions for a Newtonian fluid, and to the orders considered, the fluid behaves like a Newtonian fluid.

We note that in the case of high $Re$ , $We\sim O(1)$ , $0<\unicode[STIX]{x1D706}<1$ and $0<a<1$ , the model equation (4.37) does not include the normal stresses or the normal stress difference since they contribute effects of higher orders, $O(\unicode[STIX]{x1D716}_{1}^{3})$ . Moreover, within the present assumptions, the effects resulting from the unsteady terms in the Giesekus model are also of higher order $O(\unicode[STIX]{x1D716}_{1}^{3})$ and are therefore not included in (4.37). When $Re$ is high, the dynamics of the flow is determined by lumped inertia effects given by $P_{s}$ , $P_{1}$ , $P_{2}$ and the lumped viscoelastic effect $P_{3}$ , similarly to the Newtonian fluid case.

For relevant base rotating flows such as the Rankine vortex, the Lamb–Oseen vortex ( $\unicode[STIX]{x1D6FF}=0$ ) or the Q-vortex (see Leibovich Reference Leibovich1984), it is found that $P_{s},P_{1},P_{2}$ and $P_{3}$ are usually positive for every vortex core radius, $We\sim O(1)$ , $\unicode[STIX]{x1D706}$ and $a$ .

4.3 The steady-state case

The steady-state version of (4.37) (where $\text{d}\unicode[STIX]{x1D716}_{1}/\text{d}t=0$ ) has a real solution of $\unicode[STIX]{x1D716}_{1}$ if and only if

When $|\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}|<16/\sqrt{3\unicode[STIX]{x03C0}^{3}}(\sqrt{P_{1}P_{3}}/P_{2})\sqrt{x_{0}/Re}$ , (4.37) has no real solution, indicating that no near-columnar steady solution exists at inlet swirl levels close to the critical state. The states where $|\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}|=16/\sqrt{3\unicode[STIX]{x03C0}^{3}}(\sqrt{P_{1}P_{3}}/P_{2})\sqrt{x_{0}/Re}$ or where

are special states which are saddle fold bifurcation points of solutions of the steady and axisymmetric viscoelastic fluid flow equations of motion. The limit swirl ratios $\unicode[STIX]{x1D714}_{c1}$ and $\unicode[STIX]{x1D714}_{c2}$ may be considered as the corrections to the inviscid critical swirl as a result of the viscoelastic effects. In the case of a Newtonian fluid, the inlet swirl level $\unicode[STIX]{x1D714}_{c1}$ is identified as the first limit point found in the computations of Beran & Culick (Reference Beran and Culick1992).

We also find that near the critical state and when (4.40) is satisfied, there exist two solutions of $\unicode[STIX]{x1D716}_{1}$ ,

Here, $\unicode[STIX]{x1D716}_{1}^{+}$ corresponds to the $+$ sign in front of the square root term in (4.42) and $\unicode[STIX]{x1D716}_{1}^{-}$ to the $-$ sign in (4.42). Equations (4.1), (4.32) and (4.42) give the asymptotic expansion of the axial velocity $w(x,y)$ of the two possible steady states at inlet swirl levels near the critical swirl $\unicode[STIX]{x1D714}_{c}$ ,

Equation (4.43) illustrates several interesting points. For a high- $Re$ flow, the transcritical bifurcation of solutions of the inviscid-limit equations at $\unicode[STIX]{x1D714}_{c}$ breaks into two branches of solutions of the viscoelastic fluid flow equations with a finite gap between the two branches, the size of which is $\unicode[STIX]{x1D714}_{c2}^{2}-\unicode[STIX]{x1D714}_{c1}^{2}=32/\sqrt{3\unicode[STIX]{x03C0}^{3}}\sqrt{P_{1}P_{3}}/P_{2}\sqrt{x_{0}/Re}$ .

A representative bifurcation diagram of swirling flow behaviour of a viscoelastic fluid in terms of the centreline axial velocity at the pipe outlet, $w(x_{0},0)$ , as a function of the inlet swirl ratio, $\unicode[STIX]{x1D714}$ , is shown in figure 6. Here, the inlet flow is a Q-vortex (3.9) with an axial jet where $\unicode[STIX]{x1D6FD}=4$ , $\unicode[STIX]{x1D6FF}=1$ and the pipe length is $x_{0}=6$ . The axial velocity at the outlet centreline is calculated from $w^{\pm }(x_{0},0)=w_{0}(0)+((\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}P_{2}\pm \sqrt{(\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA})^{2}P_{2}^{2}-256(x_{0}/Re)P_{1}P_{3}/3\unicode[STIX]{x03C0}^{3}})/(16P_{1}/3\unicode[STIX]{x03C0}))$ . Moreover, in figure 6, $Re=1\,000$ , $We=1.5$ , $\unicode[STIX]{x1D706}=0.5$ and $a=0.5$ . For this flow, $\unicode[STIX]{x1D714}_{c}=1.5216$ , $\unicode[STIX]{x1D714}_{c1}=1.3111$ and $\unicode[STIX]{x1D714}_{c2}=1.7320$ . It can be seen that steady near-columnar flow solutions of the viscoelastic fluid flow equations exist for $0<\unicode[STIX]{x1D714}<\unicode[STIX]{x1D714}_{c1}$ and $\unicode[STIX]{x1D714}>\unicode[STIX]{x1D714}_{c2}$ . For the branch with $0<\unicode[STIX]{x1D714}<\unicode[STIX]{x1D714}_{c1}$ , we find that the flow is decelerated along the pipe centreline and the axial velocity reaches a minimum at the pipe outlet. For the branch where $\unicode[STIX]{x1D714}>\unicode[STIX]{x1D714}_{c2}$ , the flow is accelerated along the pipe centreline and the axial velocity reaches a maximum at the pipe outlet. The gap between the two branches of solutions disappears and a transcritical bifurcation dominates the flow behaviour in the inviscid-limit flow as $Re$ tends to infinity (see figure 6).

The meaning of the gap between the branches of solutions around the critical swirl is that no near-columnar axisymmetric equilibrium state can exist in a finite range of incoming swirl around the critical level, $\unicode[STIX]{x1D714}_{c1}<\unicode[STIX]{x1D714}<\unicode[STIX]{x1D714}_{c2}$ . The flow must develop large disturbances in this swirl range. Beyond this gap, when $0<\unicode[STIX]{x1D714}<\unicode[STIX]{x1D714}_{c1}$ or $\unicode[STIX]{x1D714}>\unicode[STIX]{x1D714}_{c2}$ , there is a range of inlet swirl where two equilibrium states may co-exist under the same inlet/outlet conditions. One, given by the positive root in (4.43), describes a near-columnar flow for which $\unicode[STIX]{x1D716}_{1}$ is small in size, and the other, given by the negative root in (4.43), describes a relatively large-amplitude solitary-wave disturbance that either decelerates the axial flow when $0<\unicode[STIX]{x1D714}<\unicode[STIX]{x1D714}_{c1}$ or accelerates the flow when $\unicode[STIX]{x1D714}>\unicode[STIX]{x1D714}_{c2}$ .

Let $\unicode[STIX]{x1D713}(x,y)=\unicode[STIX]{x1D713}_{0}(y)+\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D713}_{1}(x,y)$ be the stream function field of the steady-state flow, where $\unicode[STIX]{x1D713}_{0}(y)=\int _{0}^{y}w_{0}(y^{\prime })\,\text{d}y^{\prime }$ . Figure 7 presents contours of $\unicode[STIX]{x1D713}(x,y)$ for four representative steady states along the branches of the viscoelastic flow solutions shown in figure 6. Figure 7(a) shows the near-columnar state (with a small disturbance) and figure 7(b) describes the large-disturbance state, both found at $\unicode[STIX]{x1D714}=1.26$ . In the near-columnar state, the contours are nearly parallel to the centreline, while in the large-disturbance state, the flow is decelerated around the centreline, as reflected by the diverging contours around the centreline. Figures 7(c) and 7(d) show respectively the near-columnar state and the large-disturbance state found at $\unicode[STIX]{x1D714}=1.8$ . In the large-disturbance state, the flow is accelerated around the centreline, as reflected by the contracting contours around the centreline.

Figure 6.

A representative bifurcation diagram of steady-state solutions of near-critical swirling flows of viscoelastic fluids in a pipe for a Q-vortex (3.9) with $\unicode[STIX]{x1D6FD}=4$ , $\unicode[STIX]{x1D6FF}=1$ and $x_{0}=6$ , where $Re=1000$ , $We=1.5,$ $\unicode[STIX]{x1D706}=0.5$ and $a=0.5$ . Also shown for reference is the inviscid-limit case, where $Re\rightarrow \infty$ .

Figure 7.

Contours of $\unicode[STIX]{x1D713}(x,y)$ of representative flow states at (a) $\unicode[STIX]{x1D714}=1.26$ , the near-columnar state; (b) $\unicode[STIX]{x1D714}=1.26$ , the large-disturbance state; (c) $\unicode[STIX]{x1D714}=1.8$ , the near-columnar state; (d) $\unicode[STIX]{x1D714}=1.8$ , the large-disturbance state. Here, $\unicode[STIX]{x1D6FD}=4$ , $\unicode[STIX]{x1D6FF}=1$ and $x_{0}=6$ . Moreover, $Re=1000$ , $We=1.5,$ $\unicode[STIX]{x1D706}=0.5$ and $a=0.5$ .

The above results also show that the effect of the non-Newtonian fluid properties on the swirling flow behaviour can be lumped into an effective Reynolds number,

Essentially, the effective Reynolds number provides, within the details of the present theory, a simple universal similarity parameter to link between the behaviour of a swirling flow of a non-Newtonian fluid and the respective behaviour of a Newtonian fluid. The use of an effective Reynolds number can be applied to any model for the stresses of a non-Newtonian fluid.

Figure 8 presents three lines: (i) $Re_{eff}/Re$ as a function of $We$ with fixed values of $\unicode[STIX]{x1D706}=0.5$ and $a=0.5$ (the solid line); (ii) $Re_{eff}/Re$ as a function of $\unicode[STIX]{x1D706}$ with fixed values of $We=1.5$ and $a=0.5$ (the dash–dot line); (iii) $Re_{eff}/Re$ as a function of $a$ with fixed $We=1.5$ and $\unicode[STIX]{x1D706}=0.5$ (the dashed line). Here, we use a Q-vortex (3.9) with $\unicode[STIX]{x1D6FD}=4$ , $\unicode[STIX]{x1D6FF}=1$ and $x_{0}=6$ . It is found that $Re_{eff}/Re$ always increases above one with increase of $We$ and $a$ from zero and with decrease of $\unicode[STIX]{x1D706}$ from one. The effective Reynolds number is always above the Reynolds number of the flow of a viscoelastic fluid due to the reduction in size of the shear stresses. The most dominant effect is the significant increase of $Re_{eff}$ on decreasing $\unicode[STIX]{x1D706}$ from one to zero, which is attributed to the increase of corresponding shear-thinning effects. This behaviour causes the increase of the critical swirl $\unicode[STIX]{x1D714}_{c1}$ and the decrease of $\unicode[STIX]{x1D714}_{c2}$ with the increase of $We$ (fluid elasticity) and $a$ (mobility) from zero or with the decrease of $\unicode[STIX]{x1D706}$ from one (increase of shear thinning).

Figure 8.

Plots of $Re_{eff}/Re$ as a function of $We$ with fixed values $\unicode[STIX]{x1D706}=0.5$ and $a=0.5$ (the solid line), $Re_{eff}/Re$ as a function of $\unicode[STIX]{x1D706}$ with fixed $We=1.5$ and $a=0.5$ (the dash–dot line) and $Re_{eff}/Re$ as a function of $a$ with fixed $We=1.5$ and $\unicode[STIX]{x1D706}=0.5$ (the dashed line) for a near-critical swirling flow of viscoelastic fluids in a pipe. Here, we use a Q-vortex (3.9) with $\unicode[STIX]{x1D6FD}=4$ , $\unicode[STIX]{x1D6FF}=1$ and $x_{0}=6$ .

4.4 The stability of steady states

We use the ordinary differential equation (4.37) for the evolution of $\unicode[STIX]{x1D716}_{1}(t)$ to study the nonlinear stability of the steady-state solutions determined in § 4.3 in response to an initial perturbation at $t=0$ given by $\unicode[STIX]{x1D716}_{1}(0)=\unicode[STIX]{x1D716}_{10}$ . When $0<\unicode[STIX]{x1D714}<\unicode[STIX]{x1D714}_{c1}$ or $\unicode[STIX]{x1D714}>\unicode[STIX]{x1D714}_{c2}$ , the solution of (4.37) gives

It should be noted that $\unicode[STIX]{x1D716}_{1}^{+}-\unicode[STIX]{x1D716}_{1}^{-}=\sqrt{(\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA})^{2}P_{2}^{2}-256(x_{0}/Re)P_{1}P_{3}/3\unicode[STIX]{x03C0}^{3}}/(8P_{1}/3\unicode[STIX]{x03C0})>0$ . We find that as $t$ tends to infinity, $\unicode[STIX]{x1D716}_{1}$ tends to $\unicode[STIX]{x1D716}_{1}^{+}$ when $\unicode[STIX]{x1D716}_{10}>\unicode[STIX]{x1D716}_{1}^{-}$ . On the other hand, when $\unicode[STIX]{x1D716}_{10}<\unicode[STIX]{x1D716}_{1}^{-}$ , $\unicode[STIX]{x1D716}_{1}$ evolves to growing negative values within a finite time, indicating the formation of a separation zone around the centreline, on the way to establishing a vortex breakdown state. This means that the steady-state solution with $\unicode[STIX]{x1D716}_{1}^{+}$ is nonlinearly asymptotically stable, while the steady-state solution with $\unicode[STIX]{x1D716}_{1}^{-}$ is unstable and cannot be realized from the long-term flow dynamics.

Figure 9 demonstrates the evolution of $\unicode[STIX]{x1D716}_{1}(t)$ starting from various initial conditions $\unicode[STIX]{x1D716}_{10}$ for the case of a Q-vortex (3.9) with $\unicode[STIX]{x1D6FD}=4$ , $\unicode[STIX]{x1D6FF}=1$ , $\unicode[STIX]{x1D714}=1.26$ and $x_{0}=6$ , where $Re=1000$ , $We=1.5,$ $\unicode[STIX]{x1D706}=0.5$ and $a=0.5$ . In this case, the steady-state solutions are $\unicode[STIX]{x1D716}_{1}^{+}=-0.5207$ (the dash–dot line) and $\unicode[STIX]{x1D716}_{1}^{-}=-2.0440$ (the dashed line). It can be seen that, indeed, starting from all $\unicode[STIX]{x1D716}_{10}>\unicode[STIX]{x1D716}_{1}^{-}$ , the solution of $\unicode[STIX]{x1D716}_{1}(t)$ approaches $\unicode[STIX]{x1D716}_{1}^{+}$ as $t$ tends to infinity. However, starting from all $\unicode[STIX]{x1D716}_{10}<\unicode[STIX]{x1D716}_{1}^{-}$ , the solution of $\unicode[STIX]{x1D716}_{1}(t)$ increases in negative values as $t$ increases.

When $\unicode[STIX]{x1D714}_{c1}<\unicode[STIX]{x1D714}<\unicode[STIX]{x1D714}_{c2}$ , the solution of (4.37) gives

Here, $C$ is a constant of integration and is a function of $\unicode[STIX]{x1D716}_{10}$ . As $t$ approaches a finite value $t^{\ast }$ for which the argument of the tangent function on the right-hand side of (4.46) tends to $-\unicode[STIX]{x03C0}/2$ , $\unicode[STIX]{x1D716}_{1}$ evolves to growing negative values. This indicates again, in the range of swirl $\unicode[STIX]{x1D714}_{c1}<\unicode[STIX]{x1D714}<\unicode[STIX]{x1D714}_{c2}$ , the continuous growth of the decelerating flow perturbation along the centreline in a finite time to form a separation zone around the centreline, on the way to establishing a vortex breakdown state. The present small-disturbance theory does not apply anymore after such a situation is established.

Figure 9.

Solutions of $\unicode[STIX]{x1D716}_{1}$ as a function of time $t$ , starting from various initial conditions $\unicode[STIX]{x1D716}_{10}$ , for a Q-vortex (3.9) with $\unicode[STIX]{x1D6FD}=4$ , $\unicode[STIX]{x1D6FF}=1$ , $\unicode[STIX]{x1D714}=1.26$ and $x_{0}=6$ , where $Re=1000$ , $We=1.5,$ $\unicode[STIX]{x1D706}=0.5$ and $a=0.5$ . Also shown for reference are the steady-state solutions for this case, a dotted line for $\unicode[STIX]{x1D716}_{1}^{+}=-0.5207$ and a dashed line for $\unicode[STIX]{x1D716}_{1}^{-}=-2.0440$ .

The stability analysis shows that the branch of the near-columnar flow solutions is asymptotically stable to small axisymmetric disturbances when $0<\unicode[STIX]{x1D714}<\unicode[STIX]{x1D714}_{c1}$ and is unstable when $\unicode[STIX]{x1D714}>\unicode[STIX]{x1D714}_{c2}$ . Moreover, the branch of large-disturbance solutions that folds backward at $\unicode[STIX]{x1D714}_{c1}$ is unstable to axisymmetric disturbances, while the branch of large-disturbance solutions that folds forward at $\unicode[STIX]{x1D714}_{c2}$ is asymptotically stable. The limit (fold) points $\unicode[STIX]{x1D714}_{c1}$ and $\unicode[STIX]{x1D714}_{c2}$ found in the present analysis are also points of exchange of stability of the respective branches of steady-state solutions. This is indicated in figure 6 by the solid lines which represent the stable states and the dashed lines which represent the unstable states. As $Re$ tends to infinity, the linear stability results of Wang & Rusak (Reference Wang and Rusak1996, Reference Wang and Rusak1997a ) for the inviscid-limit solutions are recovered. Moreover, the flow simulations of Beran (Reference Beran1994) and Lopez (Reference Lopez1994), based on the unsteady and axisymmetric Navier–Stokes equations, demonstrate the same stability characteristics of near-critical swirling flows for a Newtonian fluid.