JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

43, 1, pp. 135-155, Warsaw 2005

NUMERICAL INVESTIGATIONS OF THE NATURE OF THE
FIRST BIFURCATION FOR THE FLOW IN AN ANNULAR

ROTOR/STATOR CAVITY

Ewa Tuliszka-Sznitko

Artur Zieliński

Institute of Thermal Engineering, Technical University of Poznań

e-mail: sznitko@sol.put.poznan.pl

Direct numerical simulation is performed to study a transitional flow in
an annual rotating cavity of different aspect ratios L=4.0, 5.0 and diffe-
rent curvature parameters Rm = (R1 +R0)/(R1 −R0) = 1.5−6.0. This
paper reports on the influence of the curvature parameter Rm and end-
wall boundary layers on the nature of the first bifurcation to unsteadiness
and on instability structures in the rotor and stator boundary layers. For
all considered end-wall boundary conditions, we have observed clearly su-
percritical transition to unsteadiness for larger Rm and direct transition
from a steady flow to chaotic one for small Rm. A spectral collocation
method based on the Chebyshev polynomial is used for solving the incom-
pressible Navier-Stokes equations. The time scheme is semi-implicit and
second-order accurate; it corresponds to a combination of the second-order
backward differentiation formula for the viscous diffusion terms and the
Adams-Bashforth scheme for the non-linear terms.Themethod uses a pro-
jection scheme tomaintain the incompressibility constraint.

Key words: rotating cavity, direct method, laminar-turbulent transition

1. Introduction

Flow in a rotating cavity is very important from both theoretical and
practical points of view. Typical industrial configurations are cavities between
compressors and turbine disks. From the theoretical point of view, it is al-
so important to understand the sequence of consecutive bifurcations through
which the initially steady flow undergoes transition to turbulence as the Rey-
nolds number increases. Numerous works have been recently devoted to the



136 E.Tuliszka-Sznitko, A.Zieliński

investigation of instabilities associated with flow in a single rotating disk (Ko-
bayashi et al., 1980; Lingwood 1995, 1996, 1997) and to flow between different
rotating disks (Daube et al., 2001; Itoh, 1991; Tuliszka-Sznitko and Soong,
2000; Tuliszka-Sznitko et al. 2002; Serre et al., 2001, 2004; Serre andPulicani,
2001; Cousin-Rittemard, 1996; Cousin-Rittemard et al., 1998). In the case of
rotor/stator flows, at high rotational rates, the flow consists of two boundary
layers (on the rotating and stationary disk) separated by an inviscid rotating
core. The transition process in both boundary layers is related to type I and
type II generic linear instabilities. Type I instability is due to the presence of
an inflection point in the boundary layer velocity profile. The mechanism for
type II instability is related to combined effects of Coriolis and viscous forces.
Faller (1963), Caldwell and Van Atta (1970) investigated experimentally ty-
pe I and type II instabilities in the Ekman flow and found reasonably good
agreement with the linear stability theory. Savas (1987) studied experimen-
tally unsteady uniformly rotating flow over a stationary disk and found both
rings and spiral structures recognised as type II and type I instability, respec-
tively. In the rotor/stator cavity flow and in flows around a single rotating
disk, experimental results exhibit similar instability structures. However, the
confinement of the geometry (rotor/stator geometry) has effect on the critical
Reynoldsnumber.Large influenceof the end-wall boundaryconditions alsohas
been reported. The influence of the attachment of the end-wall to the rotor or
to the stator has been investigated bymany authors, Adams and Szeri (1982),
Dijkstra andHeijst (1983), Cousin-Rittemard (1996), Cousin-Rittemard et al.
(1998), Oliveira et al. (1991).
In the present paper, we focused attention on the first bifurcation of the

flow in an annular rotor/stator cavity. Calculations are performedusing highly
precise spectral collocationmethod based on theChebyshev polynomials. The
study of the first stage of transition to turbulence using DNS is vital because
this method allows for an accurate description of the instability mechanisms
which are known to play an important role in the breakdown process to tur-
bulence. The nature of the first bifurcation of the flow in a cylindrical cavity
of the aspect ratio L=5 and Rm =1was investigated in detail in the paper
by Serre et al. (2004). These 3D calculations showed, for the first time, that
the transition to unsteadiness for a high aspect ratio (L = 5) in the cylin-
drical cavity is supercritical. Daube et al. (2001) performed 2D calculations
to investigate nature of first bifurcations in annular cavity of the aspect ratio
L=5 and different curvature Rm =1.5−99.0. He found that for larger Rm,
the transition is supercritical and for small Rm ≈ 1.5 is subcritical. In the
present paper, we performed calculations for L = 5 and 4, for different Rm
and for different end-wall boundary conditions to analyse the influence of end-



Numerical investigations of the nature of the first bifurcation... 137

wall conditions on the nature of first bifurcations and instability structures.
We restricted ourselves to a 2D flow to compare results obtained for linear
end-wall profiles to those obtained by Daube et al. (2001).
In the following Sections 2, 3 and 4, we formulated the problemandbriefly

described the numerical method. The results are presented in Section 5, and
concluding remarks are given in Section 6.

2. Geometrical model

The geometrical model is a rotor-stator annular cavity characterized by
the aspect ratio L=(R1−R0)/2h and the curvature parameter Rm =(R1+
R0)/(R1 −R0). The disks are bounded by cylinders of height 2h, which can
be stationary or rotating (Fig.1).

Fig. 1. (a) A scheme of the annular rotating cavity, (b) location of monitoring points
in the radial direction

The rotor rotates with uniform angular velocity Ω = Ωez, where ez de-
notes the unit vector. The origin of the z-axis is located at the mid-height
between the disks. The governing parameters are:

a) the Reynolds number based on the external radius of the disks is
ReR =R

2
1Ω/ν̂

b) the Reynolds number based on the height of the cylinder is
Reh =(2h)

2Ω/ν̂

c) the local Reynolds number based on the viscous scale δ=
√
ν̂/Ω is:

Reδ =
δr∗Ω

ν̂
=

√
r∗2Ω

ν̂

where ν̂ is kinematic viscosity.



138 E.Tuliszka-Sznitko, A.Zieliński

The asterisk denotes a dimensional value. The Reynolds number
ReR = R

2
1Ω/ν̂ is the squared upper bound of the local Reynolds number

Reδ =
√
r∗2Ω/ν̂. The local Reynolds number is used to discuss instability

thresholds and characteristic parameters of instability waves.

3. Mathematical model

The governing equations are 2D Navier-Stokes equations written in
velocity-pressure formulation together with the continuity equation. The equ-
ations are written in the cylindrical polar coordinate system (r,ϕ,z), with
respect to the stationary frame of reference

1

L

∂u

∂r
+

u

L(Rm+r)
+
∂w

∂z
=0

∂u

∂t
+L(Rm+1)Au=−(Rm+1)

∂P

∂r
+
L2(Rm+1)

2

ReR

[
∆u−

u

L2(Rm+r)2

]

(3.1)

∂v

∂t
+L(Rm+1)Av=

L2(Rm+1)
2

ReR

[
∆v−

v

L2(Rm+r)2

]

∂w

∂t
+L(Rm+1)Aw=−L(Rm+1)

∂P

∂z
+
L2(Rm+1)

2

ReR
∆w

where

Au=
1

L
u
∂u

∂r
+w
∂u

∂z
−

v2

L(Rm+r)

Av=
1

L
u
∂v

∂r
+w
∂v

∂z
−

uv

L(Rm+r)
(3.2)

Aw=
1

L
u
∂w

∂r
+w
∂w

∂z

The cylindrical Laplacian operator for a two-dimensional flow is defined in
the following way

∆=
1

L2
∂2

∂r2
+

1

L2(Rm+r)

∂

∂r
+
∂2

∂z2
(3.3)

where t is time, (u,v,w) are velocity V components in the r, ϕ and z di-
rections, respectively, P is pressure. The scales for dimensionless variables of



Numerical investigations of the nature of the first bifurcation... 139

time and velocity are Ω−1 and ΩR1, respectively. The dimensionless axial
co-ordinate is z = z∗/h; z ∈ [−1,1]. The radius co-ordinate is normalised to
obtain the domain [−1,1] requested by the spectralmethodbased on theChe-
byshevpolynomials: r= [2r∗−(R1+R0)]/(R1−R0).Theboundaryconditions
are as follows: no slip boundary conditions at any rigid wall u=w = 0. For
the azimuthal velocity component, the boundary conditions are: v=0 on the
stator and v = (Rm + r)/(Rm +1) on the rotating disk. On the stationary
end-walls, the azimutal velocity equals zero v = 0. However, because of the
singularity of the azimuthal velocity at the junction between the stationary
end-walls and rotating disk, this boundary condition must be modified. The
singularity expresses a physical situation, inwhich there is a thin gap between
the edge of the rotating disk and the stationary end-walls. To eliminate this
singularity, different azimuthal velocity profiles are used for r=±1: the linear
profile v = (1+ z)(Rm + r)/[2(Rm +1)] (Serre and Pulicani, 2001) and the
exponential profile v = [(Rm + r)/(Rm +1)]exp[(z− 1)/0.006] (Serre et al.,
2004).
We proceed as follows: the rotation of the rotor is increased step by step

with a small increment ∆ReR =2000. The solution obtained for smaller ReR
is then used as the initial condition for computation with a higher Reynolds
number. In the first iteration, as the initial condition we use a flowwhich cor-
responds to lack ofmotion in themeridional plane and to the linear azimuthal
velocity profile: u=0, v=(1+r)(z+1)/4, w=0. Computations start with
a sufficiently low Reynolds number, e.g. ReR = R

2
1Ω/ν̂ = 3000 to obtain a

steady flow.The steady flow solutions are expectedwhen the convergence rate
becomes smaller than 10−7

1

∆T
|V n+1−V n| ¬ 10−7 (3.4)

where ∆T is the time step.

4. Direct numerical simulation (DNS)

In this paper, we restrict ourselves to a 2D flow. A numerical solution is
based on the spectral collocation method (Canuto et al., 1988). The Gauss-
Lobatto collocation points are used

ri =cos
iπ

N
zj =cos

jπ

M

i=0,1, . . . ,N
j =0,1, . . . ,M
(ri,rj)∈ [−1,1]× [−1,1]

(4.1)



140 E.Tuliszka-Sznitko, A.Zieliński

where N and M are the number of collocation points in the radial and axial
directions, respectively. The solution Ψ = (u,w,v,P) of equations (3.1) is
approximated bymeans of an expansion of the Chebyshev polynomials in the
r and z directions (Serre et al., 2001; Serre and Pulicani, 2001)

Ψ(r,z)=
N∑

n=0

M∑

m=0

ΨnmTn(r)Tm(z) for r,z ∈ [−1,1] (4.2)

The time scheme is semi-implicit and second-order accurate. It corresponds
to a combination of the second-order backward differentiation formula for the
viscous diffusion term and the Adams-Bashforth scheme for the non-linear
terms. Themethod uses the projection scheme tomaintain the incompressibi-
lity constraint. Details are described in Serre et al. (2001), Serre and Pulicani
(2001).

4.1. Predictor

The first task is to find the pressure distribution. The pressure predictor
Pp is computed from the pressure elliptic equation derived from continuity
equation (3.1)1 and Navier-Stokes equations (3.1)2 -(3.1)4

∆Pp =−DivN(V ) (4.3)

where N(V ) denotes the non-linear terms of Navier-Stokes equations

N(V )= [(N(V ))r,(N(V ))ϕ,(N(V ))z]
>

For the numerical purpose and 2Dflow,wewrite this equation in the following
way

1

L2
∂2Pp

∂r2
+

1

L2(Rm+r)

∂Pp

∂r
+
∂2Pp

∂z2
=

=−2
(1
L

∂

∂r
Aun+

1

L

1

Rm+r
Aun+

∂

∂z
Awn
)
+ (4.4)

+
(1
L

∂

∂r
Aun−1+

1

L

1

Rm+r
Aun−1+

∂

∂z
Awn−1

)

where n, and n−1 denote values obtained from the two previous consecutive
iterations. Equation (4.4) is solved with the following boundary conditions



Numerical investigations of the nature of the first bifurcation... 141

∂Pp

∂r
=−L(2Aun−Aun−1)+

L2(Rm+1)

ReR
·

·
[
2
(
∆un−

un

L2(Rm+r)
2

)
−
(
∆un−1−

un−1

L2(Rm+r)
2

)]

(4.5)

∂Pp

∂z
=−(2Awn−Awn−1)+

L(Rm+1)

ReR
(2∆wn−∆wn−1)

The velocity predictor is calculated from the following equations

∆up−up
1

L2(Rm+r)2
−

ReR
L2(Rm+1)2

3

2

up

∆t
=

ReR
L2(Rm+1)2

·

·
[−4un+un−1

2∆t
+L(Rm+1)(2Au

n−Aun−1)+(Rm+1)
∂Pp

∂r

]

∆vp−vp
1

L2(Rm+r)2
−

ReR
L2(Rm+1)2

3

2

vp

∆t
=

(4.6)

=
ReR

L2(Rm+1)2

[−4vn+vn−1

2∆t
+L(Rm+1)(2Av

n−Avn−1)
]

∆wp−
ReR

L2(Rm+1)2
3

2

wp

∆t
=

ReR
L2(Rm+1)2

·

·
[−4wn+wn−1

2∆t
+L(Rm+1)(2Aw

n−Awn−1)+L(Rm+1)
∂Pp

∂z

]

The boundary conditions for the above equations are as follow

up =un−1 vp = vn−1 wp =wn−1 (4.7)

4.2. Corrector

Corrections of the pressure and velocity are calculated from the following
equations

3

2∆t
(un+1−up)=−(Rm+1)

(∂Pn+1

∂r
−
∂Pp

∂r

)

3

2∆t
(vn+1−vp)= 0 (4.8)

3

2∆t
(wn+1−wp)=−(Rm+1)L

(∂Pn+1

∂z
−
∂Pp

∂z

)

and
1

L

∂un+1

∂r
+

un+1

L(Rm+r)
+
∂wn+1

∂z
=0 (4.9)



142 E.Tuliszka-Sznitko, A.Zieliński

The boundary conditions for equations (4.9) are

V
n+1 ·n=V p ·n (4.10)

We calculate un+1, vn+1,wn+1, Pn+1 by introducing a new value

φ=
2∆t

3
(Pn+1−Pn) (4.11)

to equations (4.8)-(4.9). Finally, we obtain

1

L2
∂2φ

∂r2
+

1

L2(Rm+r)

∂φ

∂r
+
∂2φ

∂z2
=

(4.12)

=
1

L(Rm+1)

[1
L

∂up

∂r
+

up

L(Rm+r)
+
∂wp

∂z

]

The boundary condition for this equation is

Gradφ ·n=0 (4.13)

The corrected pressure and velocity are as follow

Pn+1 =Pp+
3

2∆t
φ un+1 =up− (Rm+1)

∂φ

∂r

vn+1 = vp wn+1 =wp− (Rm+1)L
∂φ

∂z

(4.14)

At every new time level [(n+1)∆t], each flowvariable Ψ(up,vp,wp,Pp,φ)
is a solution to a 2D equation of the following form (Serre andPulicani, 2001)

1

L2
∂2Ψ

∂r2
+

1

L2(Rm+r)

∂Ψ

∂r
+
∂2Ψ

∂z2
−λΨ =S (4.15)

where λ is a constant. The above equation is approximated by the spectral
collocation method, and the right-hand side of equation (4.15) is evaluated
by a standard spectral technique. Then, full diagonalization is used to solve
(4.15) (Haldenwang et al., 1984). For the annular cavity, thematrices of radial
and axial operators are diagonalizable with real eigenvalues. Details of the
technique can be found in Serre and Pulicani (2001) and Haldenwang et al.
(1984).



Numerical investigations of the nature of the first bifurcation... 143

5. Results

The main purpose of the paper is to investigate the influence of end-wall
boundary conditions on the instability of flow and on the nature of the first
bifurcation.

Numerical investigation has been performed for a rotating annular cavity
with aspect ratios L=4 and 5, andwith curvature parameters Rm =1.5−6.
Preliminary tests indicated that the spatial resolution (61×49) in (r,z) for
the considered range of Reynolds numbers constitutes a good compromise be-
tween accuracy and computational cost. The incorporated time step ∆t is
assumed from 1 · 10−3 to 5 · 10−3. Velocity fluctuations are computed with
respect to the average flow solution. The behaviour of dependent variables
has been monitored at 15 points in five different positions in the radial di-
rection N(1/6,1/3,1/2,2/3,5/6) and in three positions in the axial direction
M(9/10,1/2,1/10), where N and M are numbers of collocation points in the
radial and axial direction, respectively. The monitoring points in the radial
direction are marked by letters A, B,M, C and D in Fig.1.

The accuracy of the solution was assessed by considering different grids:
(61× 49), (111× 69) and (121× 91). We have found that solutions are al-
most identical with different meshes. To estimate precision of our results, we
compared them to 2D structures obtained by Serre et al. (2001), Serre and
Pulicani (2001) in his 3D computations (L = 5, Rm = 4, 5, ReR = 62500,
74250) and to 2D results byDaube et al. (2001).We compared some characte-
ristic features of 2D flow structures such as frequency, number of rolls across
the radial extent of the cavity, radial wavelengths and we obtained very good
agreement.

5.1. Basic state

At a small Reynolds number, for instance ReR =3000, the flow is steady
and composed of two disjoint boundary layers (one on each disk) and of the
central core flow. Bachelor showed that a rotating disk drives the fluid be-
low the disk into uniform rigid rotation by means of the viscous effect. This
uniformly rotating core gives then rise to shear layers. Exemplary profiles of
three components of the velocity obtained at themonitoring point A,M and
D are presented in Fig.2a,b,c. We can see that the fluid is pumped radially
outward along the rotating disk and radially inward along the stationary disk.
Both boundary layers are separated by the core rotating with nearly steady
rotation.



144 E.Tuliszka-Sznitko, A.Zieliński

Fig. 2. Profiles of the radial u (a), axial w (b) and azimuthal v (c) velocity
component obtained at the monitoring point A,M and D, ReR =59000,L=5,
Rm =4, exponential end-wall azimuthal velocity profiles, the internal end-wall
attached to the rotor and the outer end-wall attached to the stator, (d) azimuthal
velocity profile obtained at the monitoring point D for ReR =64000,L=5,

Rm =4, exponential end-wall azimuthal velocity profiles, both end-walls attached to
the stator

5.2. Influenceof end-wall conditionson instability structures inboundary
layers (L=5, Rm =4)

In this section, we analyze the influence of the end-wall boundary condi-
tions on the stability structures which appear in both rotating and stationary
disk boundary layers.

As it was pointed out in Section 3, to eliminate singularity at the junctions
between the stationary end-wall and rotatingdisk, the azimutal velocity profile
at the end-wall must be assumed. In our computations, we used two types of
the end-wall profile:

• linear

v=(1+z)
Rm+r

2(Rm +1)
for r=±1 (5.1)

• exponential:

– end-wall attached to the stator

v=
Rm+r

Rm+1
exp
(z−1
0.006

)
for r=±1 (5.2)



Numerical investigations of the nature of the first bifurcation... 145

– end-wall attached to the rotor

v=
Rm+r

Rm+1

[
1− exp

(−z−1
0.006

)]
for r=±1 (5.3)

Figures 3a,b,c show, in the (r∗/h,z) plane, instantaneous iso-lines of azi-
muthal velocity disturbances obtained for ReR = 60000, L= 5, Rm = 4 and
for the following end-wall boundary conditions:

a) linear profiles for the shaft and shroud

b) exponential profiles with both end-walls attached to the stator

c) exponential profiles with the shaft attached to the rotor and with the
shroud attached to the stator.

InFig.3a (case ”a”; linear profiles for the shaft and shroud, ReR =60000)
we observe 2 pairs of vortices in the boundary layer of the rotating disk pro-
pagate along the direction of the base flow (outward). In the boundary layer
of the stationary disk we can see 5 pairs of counter-rotating vortices of the
radial wave number 10 < λ∗r/δ < 29, propagating inward with the phase
speed V ∗φ/Ωr

∗ ≈ −0.1. The radial wavelength of disturbances is defined as
λ∗r = ∆r

∗/nr, where ∆r is the radial length occupied by nr rolls. Figure 3a
shows that the disturbances in the boundary layer of the rotating disk are
weaker than those in the stationary disk. The critical Reynolds number of the
transition to unsteadiness equals ReR =62500.

In case ”b” (both end-walls attached to the stator with exponential azi-
muthal end-wall profiles), we can see that only the boundary layer of the
stationary disk is disturbed and that the stationary internal end-wall does not
transport disturbances from the boundary layer of the stationary disk up to
the rotating disk. As no disturbance can be transported between the layers
from the stationary to rotating disk, the rotating disk remains stable and only
very small oscillations can be observed.
Case ”c” (the shaft attached to the rotor and the shroud attached to the

stator with exponential azimuthal end-wall profiles) is the most unstable. In
the boundary layer of the stationary diskwe observe 6 pairs of vortices propa-
gating inward. Then disturbances are transported up to the boundary layer of
the rotating disk along the rotating end-wall. Figure 2b shows that the axial
velocity component in the monitoring point D for case ”a” is positive and is
significantly different in character from the axial velocity profile in case ”b”
(Fig.2d). Disturbed by vortices coming from the stator, the boundary layer
of the rotating disk is unstable. In the boundary layer of the rotating disk we
observe 2pairs of strong counter-rotating vortices; amplitudes in theboundary



146 E.Tuliszka-Sznitko, A.Zieliński

layers of the stationary and rotating disks are of the same order. Vortices in
the rotating disk are propagated outward towards the outer stationary end-
wall, which transports disturbances down to the stationary disk. The axial
velocity component near the outer end-wall is negative (Fig.2b). The critical
Reynolds number of the first bifurcation equals 48750 and the angular frequ-
ency of disturbances in both boundary layers equals σ =2π/T ≈ 4.2 (where
T is a time period taken from the instability characteristics).

Fig. 3. Time dependent axisymmetric instability at ReR =60000,L=5,Rm =4.
The iso-lines of disturbances of the azimuthal velocity component in the (r∗/h,z)
plane obtained for different end-wall boundary conditions. From the top to the
bottom: (a) linear profiles for the shaft and shroud, (b) exponential profiles with
both end-walls attached to the stator, (c) exponential profiles with the shaft

attached to the rotor and with the shroud attached to the stator

5.3. Nature of the first bifurcation to unsteadiness

Stages of the transition to unsteadiness can by observed be analyzing the
changes of dependent variables in terms of time obtained for consecutive Rey-
nolds numbers. In the present paper, we analyze changes of the axial velocity
component. Let us consider first the classical scenario of the transition. At



Numerical investigations of the nature of the first bifurcation... 147

the beginning of the process, a change in rotation brings about a disturbance
characterized by awave packet which is rapidly damped, and the flow reaches
a steady state (Fig.4a). By increasing ReR, more andmore time is needed to
reach the steady flow and then, over the critical Reynolds number, exponen-
tial oscillatory growth of the disturbances appears (Fig.4b). Finally, for higher
ReR, the flow reaches an asymptotic finite amplitude state. However, for some
flows, a different pattern can be observed: At the beginning of the process, the
wave packet is rapidly damped and the flow reaches a steady state (Fig.6a)
but for higher ReR we do not observe oscillatory disturbances but only direct
transition from a steady to chaotic state (Fig.6c).

Fig. 4. Instability characteristics of the axial velocity component obtained for L=5
and Rm =5 in the monitoring point M, stator boundary layer, (a) ReR =12000,

(b) ReR =64000, (c) ReR =70000

Basing on the 2D results by Daube et al. (2001), obtained for linear end-
wall profiles, we expected direct transition froma steady flow to a chaotic flow
for small Rm. On the other hand, the 3D results published by Serre et al.
(2004) demonstrated that for a cylindrical cavity of the aspect ratio L = 5
(Rm =1.0), a change in the end-wall condition from the linear profile to the
exponential one, resulted in a shift from subcritical to supercritical transition.
In the present section we performed calculations for L=5 and different Rm
(including Rm close to 1) and for two different end-wall boundary conditions,
namely linear and exponential ones in order to check the influence of these
conditions on the nature of the first bifurcation in the annular geometry.



148 E.Tuliszka-Sznitko, A.Zieliński

5.3.1. Cases with exponential azimuthal end-wall velocity profiles

In this section, we consider the most unstable configuration: the internal
end-wall attached to the rotor and the outer end-wall attached to the stator.
The exponential azimuthal velocity end-wall profiles have been used. From
the instability characteristics of the axial velocity component (Fig.4) obta-
ined for L = 5 and Rm = 5 (the monitoring point M, the stator boundary
layer), we can see that, in this particular case, the transition to unsteadiness
is oscillatory. We estimated the critical Reynolds number at ReR = 48500
(σ=2π/T =5.3). Oscillatory solutions are observed in all monitoring points.
Instantaneous iso-lines of disturbances of the azimuthal velocity are presen-
ted in Fig.5. The disturbances in the boundary layer of the stationary disk
are convected downstream towards the internal end-wall where they are lifted
up to the rotating disk. Then, the disturbances are convected towards the
outer end-wall. In the boundary layers of the stationary and rotating disks,
we observed 6 and 2 vortices, respectively.

Fig. 5. Instantaneous iso-lines of disturbances of the azimuthal velocity component
obtained for Rm =5,L=5, ReR =64000, in (r

∗/h,z) plane

A different pattern is observed for Rm = 1.5. At the beginning of the
process (similarly to the case Rm =5.0) disturbances introduced by a change
in the rotation of the disk are very quickly damped and the flow reaches a
steady state (Fig.6). For higher ReR, more andmore time is needed to get a
steady flow. However, over the critical Reynolds number (46500), we observe
direct transition to unsteadiness.We did not observe any oscillatory solution.
This is a hint that the transition to unsteadiness could be subcritical. In Fig.7
we can see that only the boundary layer of the stationary disk is disturbed,
whereas the rotating disk remains undisturbed.Disturbances at the stationary
disk are amplified as they are convected towards the internal end-wall but at
r∗/h = 4 they are rapidly dumped. Very similar behaviour was observed in
the cylindrical cavity (Serre et al., 2004). Due to the existence of a stable area
near the internal end-wall (Fig.7) no disturbances are lifted up to the rotating
disk, and the boundary layer of the rotating disk remains undisturbed.



Numerical investigations of the nature of the first bifurcation... 149

Fig. 6. Instability characteristics of the axial velocity component obtained for L=5
and Rm =1.5 in the monitoring point M, stator boundary layer, (a) ReR =24000,

(b) ReR =39000, (c) ReR =48000, (d) ReR =57000

Fig. 7. Instantaneous iso-lines of disturbances of the azimuthal velocity component
obtained for Rm =1.5,L=5, ReR =57000, in (r

∗/h,z) plane

5.3.2. Cases with linear azimuthal end-wall velocity profiles

In this section we analyse a solution obtained for L=5 and Rm =1.5−6
with linear end-wall azimuthal velocity profiles. These calculations have been
undertaken to compare our results with computations performed by Daube
et al. (2001), who analyzed nature of the first bifurcations for L = 5 and
R0/R1 = 0.99− 0.1 (these correspond to Rm = 199− 1.5). Daube used in
his calculations two methods: stream-function – vorticity formulation, and in
the second version, velocity – pressure formulationwith theChebyshev spatial
approximation (more details can be found in Daube et al., 2001).



150 E.Tuliszka-Sznitko, A.Zieliński

For linear azimuthal velocity end-wall profiles, an oscillatory solution ap-
pears in our calculations only for higher values of Rm. For R0/R1 = 0.4
(Rm =2.333) in all monitoring points in the boundary layer of the stationary
disk we observed that the disturbances are amplified, then they are damped
and again amplified and damped. Finally, however, the flow reaches a stable
state, which is visible in Fig.8. These jumps of amplitudes are also observed
in cases discussed in the previous Section but they were not so rapid. For
Rm = 2.333, the critical Reynolds number of transition to unsteadiness was
established to be about ReR = 55000, which is close to Daube’s result obta-
ined for ReR = 6544. The iso-lines of disturbances of the azimuthal velocity
component obtained for ReR =65000 are presented in Fig.9.We can see that
only the boundary layer of the stationary disk is disturbed, where we obse-
rve 7 pairs of counter-rotating vortices. The angular frequency obtained for
ReR =65000 equals σ=2π/T =5.4.

Fig. 8. An instability characteristic of the axial velocity component obtained in the
monitoring point M in the boundary layer of the stationary disk; Rm =2.333,

L=5, ReR =48000

Fig. 9. Instantaneous iso-lines of disturbances the azimuthal velocity component
obtained for Rm =2.333,L=5, ReR =65000, in (r

∗/h,z) plane

Instability characteristics obtained for R0/R1 = 0.3 (Rm = 1.857) are
slightly chaotic.The criticalReynolds numberof the transition tounsteadiness
was established at ReR = 59000. In the iso-lines of azimuthal disturbances
(Fig.10) we can see dislocations of vortices; two of ten disturbances are out.



Numerical investigations of the nature of the first bifurcation... 151

Thefirst dislocationwe can see at r∗/h=5.4 and the second one at r∗/h=7.
Similar dislocations were observed by Daube et al. (2001). The dislocations
were reported for circular waves by Cousin-Rittemard (1996) in an annular
rotating cavity of L = 5 and Rm = 2 at 39000 < ReR < 56000, and in an
experimental investigation by Schouveiler et al. (1999).

Fig. 10. Instantaneous iso-lines of disturbances of the azimuthal velocity component
obtained for Rm =1.857,L=5, ReR =65000, in (r

∗/h,z) plane

Fig. 11. Instability characteristics of the axial velocity component obtained for
Rm =1.5,L=5, from the top to the bottom; (a) ReR =65000, (b) ReR =68000,

(c) ReR =64000, the initial condition from 68000

For R0/R1 = 0.2 (Rm = 1.5), the critical Reynolds number of transition
to unsteadiness lies somewhere between ReR =65000 and 68000 (see Fig.11).
Figure 11 shows that for ReR = 65000, the flow reaches a steady state at
t≈ 100, and for ReR =68000 the flow is unsteady and apparently chaotic. In



152 E.Tuliszka-Sznitko, A.Zieliński

this case, we did not observe an oscillatory solution but only direct transition
from the steady to chaotic flow. It indicates that the first bifurcation can be
subcritical. In the next step, we checked if some hysteresis could be present.
For thispurpose, the instantaneous statewhichwasobtained for ReR =68000,
was used as the initial condition for the case ReR = 64000. As can be seen
in Fig.11c, the flow reverted to a steady state at t ≈ 120, and therefore, no
hysteresis could be detected in the considered range of the Reynolds number.
This result indicates that transition to unsteadiness is subcritical, which is
in agreement with the experimental results by Schouveiler etal (1999) and
Gauthier et al. (1999).

6. Conclusions

A incompressible fluid flow in an annular rotor/stator configuration with
different aspect ratios L and curvature parameters Rm has been numerically
investigated. The efficient spectral collocation method enabled accurate inve-
stigation of the first stage of transition to turbulent flow.We focused attention
on the influence of the end-wall conditions on instability structures of the flow
and on the nature of the first bifurcation to unsteadiness. Calculations perfor-
med for L=4and Rm =5 showed dramatic influence of different approxima-
tions of end-wall azimutal velocity profiles on the instability structures in both
boundary layers. Themost unstable turns out to be the flowwith exponential
velocity end-wall profiles and with the internal end-wall attached to the rotor
and with the outer one attached to the stator. The spectral parameters of
vortices propagating in the boundary layers have been given.

Calculations performed for L = 5 and different Rm = 1.5− 6.0 with
linear azimuthal end-wall profiles showed good agreement with the results by
Daube et al. (2001) performed for the same geometry. We have found the
classical scenario of transition to unsteadiness (with an oscillatory solution)
for higher Rm. For smaller Rm, however, the first bifurcation seems to be
subcritical. For linear profiles and for smaller Rm, we observed dislocations
on the iso-lines of vortices in the boundary layer of the stationary disk.

Similar series of 2D calculations, carried out for exponential end-wall velo-
city profiles with the outer end-wall attached to the stator andwith the inner
end-wall attached to the rotor, again showed oscillatory transition to unste-
adiness at higher Rm and direct transition from a steady flow to chaotic one
at very small Rm. In the annular cavity, the change of the azimuthal end-wall



Numerical investigations of the nature of the first bifurcation... 153

profiles from linear to exponential did not result in a change of the nature of
the first bifurcation,whichwas observed bySerre et al. (2004)where a 3Dflow
in a cylindrical cavity of the aspect ratio L=5was investigated. Additionally,
our primary calculations of the 3D flow in an annular cavity have shown large
influence of three-dimensionality on the instability structures but rather small
on the nature of the first bifurcation.

Further 3D calculations are necessary to clarify the problem of the nature
of the first bifurcation. The study on this subject and also on the absolu-
te/convective character of the instability in both the rotor/stator (considered
in this paper) and rotor/rotor cavity with throughflow is in progress.

References

1. Adams M.L., Szeri A., 1982, Incompressible flow between finite discs, J.
Appl. Mechanics, 49, 1-9

2. Canuto C., Hussaini M.Y., Qarteroni A., Zang T.A., 1988, Spectral
Methods in Fluids Dynamics, Springer

3. Caldwell D.R., Van Atta C., 1970, Characteristics of Ekman boundary
layer instabilities, J. Fluid Mech., 44, 79-95

4. Cousin-Rittemard N., 1996, Contribution a l’étude des instabilités des
écoulements axisymetriques en cavité inter-disque de type rotor-stator. Docto-
rate Thesis, Université de Paris VI

5. Cousin-Rittemard N., Daube O., Le Quéré P., 1998, Sur la nature de la
premicre bifurcation des écoulements interdisques, C. R. Acad. Sci., serie IIb,
326, 359-366

6. Daube O., Le Quéré P., Cousin-Rittemard M., Jacques R., 2001, In-
fluence of curvature on transition to unsteadiness and chaos of rotor-statordisk
flows, submitted to J. Fluid Mech.

7. DijkstraD.,HeijstG., 1983,The flowbetweenfinite rotatingdiscs enclosed
by a cylinder, J. Fluid Mech., 128, 123-154

8. Faller A.J., 1963, An experimental study of the instability of the laminar
Ekman boundary layer, J. Fluid Mech., 15, 560-576

9. Gauthier G., Gondret P., Rabaud M., 1999, Axisymmetric propagating
vortices in the flow between a stationary and a rotating disk enclosed by a
cylinder, J. Fluid Mech., 386, 105-127



154 E.Tuliszka-Sznitko, A.Zieliński

10. HaldenwangP., LabrosseG.,Abboudi S.,DevilleM., 1984,Chebyshev
3D spectral 2D pseudospectral solver for the Helmoltz equation, J. Comput.
Phys., 55, 115-128

11. Itoh M., 1991, On the instability of the flow between coaxial rotating disks,
In: Boundary Layer Stability and Transition to Turbulence, ASME FED-Vol.
114, 83-89

12. Kobayashi R., Kohama Y., Takamadate Ch., 1980, Spiral vortices in bo-
undary layer transition region on a rotating disk,Acta Mechanica, 35, 71-82

13. Lingwood R., 1995, Absolute instability of the boundary layer on a rotating
disk, J. Fluid. Mech., 299, 17-33

14. LingwoodR., 1996,Experimental studyof absolute instability of the rotating-
disk boundary layer flow, J. Fluid. Mech., 314, 373-405

15. Lingwood R., 1997, Absolute instability of the Ekman layer and related ro-
tating flows, J. Fluid. Mech., 331, 405-428

16. Oliveira L.A., Pecheux J., Restivo A.O., 1991, On the flow between a
rotating and a coaxial fixed disc: numerical validation of the radial similarity
hypothesis,Theoret. Comp. Fluid Dyn., 2, 211-221

17. Savas O., 1987, Stability of Bödewadt flow, J. Fluid Mech., 183, 77-9

18. Schouveiler L., Le Gal P., Chauve M. P., Takeda Y., 1999, Spiral and
circularwaves in the flowbetweena rotating anda stationarydisc,Experiments
in Fluids, 26, 179-187

19. SerreE.,CrespodelArco,BontouxP., 2001,Annular and spiralpatters
in flow between rotating and stationary disks, J. Fluid. Mech., 434, 65-100

20. Serre E., Pulicani J., 2001,A three-dimensional pseudospectralmethod for
rotating flows in a cylinder,Computers and Fluids, 30, 491-519

21. Serre E., Tuliszka-Sznitko, Bontoux P., 2004, Coupled numerical and
theoretical study of the flow transition betweeen a rotating and stationarydisk,
Phys. of Fluids, 16, 3

22. Tuliszka-Sznitko E., Soong C-Y., 2000, Instability of non-isothermal flow
between coaxial rotating disks, Proceedings of European Congress on Compu-
tational Methods in Applied Sciences and Engineering, Barcelona

23. Tuliszka-Sznitko, Serre E., Bontoux P., 2002, On the nature of the
boundary layers instabilities in a flowbetween a rotating and a stationary disc,
C.R. Mecanique, 30, 90-99



Numerical investigations of the nature of the first bifurcation... 155

Badanie metodą symulacji numerycznej charakteru pierwszej bifurkacji
przepływu pomiędzy wirnikiem i stojanem

Streszczenie

Wpracy badana jest stabilność przepływuwwirującej, pierścieniowej przestrzeni
pomiędzywirnikiem i stojanemmetodąbezpośrednią. Zagadnienie przepływuwprze-
strzeniach pomiędzy tarczami wirującymi w różnych konfiguracjach jest istotne nie
tylko ze względówpoznawczych, ale również aplikacyjnych. Tego typu przepływywy-
stępują pomiędzy tarczami silników turbogazowych i sprężarekosiowych.Zagadnienie
może zainteresować w szczególności inżynierów zajmujących się chłodzeniem tarcz
silników turbogazowych. Z punktu widzenia teoretycznego jest bardzo interesujące
poznanie szeregu kolejnych bifurkacji prowadzących w wirujących przestrzeniach od
przepływu laminarnego do turbulentnego. Analizowany jest wpływ bezwymiarowych
parametrów geometrycznych przestrzeni, tj. współczynnika rozciągłości obszaru L,
parametru krzywizny Rm, jak i wpływwarunkówbrzegowych stawianychna pierście-
niach ograniczających przestrzeń na charakter pierwszej bifurkacji. Badane są struk-
tury niestabilnościowewystępujące wwarstwach przyściennychwirnika i stojana. Do
badań zastosowanometodę spektralnej kolokacji bazującą na szeregach Czebyszewa.
Zastosowany schemat po czasie jest kombinacją schematu wstecznego o dokładności
drugiego rzędu i schematu Adamsa-Bashfortha.

Manuscript received June 15, 2004; accepted for print October 18, 2004