Numerical study of Taylor–Couette flow with an axial flow

Computers & Fluids 33 (2004) 97–118 www.elsevier.com/locate/compfluid

Numerical study of Taylor–Couette flow with an axial flow Jong-Yeon Hwang, Kyung-Soo Yang

*

Department of Mechanical Engineering, Inha University, Incheon 402-020, Republic of Korea Received 20 April 2002; received in revised form 30 October 2002; accepted 23 January 2003

Abstract The flow between two concentric cylinders with the inner one rotating and with an imposed pressuredriven axial flow is studied using numerical simulation. This study considers the identical flow geometry and flow parameters as in the experiments of Wereley and Lueptow [Phys. Fluids 11 (12) (1999) 3637], where particle image velocimetry measurements were carried out to obtain detailed velocity fields in a meridional plane of the annulus. The objectives of this investigation are to numerically verify the experimental results of Wereley and Lueptow and to further study detailed flow fields and bifurcations related to Taylor–Couette flow with an imposed axial flow. The vortices in various flow regimes such as non-wavy laminar vortex, wavy vortex, non-wavy helical vortex, helical wavy vortex and random wavy vortex are all consistently reproduced with their experiments. It is demonstrated that Ôshift-and-reflectÕ symmetry holds in Taylor–Couette flow without an imposed axial flow. In case of Taylor–Couette flow with an imposed axial flow, one can find that the shift-and-reflect symmetry is roughly valid for the remaining velocity field after subtracting the annular Poiseuille flow. The axial flow stabilizes the flow field and decreases the torque required by rotating the inner cylinder at a given speed. Growth rate of the flow instability is defined and used in predicting the type of the vortices. The velocity vector fields obtained also reveal the same vortex characteristics as found in the experiments of Wereley and Lueptow. Ó 2003 Elsevier Ltd. All rights reserved. Keywords: Taylor–Couette flow; Axial flow; Instability; Simulation; Vortex; Shift-and-reflect symmetry

1. Introduction The flow between two concentric cylinders with the inner one rotating and the outer one stationary, called Taylor–Couette flow, has been studied by many researchers for decades. With a low rotating speed of the inner cylinder, the exact solution of the laminar velocity field consists of *

Corresponding author. Tel.: +82-32-860-7322; fax: +82-32-863-3997. E-mail address: [email protected] (K.-S. Yang).

0045-7930/04/$ - see front matter Ó 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0045-7930(03)00033-1

98

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

vr ¼ vz ¼ 0;

vh ¼ Xi ri

ro =r  r=ro ; ro =ri  ri =ro

ð1Þ

where r, h and z represent the radial, azimuthal, and axial directions of the cylindrical coordinate system, respectively. The angular velocity of the inner cylinder is denoted by Xi ; the inner and outer radii are represented by ri and ro , respectively. It should be noted that Eq. (1) is the solution under the idealization of infinitely long cylinders where endwall effects are ignored; this solution is not realized exactly in an physical experiment or simulation that include endwall effects. If the Taylor number (Ta) based on Xi goes over a critical one (Tac1 ), the flow instability caused by the curved streamlines of the main flow produces axisymmetric Taylor vortices. This fact was first noticed by Taylor (1923) in an analytical study of the related flow instability [1]. Since then, many researchers have studied the instability causing Taylor vortices [2–4]. In the early days of studying Taylor vortices, researchersÕ attention was mainly focused on determining Tac1 by experimental or analytical methods. As Ta further increases over a higher threshold value (Tac2 ), the Taylor vortices become unsteady and non-axisymmetric, called wavy vortices [5]. Davey et al. [5] analytically determined the value of Tac2 ; that was subsequently confirmed by EaglesÕ experiment [6]. Many measurements using laser doppler velocimetry (LDV) were also carried out [7,8]. However, they could not obtain detailed flow information but only local velocity fields. Even though some researchers [9,10] have performed numerical investigation of Taylor and wavy vortices, their attention was paid mainly to numerical methods. However, Marcus [11] investigated wavy vortices without an axial flow to get velocity fields and characteristics of traveling waves in detail. In the case of Taylor vortices with a fully developed axial flow, Chung and Astill [12] developed a linear stability theory for the spiral flow. Gravas and Martin [13] studied Taylor vortices in the presence of an axial flow with various ratios of the inner radius to the outer one. They reported that the axial flow stabilizes the flow field and Tac1 is increased. Takeuchi and Jankowski [14] made both experimental and numerical investigations into toroidal vortex structures of the spiral flow; especially flows between two independently rotating concentric cylinders with wide gap were considered. Meseguer and Marques [15] studied the shear instability associated with the axial flow and the centrifugal instability due to rotation using linear stability theory. Recently, measurements of velocity fields for Taylor–Couette flows with and without an imposed axial flow were performed by Wereley and Lueptow [16,17] using particle image velocimetry (PIV). Taylor–Couette flow is often observed in various types of engineering application, for example, journal bearing lubrication and cooling of rotating machinery among others. Since the axial flow not only stabilizes the flow field but also decreases the torque required by rotating the inner cylinder at a given speed as shown in the later part of this paper, understanding Taylor–Couette flow with an imposed axial flow allows us to control the related flow fields. In this study, we investigate the characteristics of the vortices in Taylor–Couette flow having the identical cross-sectional geometry as in the experiments of Wereley and Lueptow [16,17]. Technical limitations of experimental methods keep researchers from getting detailed flow information. They have been limited to obtain only temporally and spatially averaged torque on the inner cylinder or unsteady flow characteristics on a meridional plane in the flow field. Numerical simulations are suitable for efficiently studying effects of flow geometry and flow parameters, and can make up for the shortcomings of experimental techniques by providing detailed information on full three-dimensional (3D) time-dependent flow fields. The objectives of this investigation are

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

99

to numerically verify the experimental results of Wereley and Lueptow [16,17] and to further study detailed flow fields and bifurcations related to Taylor–Couette flow with an imposed axial flow.

2. Formulation The governing incompressible continuity and momentum equations are r  u ¼ 0;

ð2Þ

ou 1 þ ðu  rÞu ¼  rp þ mr2 u; ot q

ð3Þ

where u, q, p and m denote velocity, density, pressure, and kinematic viscosity of the fluid, respectively. The governing equations are discretized using a finite-volume method in a generalized coordinate system. Spatial discretization is second-order accurate. A hybrid scheme is used for time advancement; non-linear terms and cross diffusion terms are explicitly advanced by a thirdorder Runge–Kutta scheme, and the other terms are implicitly advanced by the Crank–Nicolson scheme. A fractional step method [18] is employed to decouple the continuity and momentum equations. The resulting PoissonÕs equation is solved by a multigrid method. For the details of the numerical algorithm used in the code, see [18].

3. Choice of parameters and boundary conditions There are two non-dimensional parameters of importance, namely, Ta and Re. Although Ta can be defined in several different ways, we choose Ta ¼ ri Xi d=m, where d is the difference between ro and ri . The axial Reynolds number is defined as Re ¼ wd=m, where w is the averaged axial velocity. Radius ratio, g, is selected 0.83 as in the experiments of Wereley and Lueptow [16,17]. A special care is required in selecting the size of the axial domain (H). A considerably large axial domain, H ¼ 27d–32d depending on ðTa; ReÞ, is employed for natural selection of the dominant axial wavelength of the vortex pairs (k). In each case of ðTa; ReÞ, we first estimate k using Fig. 11 in Ref. [17], and take 16 multiples of the estimated k as the size of H , ‘‘anticipating’’ appearance of 16 pairs of vortices. It turns out that exactly 16 pairs of vortices show up in all cases reported here. Since the dominant axial wavelength is not imposed a priori but naturally selected, it is estimated that the axial computational domain allows at worst 3.13% error in k. In fact, the error should be much smaller than our ‘‘conservative’’ estimate because the axial domain is ‘‘tuned’’ according to the experimental measurements of k. Temporal instability is considered in association with the axial flow; a periodic boundary condition is employed in the axial direction. For convective instability related to the spiral flow, see [14,15]. Selecting the size of the azimuthal domain needs a special care, too. Once wavy vortices are formed, traveling waves show up. In order to capture the longest wavelength spanning the whole azimuthal domain, if any, we have to consider the full azimuthal domain. Fig. 1 shows the computational grid system employed in this study. Fig. 1(a) and (b) exhibit the computational

100

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

Fig. 1. Computational domain and grid system; (a) cylinder cross-section, (b) meridional section.

domain at one cylinder cross-section and at one meridional section, respectively. The grid is a body-fitted O-grid system which is the most adequate for this flow; it has more resolution near the cylinder walls where gradients are steep. The number of grid points determined by grid-refinement study was 128  32  256 in the azimuthal, radial, and axial directions, respectively. No significant difference in the flow fields is noticed in comparison with those obtained by a 256  64  512 grid system; qualitative difference is hardly visible and the difference in growth rate of the primary instability is less than 1.42% between the two. No-slip boundary condition is imposed on the solid boundaries. An initial flow field for the case of a particular ðTa; ReÞ is constructed using the exact steady solution corresponding to the particular ðTa; ReÞ which would be present at the same ðTa; ReÞ if there were no instability present. In this study, 12 cases with various values of ðTa; ReÞ are computed. For each case, about 15,000 time steps are computed before the flow field is taken for flow analysis; the time span is corresponding to dimensionless time of 3500–6500 based on d and ri Xi or 114–212 revolutions of the inner cylinder depending on ðTa; ReÞ. The vortices appear rather quickly between ri Xi t=d ¼ 150 and 600 depending on ðTa; ReÞ. 4. Results and discussion 4.1. Taylor–Couette flow without an axial flow Wereley and Lueptow [16] reported that the value of Tac1 over which Taylor vortices are formed is 102 and that the value of Tac2 is in the range of 124 < Tac2 < 131. Wereley and LueptowÕs result

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

101

Fig. 2. Velocity vectors and magnitude contours of azimuthal velocity on an r–z plane without an axial flow, contour increment: 0:1Xri ; (a) Wereley and LueptowÕs experiment [16], Ta ¼ 124, (b) current study, Ta ¼ 123. The lower and upper lines represent the inner and the outer cylinder surfaces, respectively.

in Fig. 2(a) shows the velocity vector field on an r–z plane at Ta ¼ 124; Fig. 2(b) shows the result of present numerical study at Ta ¼ 123. The solid curves denote the magnitude contours of azimuthal velocity, while the lower and upper lines represent the inner and the outer cylinder surfaces, respectively. One can notice in both figures that on the boundaries of Taylor vortices the magnitude of the outer directional velocity is larger than that of the inner directional velocity, and that the two adjacent vortices which have outer directional velocity at their boundary (hereafter called ‘‘Vortex Pair’’) get closer to each other. Distributions of the radial and axial velocity components along the axial and radial directions on an azimuthal plane are compared with the corresponding results of Wereley and LueptowÕs experiment in Fig. 3(a) and (b), respectively. Each velocity is non-dimensionalized by the surface velocity of the inner cylinder. Radial location is denoted by n ¼ ðr  ri Þ=d, and axial location by f ¼ z=d. One can notice that the shape of each peak of the experimental result in Fig. 3(a) is different from one another; this could be explained in part by the ‘‘end effects’’ of the cylinder in the experimental setup. Since the axial direction is assumed homogeneous in our simulations, the shape of each peak is identical in the computational result. The peak corresponding to the outflow is narrower and stronger than the one corresponding to the inflow; that is the case not only in the simulation but also in the experiment. This is due to the fact that the distance between the two vortex centers of a Vortex Pair is shorter than that of two adjacent Vortex Pairs, and the outflow between the vortex centers of a Vortex Pair is stronger than the inflow between two adjacent Vortex Pairs as explained in Fig. 2. In Fig. 3(b), computation slightly overpredicts vz along the radial line; this situation may be due to the fact that experimentally measured value of vz at a given point is actually smaller than the corresponding real value of vz . This comes about because the measurement actually reflects the velocity at a region around the measurement point called ‘‘interrogation region’’, not just the measurement point. As a result, at a maximum, the value is slightly reduced due to spatial averaging because the nearby values are slightly less than the maximum [19]. One can also notice that

102

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

present Wereley and Lueptow [16]

0.05

vr /r i Ωi 0

-0.05 1

(a)

2

3

ζ

0.05

present Wereley and Lueptow [16]

vz /r i Ωi 0

-0.05

(b)

0.25

0.5

ξ

0.75

Fig. 3. (a) Distribution of the radial velocity component along an axial line through the center of the annular gap, (b) distribution of the axial velocity component along a radial line through a Taylor vortex center; n ¼ ðr  ri Þ=d, f ¼ z=d, without an axial flow, Ta ¼ 123.

the numerical results can accurately predict the peak locations and the tendency of higher magnitudes near the inner cylinder. Wavy vortices exhibit oscillations in the azimuthal direction as well as in the axial direction. This fact can be explained by a traveling wave. On a given r–z plane, oscillation is noticed as a traveling wave passes through that plane. Therefore, time-dependent characteristics on that plane must be consistent with those along the azimuthal direction. Fig. 4(a) shows time-dependent results of Wereley and Leuptow [16] on a typical r–z plane at Ta ¼ 131, where each frame is taken at a fixed azimuthal location during passage of one wavelength of the traveling wave. Fig. 4(b) shows our results at Ta ¼ 129, where each frame is taken along the azimuthal direction at a given time and at the same phase as the corresponding frame of Fig. 4(a). The symbol ( ) denotes the centers of the vortices. In order to show the computed flow field more clearly, velocity vectors were slightly magnified in Fig. 4(b). The two figures are consistent with each other as shown, confirming the traveling wave found in Wereley and LueptowÕs experiment. The averaged oscillation distance of the centers of the vortices (1:1d) in Fig. 4(b) is slightly longer than that (0:9d) in Fig. 4(a). In addition, the vortex centers oscillate slightly in the radial direction consistent with experimental results [16]. If we compare each frame of Fig. 4(b) with that of half period later, a symmetry rule can be identified. For instance, comparing the second frame with the sixth one in Fig. 4(b), one can notice that vortex 1 and vortex 2 in the second frame are shifted to the counterparts in the sixth

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

103

Fig. 4. Velocity vectors on r–z planes; (a) Wereley and LueptowÕs experiment [16]: Each frame is taken on a fixed azimuthal plane during passage of one wavelength of the traveling wave (from top to bottom), Ta ¼ 131. (b) Current study: at a given instant, each frame is taken along the azimuthal direction at the same phase as the corresponding frame of Fig. 4(a), Ta ¼ 129. The lower and upper lines of each frame represent the inner and the outer cylinder surfaces, respectively.

frame and that vortex 2 in the sixth frame is a reflection of vortex 1 in the second frame with respect to the boundary between vortex 1 and vortex 2. This is called Ôshift-and-reflectÕ symmetry which can be expressed as follows; vr ðr; h; z; tÞ ¼ vr ðr; h þ Dh; z; tÞ;

ð4Þ

vh ðr; h; z; tÞ ¼ vh ðr; h þ Dh; z; tÞ;

ð5Þ

vz ðr; h; z; tÞ ¼ vz ðr; h þ Dh; z; tÞ;

ð6Þ

where z ¼ 0 corresponds to a vortex boundary, and Dh represents the angle corresponding to one half of azimuthal wavelength. Marcus [11] claimed that shift-and-reflect symmetry is adequate to

104

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118 1

ξ

0.8 0.6 0.4 0.2

(a) 0 1.5

2

2.5

ζ

3

1

ξ

0.8 0.6 0.4 0.2

(b)

0

1.5

2

ζ 2.5

Fig. 5. Velocity vectors on an r–z plane, n ¼ ðr  ri Þ=d, f ¼ z=d; (a) Reflected velocity vectors of vortex 1 and vortex 2 in the second frame of Fig. 4(b) with respect to their boundary, (b) velocity vectors corresponding to vortex 1 and vortex 2 in the sixth frame of Fig. 4(b). Solid lines represent magnitude contours of azimuthal velocity, contour increment: 0:14Xri .

be implemented in numerical simulation. Fig. 5(a) shows vortex 1 and vortex 2 in the second frame of Fig. 4(b) after reflection with respect to their boundary and Fig. 5(b) shows the corresponding vortices in the sixth frame of Fig. 4(b). Fig. 5 confirms that shift-and-reflect symmetry is correct, based on our calculation in which such a symmetry was not imposed. In order to quantify the ‘‘degree’’ of shift-and-reflect symmetry, one can define a root-mean-squared ratio (c) as follows; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jua  ub j2 ð7Þ c ¼ qffiffiffiffiffiffiffiffiffi : 2 jub j Here, ua and ub represent two velocity fields on r–z planes which are one half of azimuthal wavelength apart, and ua is the ‘‘reflected’’ field. For example, ua and ub are corresponding to the velocity fields in Fig. 5(a) and (b), respectively. The overline denotes averaging over all the grid points on an r–z plane containing a pair of vortices. In Fig. 5, it turns out that c ¼ 0:00262. 4.2. Taylor–Couette flow with an axial flow Since the axial flow stabilizes the flow field, both Tac1 and Tac2 increase in the presence of the axial flow [17]. Fig. 6(a) shows the velocity field at an r–z plane at Ta ¼ 123 and Re ¼ 4:9, corresponding to non-wavy Taylor vortices. One can find that the axial flow makes the fluid flow up and down alternately around the vortices. Wereley and Lueptow [17] called this a ÔwindingÕ flow.

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

105

1

ξ 0.5

0

(a)

1

2

ζ

3

1

2

ζ

3

1

ξ 0.5

0

(b)

Fig. 6. Velocity vectors and magnitude contours of azimuthal velocity on an r–z plane at Ta ¼ 123 and Re ¼ 4:9, n ¼ ðr  ri Þ=d, f ¼ z=d, contour increment: 0:11Xri ; (a) including the axial velocity profile, (b) with the axial velocity profile removed.

The vortices shrink due to the axial flow and the vortices located toward the inner cylinder is slightly larger than those toward the outer cylinder. Fig. 6(b) reveals the remaining velocity field after subtracting the exact laminar axial velocity profile of annular Poiseuille flow at Re ¼ 4:9, wðrÞ ¼ K

ðro2  ri2 Þ lnðro =rÞ  Kðr2  ro2 Þ; lnðri =ro Þ

ð8Þ

where K ¼ ð1=4lÞðdp=dzÞ and l is the viscosity coefficient of the fluid. Fig. 6(b) is very similar to Fig. 2. Thus one can speculate that the Taylor–Couette flow with an imposed axial flow (Fig. 6(a)) is nearly a linear superposition of non-wavy Taylor vortex flow at Ta ¼ 123 and annular Poiseuille flow at Re ¼ 4:9, as reported in [17]. In order to make an estimate of ‘‘legitimacy’’ of the speculation, one can substract wðrÞ from the axial velocity component at Ta ¼ 123 and Re ¼ 4:9, and compare the remainder with that of Ta ¼ 123 and Re ¼ 0. In Fig. 7, the radial distributions of the two are presented in the same format as Fig. 3(b). Even though there exists a slight quantitative difference in the ‘‘peak’’ magnitudes, it is not difficult to expect a qualitative resemblance between the two flow fields. The reduction in ‘‘peak’’ magnitudes implies reduction of circulation around a Taylor vortex; this is consistent with the notion that axial flow stabilizes Taylor–Couette flow. Fig. 8 shows the velocity field on an r–z plane at Ta ¼ 129 and Re ¼ 13:1, corresponding to non-wavy helical vortices. Vortices are hidden by the relatively strong axial flow; no vortices are observed in Fig. 8(a). Helical vortices do not occur without an axial flow. To visualize ‘‘helical-type’’ structure of Taylor vortices, Fig. 9 presents instantaneous contours of azimuthal velocity component on a portion of the center surface of the annulus for various types of Taylor vortices for comparison. The regions of higher (lower) azimuthal velocity than the median is marked by black (white). Fig. 9(a) represents the characteristics of ‘‘laminar vortex’’ (LV) at Ta ¼ 123 and Re ¼ 0. Traveling waves are not seen and the vortices are stationary. Traveling waves are observed in Fig. 9(b) for

106

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

Ta=123, Re=0 Ta=123, Re=4.9

0.05

vz /r i Ωi 0

-0.05

0

0.25

0.5

ξ

0.75

1

Fig. 7. Distribution of the axial velocity component along a radial line through a Taylor vortex center; n ¼ ðr  ri Þ=d, without an axial flow (Ta ¼ 123 and Re ¼ 0) and with the exact laminar axial velocity profile substracted (Ta ¼ 123 and Re ¼ 4:9).

Fig. 8. Velocity vectors and magnitude contours of azimuthal velocity on an r–z plane at Ta ¼ 129 and Re ¼ 13:1, n ¼ ðr  ri Þ=d, f ¼ z=d, contour increment: 0.11Xri ; (a) including the axial velocity profile, (b) with the axial velocity profile removed.

the case of Ta ¼ 139 and Re ¼ 4:9; this type of Taylor vortex is called ‘‘wavy vortex’’ (WV). Two wavelengths are seen in Fig. 9(b). Fig. 9(c) shows the features of ‘‘helical vortex’’ (HV) appearing at Ta ¼ 129 and Re ¼ 24. Traveling waves do not occur in this case, and the helix angle is approximately 2.9°. Traveling waves do appear at a higher Ta than HV; two wavelengths are contained in Fig. 9(d) for the case of Ta ¼ 167 and Re ¼ 13:1. This type of Taylor vortex is called ‘‘helical wavy vortex’’ (HWV). Fig. 9(e) reveals the features of ‘‘random wavy vortex’’ (RWV) at Ta ¼ 215 and Re ¼ 24. Approximately four wavelengths of dominant traveling waves are identified in the figure. It turned out, however, that unlike the other types of Taylor vortices, circulation around an RWV significantly varies in h within one wavelength as explained later. Nevertheless, RWV maintains quasi-periodicity both in h and in z.

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

107

Fig. 9. Instantaneous contours of azimuthal velocity component on a portion of the center surface of the annulus; (a) laminar vortex (LV) at Ta ¼ 123 and Re ¼ 0, (b) wavy vortex (WV) at Ta ¼ 139 and Re ¼ 4:9, (c) helical vortex (HV) at Ta ¼ 129 and Re ¼ 24, (d) helical wavy vortex (HWV) at Ta ¼ 167 and Re ¼ 13:2, (e) random wavy vortex (RWV) at Ta ¼ 215 and Re ¼ 24. The regions of higher (lower) azimuthal velocity than the median is marked by black (white).

More detailed description on the ‘‘helical’’ structure is given in Figs. 10 and 11 for the case of Ta ¼ 129 and Re ¼ 24 (HV). Fig. 10 shows instantaneous velocity vectors of helical vortices on a portion of an r–z plane in the flow domain along with arrows pointing the corresponding vortex on the other side. Solid and dotted lines represent foreground and background pointers, respectively. The helical vortices form a negative angle with respect to the axis of rotation, and they are shifted by a vortex pair after one revolution. Three-dimensional vortical structures of HV are presented in Fig. 11 for the case of Ta ¼ 129 and Re ¼ 24. Fig. 11(a) shows contours of azimuthal component of vorticity (xh ) in a portion of the flow domain; regions of positive (negative) xh are denoted by dark (bright) color. Nevertheless, the helical structure is clearly seen. See the arrows which point corresponding vortices. To highlight vortex cores only, positive values of J  0:25D2 are plotted in Fig. 11(b) where J represents Jacobian of shear-stress vectors and D denotes divergence of shear-stresses. This method is based on the topological consideration that a vortex core constitutes a focus which satisfies J > 0:25D2 [20]. Vortex cores are clearly visualized by this method (Fig. 11(b)), even though their rotational direction cannot be distinguished.

108

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118 O.C.

I.C.

axis of rotation

I.C.

O.C.

Fig. 10. Instantaneous velocity vectors of helical vortices on a portion of an r–z plane in the flow domain along with arrows pointing the corresponding vortex on the other side; Ta ¼ 129 and Re ¼ 24 (HV). Solid and dotted lines represent foreground and background pointers, respectively. OC and IC denote Ôouter cylinderÕ and Ôinner cylinderÕ, respectively.

Fig. 11(a) and (b) are complementary to each other; one can grasp 3D structures of helical vortices using those two figures. A Ta–Re map is presented in Fig. 12, where solid lines, taken from [17] denote the approximate boundaries of various types of vortices in Taylor–Couette flow. As Re increases, so does Tac2 . For example, the four data points corresponding to Ta ¼ 139 show that Hopf bifurcation to wavy vortices occurs between Re ¼ 13:1 and Re ¼ 24, while the four data points at Ta ¼ 129 imply a lower Re for the bifurcation. Random wavy vortices do not occur without an imposed axial flow. It should be also noted that all the numerical data points are in accordance with the Ta–Re map. Fig. 13 shows velocity vectors and magnitude contours of azimuthal velocity on r–z planes at Ta ¼ 139 and Re ¼ 4:9 for one azimuthal wavelength of the traveling wave. Each frame was taken at a given instant at equal circumferential distance in the azimuthal direction. In Fig. 13(a), the axial velocity profile is included; it is removed in Fig. 13(b). The axial flow is ‘‘winding’’ in Fig. 13(a). Oscillation of the vortices is clearly seen as in Fig. 4(b) after the axial velocity profile is subtracted from the vectors. This spatial observation reveals the oscillation pattern of the wavy vortices better than the temporal approach reported in [17]. They measured time-dependent velocity fields on a given meridional plane; their measurements include not only the oscillation pattern but also the translation movement along the axis due to the imposed axial flow (see Fig. 7

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

109

Fig. 11. Three-dimensional structures of helical vortices; (a) contours of azimuthal component of vorticity (xh ) in a portion of the flow domain, dark color: positive xh , bright color: negative xh , (b) contours of positive J  0:25D2 where J represents Jacobian of shear-stress vectors and D denotes divergence of shear-stresses.

in [17]). For direct comparison of the numerical result with the experimental one, a time-sequence of the same case is shown on a typical r–z plane in Fig. 14. Both the ‘‘winding’’ axial flow (Fig. 14(a)) and the non-uniform translation of the vortex centers (Fig. 14(b)) are clearly seen and well match Fig. 7 in [17]. Especially, the retrograde motion of the vortex centers is also identified as noticed in Fig. 7(b) in [17]. See the retrograde axial motion of the center of vortex 1 in frames 4–5 of Fig. 14(b). In the case with an imposed axial flow, shift-and-reflect symmetry is expressed as follows; vr ðr; h; z; tÞ ¼ vr ðr; h þ Dh; z; tÞ;

ð9Þ

110

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118 25

HV

20

HWV

CP

Re

RWV

15 10 5

LV

WV

0 100

150

Ta

200

Fig. 12. Flow types indicated on the Ta–Re plane, non-vortical Couette–Poiseuille flow (CP, M), non-wavy laminar vortex (LV, ), wavy vortex (WV, j), non-wavy helical vortex (HV, ), helical wavy vortex (HWV, ), random wavy vortex (RWV, }). Notations and solid lines are taken from [17].



vh ðr; h; z; tÞ ¼ vh ðr; h þ Dh; z; tÞ;

ð10Þ

vz ðr; h; z; tÞ  wðrÞ ¼ ðvz ðr; h þ Dh; z; tÞ  wðrÞÞ:

ð11Þ

As shown in Fig. 15, shift-and-reflect symmetry is roughly valid in the case with an imposed axial flow; the symmetry is less perfect than that in the case without an axial flow (c ¼ 0:0294 in Fig. 15). It should be noted that unlike Taylor–Couette flow without an imposed axial flow, this symmetry is applicable only to the remaining flow field after subtracting the annular Poiseuille flow from the computed flow field. In the following two figures, are shown velocity vectors for HWV flow at Ta ¼ 167, Re ¼ 13:1 (Fig. 16), and for RWV flow at Ta ¼ 215, Re ¼ 24:0 (Fig. 17), respectively. Each frame is taken on a typical azimuthal plane during passage of one wavelength of the traveling wave for direct comparison with Wereley and LueptowÕs result [17]. In Fig. 16, the vortex centers persistently move downstream without the retrograde motion (Fig. 16(b)). Fig. 17 shows that the behaviors of the vortices are somewhat ‘‘random’’, and some of the vortices sometimes almost disappear (the 2nd and 6th frames of Fig. 17(b)). This kind of significant variation in vortex strength does not occur in the other types of Taylor vortex. All these features are also experimentally observed. See Fig. 8 for HWV flow and Fig. 9 for RWV flow in [17], respectively. It should be noted that the behaviors of RWV are random in one wavelength of the traveling wave, but quasi-periodic in the azimuthal and the axial directions from the viewpoint of the entire flow domain (Fig. 9(e)). Therefore, velocity vectors for a given r–z plane reveal a random behavior ‘‘quasi-periodically’’ in time. Fig. 18 exhibits the speed of the center of a vortex in Taylor–Couette flow in the axial direction (wvortex ) at various Ta=Re, where wvortex is non-dimensionalized with the mean axial velocity. The theoretical value [21] was obtained for the wavy vortices in the range of 1:6 < Re < 20. The theoretical, experimental, and computational results are all consistent with one another. Furthermore, the tendency persists even beyond Ta=Re ¼ 20.

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

(a)

111

(b)

Fig. 13. Velocity vectors and magnitude contours of azimuthal velocity on r–z planes at Ta ¼ 139 and Re ¼ 4:9 for one azimuthal wavelength of the traveling wave (from top to bottom), contour increment: 0:14Xri . Each frame was taken at a given instant at equal circumferential distance in the azimuthal direction; (a) including the axial velocity profile, (b) with the axial velocity profile removed. The lower and upper lines of each frame represent the inner and the outer cylinder surfaces, respectively.

Phase speeds of the traveling waves of various wavy vortices are presented in Fig. 19; they are normalized by corresponding angular velocity of the inner cylinder. One can notice that the normalized phase speed is almost constant at c=Xi ¼ 0:42 regardless of Ta or Re. In order to quantify the flow instability, we define a non-dimensional variable (Vcl ) on the h–z surface at the center of the annulus (r ¼ ri þ 0:5d); Z H Z 2p 1 jvr ðh; r ; zÞjr dh dz; ð12Þ Vcl ¼ Xi ri 2pr H 0 0

112

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

Fig. 14. Velocity vectors for wavy vortex flow on an r–z plane at Ta ¼ 139, Re ¼ 4:9. Each frame is taken on a typical azimuthal plane during passage of one wavelength of the traveling wave (from top to bottom); (a) including the axial velocity profile, (b) with the axial velocity profile removed. The lower and upper lines of each frame represent the inner and the outer cylinder surfaces, respectively.

which is a quantitative measure of absolute flow rate through the center surface. As flow instability induces vr , the onset of bifurcation can be detected by exponential temporal growth of Vcl . Fig. 20 shows time history of Vcl at various ðTa; ReÞ in a linear-log scale, where time (t) is nondimensionalized with d=Xi ri . In the case of higher Ta for a fixed Re, Vcl grows more rapidly. However, the slopes of the curves for Ta ¼ 123, Re ¼ 4:9 and for Ta ¼ 129, Re ¼ 13:1 are almost identical, which means that a stronger axial flow stabilizes the flow field. Growth rate(r) of Vcl is defined as Vcl  expðrtÞ. Fig. 21 reveals growth rates of Vcl in various cases of Ta and Re, computed from the sloped portion of the curves in Fig. 20; the solid line denotes an approximate boundary between non-wavy and wavy vortices. The smaller the value of Re is, the larger the value of r is at a given Ta. The simulation results corresponding to Ta ¼ 139 tell us that wavy vortices do not show up in the case of stronger axial flow (Re ¼ 24:0). The type of vortices at particular values of Ta and Re can be predicted using Fig. 21 and the corresponding growth rate of Vcl without a ‘‘full’’ simulation. In other words, a big saving in computing time can be achieved because the growth rate of Vcl can be computed in an early stage of full simulation.

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

113

1

ξ

0.8 0.6 0.4 0.2 0

5

(a)

5.25

5.5

5.75

6

6.25

ζ

6.5

6.75

1

ξ

0.8 0.6 0.4 0.2 0

1

1.25

1.5

(b)

1.75

2

2.25

ζ

2.5

2.75

3

Fig. 15. Velocity vectors and magnitude contours of azimuthal velocity on an r–z plane, n ¼ ðr  ri Þ=d, f ¼ z=d, contour increment: 0:13Xri ; (a) reflected velocity vectors of vortex 1 and vortex 2 in the first frame of Fig. 13(b), (b) velocity vectors of vortex 1 and vortex 2 in the fifth frame of Fig. 13(b).

The torque coefficient (CM ) is defined as 2

CM ¼ Ti =½0:5pqðXi ri Þ ri2 H;

ð13Þ

where Ti is the torque required by rotating the inner cylinder at a given speed. Fig. 22 shows CM at various Ta and Re, where the solid line denotes the exact values [22] for non-vortical Couette– Poiseuille flow. One can see that the axial flow significantly reduces CM . In order to elucidate this point more clearly, instantaneous wall shear-stresses (s) along a typical axial line on the inner cylinder surface are presented in Fig. 23, where s is normalized by the wall shear stress (so ) obtained analytically by assuming that no vortices are present. It can be noticed that s=so is larger than unity in most of the range, and this tendency is more evident as Ta increases for a fixed Re. This tells us that the vortices in Taylor–Couette flow actually increase the wall shear stress. This is also consistent with a well-known fact that vortices increase the shear stress. Furthermore, the peaks and valleys in each curve are corresponding to the inner and outer directional ‘‘boundary’’ flows of neighboring vortices, respectively; for example, compare the case of Ta ¼ 139, Re ¼ 4:9 in Fig. 23 with Fig. 15(b). In Fig. 22, it is not surprising that the torque coefficient decreases in the range of Ta > 139 as Ta increases. This is due to the particular way of normalization of Ti ; the denominator of CM contains a square of Xi term. Axial flow, however, has the opposite effect on s=so . Compare the case of Ta ¼ 123, Re ¼ 0 with that of Ta ¼ 123, Re ¼ 4:9 in Fig. 23.

114

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

Fig. 16. Velocity vectors for helical wavy vortex flow on an r–z plane at Ta ¼ 167, Re ¼ 13:1. Each frame is taken on a typical azimuthal plane during passage of one wavelength of the traveling wave (from top to bottom); (a) including the axial velocity profile, (b) with the axial velocity profile removed. The lower and upper lines of each frame represent the inner and the outer cylinder surfaces, respectively.

5. Conclusion In this study, numerical simulations of Taylor–Couette flow with an imposed axial flow were carried out for the identical flow geometry and flow parameters as in the recent PIV measurements of Wereley and Lueptow [17]. The numerical results obtained are in excellent agreement with their experimental results both qualitatively and quantitatively. The vortices in various flow regimes such as non-wavy laminar vortex, wavy vortex, non-wavy helical vortex, helical wavy vortex and random wavy vortex are all consistently reproduced with their experiments. It is confirmed that the shift-and-reflect symmetry is valid for Taylor–Couette flow without an imposed axial flow. The symmetry is roughly valid for Taylor–Couette flow with an imposed axial flow, if the annular Poiseuille flow at the corresponding Reynolds number is subtracted from the computed flow field.

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

115

Fig. 17. Velocity vectors for random wavy vortex flow on an r–z plane at Ta ¼ 215, Re ¼ 24:0. Each frame is taken on a typical azimuthal plane during passage of one wavelength of the traveling wave (from top to bottom); (a) including the axial velocity profile, (b) with the axial velocity profile removed. The lower and upper lines of each frame represent the inner and the outer cylinder surfaces, respectively.

It is revealed that an axial flow stabilizes the flow field as expected. A traveling wave along the azimuthal direction can be detected from instantaneous velocity fields; it is verified that the temporal measurements of Wereley and Lueptow on a given meridional section are valid. In various cases of Taylor and Reynolds numbers, detailed flow information is obtained and the type of the vortices is identified. The computed speed of vortex translation is in consistent with the experimental and theoretical values, and the normalized phase speed of traveling waves is almost constant in the range of Taylor and Reynolds numbers considered in this study. Growth rate of the flow instability is defined and used in predicting the type of the vortices. It is also shown that the torque coefficient on the inner cylinder surface is significantly reduced by the axial flow. The results obtained in this study demonstrate that the numerical simulation employed for the current investigation is suitable and economical.

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118 2

Wereley and Lueptow [17] present theory [21]

1.8 1.6

wvortex /w

116

1.4 1.2 1 0.8 5

10

15

20

25

30

Ta /Re

Fig. 18. Vortex translation speed as a function of Ta=Re. Here, w is the mean axial velocity.

1.2

1

Re=4.9 Re=13.2 Re=24.0 Re=0

c/ Ω i 0.8

0.6

0.4

0.2

0 100

120

140

160

180

200

220

240

Ta Fig. 19. Normalized phase speeds of traveling waves.

10

Vcl

0

10

-1

10

-2

10

-3

10

-4

10

-5

10

-6

10

-7

10

-8

10

-9

Ta=123, Re=4.9 Ta=129, Re=13.1 Ta=139, Re=4.9 Ta=167, Re=13.2 Ta=215, Re=24.0

500

tri Ω i /d

1000

Fig. 20. Time history of Vcl .

1500

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

117

Re=0 Re=4.9 Re=13.1 Re=24.0

0.1

wavy vortex

σ

0.05

non-wavy vortex

120

140

Ta

160

180

Fig. 21. Growth rates of Vcl in various cases of Ta and Re.



Fig. 22. Torque coefficient (CM ) for the inner cylinder vs. Ta at various Re, Re ¼ 0:0 (j), Re ¼ 4:9 (M), Re ¼ 13:1 ( ), Re ¼ 24:0 (}). The linear theory is taken from [22].

Fig. 23. Instantaneous wall shear-stresses along an axial line on the inner cylinder surface.

118

J.-Y. Hwang, K.-S. Yang / Computers & Fluids 33 (2004) 97–118

Acknowledgement This work was supported by grant no. R01-2002-000-00060-0 from the Basic Research Program of the Korea Science & Engineering Foundation.

References [1] Taylor GI. Stability of a viscous liquid contained between two rotating cylinders. Phil Trans R Soc 1923;A223:289– 343. [2] Stuart JT. On the non-linear mechanics of hydrodynamic stability. J Fluid Mech 1958;4:1–21. [3] Davey A. The growth of Taylor vortices in flow between rotating cylinders. J Fluid Mech 1962;14:336–68. [4] DiPrima RC. Vector eigenfunction expansions for the growth of Taylor vortices in the flow between rotating cylinders. In: Ames WF, editor. Nonlinear Partial Differential Equations. Academic; 1967. p. 19–42. [5] Davey A, DiPrima RC, Stuart JT. On the instability of Taylor vortices. J Fluid Mech 1968;31:17–52. [6] Eagles PM. On the stability of Taylor vortices by fifth-order amplitude expansions. J Fluid Mech 1971;49:529–50. [7] Brandstater A, Swinney HL. Strange attractors in weakly turbulent Couette–Taylor flow. Phys Rev A 1987; 35:2207–20. [8] Wereley ST, Lueptow RM. Azimuthal velocity in supercritical circular Couette flow. Exp Fluids 1994;18:1–9. [9] Jones CA. The transition to wavy Taylor vortices. J Fluid Mech 1985;157:135–62. [10] Schr€ oder W, Keller HB. Wavy Taylor-vortex flows via multigrid-continuation methods. J Comput Phys 1990;91:197–227. [11] Marcus PS. Simulation of Taylor–Couette flow, Part 2, Numerical results for wavy-vortex flow with one traveling wave. J Fluid Mech 1984;146:65–113. [12] Chung KC, Astill KN. Hydrodynamic instability of viscous flow between rotating coaxial cylinders with fully developed axial flow. J Fluid Mech 1977;81:641–55. [13] Gravas N, Martin BW. Instability of viscous axial flow in annuli having a rotating inner cylinder. J Fluid Mech 1978;86:385–94. [14] Takeuchi DI, Jankowski DF. A numerical and experimental investigation of the stability of spiral Poiseuille flow. J Fluid Mech 1981;102:101–26. [15] Meseguer A, Marques F. On the competition between centrifugal and shear instability in spiral Poiseuille flow. J Fluid Mech 2002;455:129–48. [16] Wereley ST, Lueptow RM. Spatio-temporal character of non-wavy and wavy Taylor–Couette flow. J Fluid Mech 1998;364:59–80. [17] Wereley ST, Lueptow RM. Velocity field for Taylor–Couette flow with an axial flow. Phys Fluids 1999;11(12):3637–49. [18] Rosenfeld M, Kwak D, Vinokur M. A fractional step solution method for the unsteady incompressible Navier– Stokes equations in generalized coordinate systems. J Comput Phys 1994;94:102–37. [19] Lueptow RM. Personal communication, 2001. [20] Hunt JCR, Abell CJ, Peterka JA, Woo H. Kinematical studies of the flows around free or surface-mounted obstacles; applying topology to flow visualization. J Fluid Mech 1978;86:179–200. [21] Recktenwald A, Lucke M, Muller HW. Taylor vortex formation in axial through-flow: Linear and weakly nonlinear analysis. Phys Rev E 1993;48:4444–54. [22] Schlichting H. Boundary layer theory. McGraw-Hill; 1979.