On the flow between two rotating shrouded discs

Computers&FluidsVol. 14, No. 3, pp. 183-196, 1986

0045-7930/86 $3.00+ 0.00 PergamonJournalsLtd

Printed in Great Britain

ON

THE

FLOW BETWEEN TWO SHROUDED DISCS

ROTATING

P. W. DUCK Department of Mathematics, University o f Manchester, Great Britain

(Received 1 August 1984; in revised form 19 March 1985) Abstract--Numerical solutions to the full Navier Stokes equations for the flow enclosed between two rotating discs and a sidewall are presented. The sidewall remains fixed throughout, while we consider examples in which: (i) the bottom disc is rotating and the top disc is fixed; (ii) the two discs rotate in the same sense; (iii) the two discs rotate in opposite senses. We present solutions up to a Reynolds number o f 1200, using second order central finite differencing. 1. I N T R O D U C T I O N

The flows involving rotating discs have received considerable attention. The most fundamental of these, namely that of the flow above a single, infinite rotating disc, was treated originally by von Karman [1], who showed how the Navier-Stokes equations could be reduced, in tiffs case, to a set of ordinary differential equations, using a similarity transform. Batcbelor [2] extended tiffs work to the flow between two rotating discs and argued that in the limit of large Reynolds number, boundary layers develop on each disc, while the core of the flow will rotate with constant angular velocity. Stewartson [3] studied the same problem, and argued that in this limit, the main body of the fluid, outside of the boundary layers, did not rotate, in direct contrast with Batchelor's [2] suggestion. More recent works on this problem (both experimental and numerical) have since shown that both suggestions are applicable, the point in parameter space determining which particular model is taken. A number of authors have found multiple solutions of this problem, for example Mellor, Chapple and Stokes [4], Roberts and Shipman [5], and Nguyen, Ribault and Florent [6]. In fact, multiple solutions have been found, more recently, to occur even for the original yon Karman [1] problem, by Zandbergen [7] and Dijkstra [8]. From the practical point of view, however, the most important class of problem involves confined disc flows. These are particularly relevant in chemical engineering, for example in the use of gas centrifuges. However, problems of this class are no longer amenable to similarity transformations, but instead involve solutions of partial differential equations (i.e., the Navier-Stokes equations) and hence require fully numerical solutions. The problem of one rotating disc and one fixed disc, enclosed by a fixed shroud or sidewall, has been studied by Pao [9] and by Tomlan and Hudson [10], while Lehmkuhl and Hudson [ 11] have studied this problem experimentally. A number of authors have also studied the case of one fixed disc, one rotating disc with the sidewall also rotating (with the same rate as the rotating disc), including Pao [12] and Lugt and Haussling [13]. These latter authors also considered the case of counter-rotating discs (the sidewall rotating with the same rate as one of the discs). All these studies considered configurations of aspect ratio unity. Recently Dijkstra and van Heijst [ 14] have considered this last problem for small aspect ratio, and verified their numerical results experimentally. In this paper we consider the related problem of a rotating upper and lower disc, but with a fixed sidewall (we shall also include two computations involving a fixed shroud, one fixed disc, and one rotating disc in order to verify our numerical scheme). We study a number of cases including counter-rotation of the discs, and also examples in which the discs rotate in the same sense. 2. E Q U A T I O N S

OF

MOTION

AND

NUMERICAL

SCHEME

Consider a closed circularcylinder of height H and radius R, which is completely flied with incompressible fluid of density p*. The bottom surface rotates about the axis with 183

184

P . W . DUCK

rotational rate ftt s -t, and the upper surface rotates about the axis with rotational rate ft2 s -t. We assume the flow to be steady, laminar, and axisymmetric, and take r = r * / R and z = z * / R to be the dimensionless radial and vertical coordinates, respectively. The equations of motion determining the flow may then be written (see, for example, Tomlinson and Hudson [ 10]) Ou v2 0 u Op 1 {O[ld 1 02ul u -- - -- + w . . . . + (ru) + Or r Oz Or Ree g ~r Oz2J ' ov uv ov 1 fat'lO q 0:v~ u ~rr + ~r + w --0z = --Rel o r L r ~ r I ( ~r +v ~) z 2 ,

Ow Ow Op+ 1 {10[rOW ] 02w t u -~r + w O--'z = - O--z Ree r Or k Or] + Oz2J '

(2.1)

(2.2) (2.3)

together with 1 O Ow r Or (ru) + -~z = 0.

(2.4)

Here u = u * / R f t t , v = v * / R f t t , w = w * / R f t l denote radial, azimuthal and vertical dimensionless velocities, respectively, while p = p * / p * R 2 ft21 denotes dimensionless pressure. Then the (rotational) Reynolds number is defined by Re = R 2 f t t / v .

(2.5)

The boundary conditions to be applied, generally, are those corresponding to no slip and symmetry, namely Ow

u=v=--=O Or

on

r=0,

u=v=w=O

on

r=l,

u = w = O,

v = ar

on

z=A,

u=w=O,

v=r

on

z=O.

(2.6)

Here we have written A = H / R (the aspect ratio) and a = ft2/ftt, the ratio of the two angular velocities. However, from a computational standpoint it is found convenient to eliminate the pressure from (2.1)-(2.4) to yield (see, for example, Goldstein [ 15]) 2fl Oft r 2 0z +

10(~b, ~') -

-

r O(r, z)

2~" O~b +

r 2 Oz

1 Re

= --

D2~'

1 0(~b, ft) = 1 D2ft ' r 0(r, z) Re g"= D2~b.

(2.7) (2.8) (2.9)

~k is the streamfunction defined by

l O¢ U "= - - - - ~ r OZ '

1o¢ w = - --

(2.10a)

r :Jr '

ft is defined by v = ft/r.

(2.10b)

Note that fl is not the angular velocity, g"is related to the vorticity w by = -wr.

(2.10c)

How between two rotating shrouded discs

185

Further note that D2mr0

(1 0 )

d2

(2.11)

Or

The appropriate boundary conditions are then

o6

on

z = 0

and

z = A,

6=~r=0

on

r = 1

and

r = 0,

[~ = r 2

on

z = 0,

[2 = tzr 2

on

z = A,

fl=0

on

r= 1

and

r -- 0.

6=~z=0

o6

(2.12)

Equations (2.7)-(2.9) were then approximated using second order, central, finite differeneing. If Az and Ar denote vertical and radial step sizes, respectively, then we made the approximation ~zz(r, z)

6(r, z + Az) - 6(r, z - Az) 2Az + O((Az) 2)

(2.13)

(and similarly for o~lOr, and analogous derivatives of fl and D and 026

6(r, z + Az) -- 26(r, z) + 6(r, z -- Az)

Oz2(r, z)

( az) 2

+ O((Az)2).

(2.14)

The manner in which the conditions at the edge of the domain are imposed is not unique. One method involves utilising the zero normal derivative condition directly (by, for example, employing one sided differencing). However, we preferred to use an approach analogous to that used by (for example) Burggraf [ 16], namely by specifying 6 and fl along the edges, together with ~'(1, z) = 26(1 - Ar, z)/(Ar) 2, ~'(r, 0) = 26(r, Az)/(Az) 2, ~'(r, A) = 26(r, a - Az)/(Az) 2, ~'(0, z) = 0.

(2.15)

The first three of these expressions were obtained by combining (2.9) with the zero normal derivative requirement, while the final condition is determined by demanding symmetry/antisymmetry about r = 0, together with boundedness along r = 0. Notice that strictly fl(0, 1) and fl( l, l) are in general undefined due to the discontinuity in angular velocity between the angular velocities of the discs, and the fixed shroud (however, see next section). Fortunately (from the numerical point of view) the value of fl at these points is not required in our scheme. We chose to solve the system along lines z = constant, and the result is a system of blocks centered on the diagonal (with each block comprising 9 X 3 elements), and as such is amenable to standard Gaussian elimination routines. Each sweep of the domain extended from r = Ar, up to r = 1 - Ar, and was repeated until convergence was achieved (typically maximum computed dependent variable change < 10-6). Note that all the nonlinear terms were included on the right-hand side of the difference equations (this amounts to a form of Picard's method). In some cases, at larger Reynolds numbers, a certain amount of underrelaxation was required (details of which are given in Section 3).

186

P.W. DUCK

3. RESULTS 3,1 Rotating lower disc, fixed shroud and upper disc (a = O) The case of one rotating disc, with the shroud and other disc fixed (corresponding to a = 0) has been considered by Tomlan and Hudson [10] up to Re = 100, by Pao [9] up to Re = 400 and by Lugt and Haussling [13] up to a Reynolds number of 1000. In all cases the aspect ratio A was set equal to unity. Following these previous works we also restricted attention to the A = 1 configuration, and these previous results provided a useful check on our numerical results. All the results shown here were obtained using a uniform mesh of Ar = Az = 0.0125, and repeated on a control grid of Ar = Az = 0.025 to assess accuracy. Figure l (a) shows the flow pattern in the axial plane for Re = 400, and may be compared to Fig. 7 of [9] and Fig. 13 of [ 13], The agreement is satisfactory, although there does exist one discrepancy with [9], namely a very small separating streamline in the top right-hand corner of the present results, which is not shown in the results of [9]. A similar feature was seen to occur in [13], and consistently throughout our computations, and is also visible in the Re = 100 results of Tomlan and Hudson [10]. Further, this small recirculating re#on was also present (and of similar dimension) in the control computation on the coarse grid. Figure l(b) shows lines of constant angular velocity (v/r) at Re = 400, in the radial plane, while Fig. l(c) shows the radial distribution of angular velocity at a number of z stations. The most noticeable feature of these results [confirmed by Fig. l(b)] is the increase in angular velocity at the join of the fixed shroud and the rotating disc (i.e., r --- 1, z "~ 0). Note that strictly at z = 0, r = 1 the angular velocity is undetermined due to the discontinuity in boundary conditions (as noted previously). In fact, it is a simple matter to study the local behaviour of f~ (and hence of the angular velocity) in this region. From eqn (2.8), we may expect that the motion will be dominated by the viscous terms in this corner region, and hence the governing equation for f~ in this region is (approximately)

0 .8

6

4

.2

I .2

I .4

I .6

I .8

0

t.

2

I 4

6

.B

.6

8

t

r

(b)

(0)

vlr I.

.8

.6 Zmj 4

.2

I 2

I 4

I 6 (C)

I .9

I 1

f

0

2

4

I.

r

(d|

Fig. 1. a = 0, Re = 400. (a) Streamlines. (b) Lines of constant angularvelocity. (c) Angularvelocity distributions. (d) Lines of constant vorticity.

Flow betweentwo rotatingshroudeddiscs 02fl 1 0fl 1 02fl ---T + + 7 a , = o,

187 (3.1)

where P = {z 2 + ( r -

1)2} 1/2,

( F a 1),

~p = t a n - l { z / ( r -

1)},

~/2 < $ < ~ .

(3.2)

The boundary conditions to be applied are fl = 0

or

~ = 7r/2,

fl= 1

or

$=lr.

(3.3)

n = (2/~)(@ - ~/2).

(3.4)

The solution of this is quite simply

Figure l(d) shows lines of constant vorticity (w) for Re = 400, and indicates clearly the presence of some form of boundary layer on the lower (rotating) disc, and also on the shroud. Finally, it is worth mentioning that this particular computational task took approximately 1 hour of computational time on a CDC 7600. The next set of results shown is for Re = 700. Above this value of Reynolds number our scheme very rapidly lost numerical stability, and for this case an underrelaxation parameter of 0.1 was employed, although no optimisation of this value was carried out. (Lugt and Haussling [ 13] present "almost steady" solutions for the higher Reynolds number of 1000, obtained using a time marching scheme to seek the steady state solution, a technique which is well known to be more stable, but also much more time consuming; these calculations [13] were incidentally carried out on a coarser, albeit stretched, grid.) Figure 2(a) shows the flow pattern in the axial plane. Again the small separated flow region described previously is present in the top right-hand corner of the axial plane. There is in fact a very slight reduction in size, compared with the previous (Re = 400) example-see Fig. 1(a). Generally we see the axial flow is weaker than the Re = 400 case, as least in the core of the flow. Comparing the angular velocities (v/r) for the two Reynolds numbers [Figs. 2(b) and l(b)], the larger Reynolds number example exhibits generally smaller angular velocities in the bulk of the core, however close to the axis of rotaiton, in the upper half of the plane, the angular velocity is seen to increase with an increase in Reynolds number. This feature is confirmed by the angular velocity distributions at selected z stations shown in Fig. 2(c). These distributions again exhibit a peak as z ---. 0 and as r --, 1, these peaks being slightly closer to the fixed shroud than in the corresponding Re = 400 example. Figure 2(d) shows the distributions of angular velocity at selected r stations. Close to the moving disc the angular velocity of the fluid is fairly independent of r. However it is by no means conclusive that this an example of the Batchelor [2] model, which predicts that the bulk should rotate at an angular velocity of 0.3. Indeed, combined, Figs. l(c), 2(c), 2(d) show no firm signs of approaching this limit for r small. The vorticity lines for Re = 700 are shown in Fig. 2(e), and again point to a boundary layer on the rotating disc and also the shroud. In addition, there are signs of a boundary layer developing on the upper (fixed) disc. Notice, also, that these observed boundary layers are thinner than in the previous Re = 400 calculations. As an indication of the effect of grid size on solutions we present, in Fig. 2(f), lines of constant angular velocity for the Re = 700 case, obtained using our control grid. The agreement with our finer grid results [Fig. 2(b)] is seen to be satisfactory, with the largest discrepancy being seen in the region of the join of the rotating surface and the shroud (as might well be expected). One final point is that the most effective means of generating results at increasingly large Reynolds numbers was by gradually incrementing the Reynolds number. This was the means by which we generated the Re = 400 and 700 results. We also adopted this technique for the a ~ 0 cases to be described next.

188

P. W. DUCK l

-.0 I O

I

2

4

I

1

6

8

I 2

r

(0)

I 4

I .6

I 8

t

(b) v/r

v/r

.8

.8

6

6

~

z ".7

r=O /

r. 1 r=2

r=.5

r=5

r=o,.

r = 2 r-.7

~ "2

2

0

2

2

.4

6

8

t.

r

r-.7 r'.9

I

I

I

I

I

2

4

6

.B

1

I .8

1

Z

(d)

(C) z

z

I

25

18

~

6 1

.2

t .~

' 2

I ,

' 4

I .6

, 8

1

0

I 2

I 4

I .6

iw, r

(e)

Fig. 2. a = 0, Re = 700. (a) Streamlines. (b) Lines of constant angular vorticity. (c) Angular velocity distributions (z = constant). (d) Angular velocity distributions (r = constant). (e) Lines of constant vorticity. (f) Lines of constant angular velocity (Ar = ~ z = 0.025).

3.2 Rotating upper and lower discs, fixed shroud (o~ ~ O) We next turn our attention to a second class o f problem, namely that o f a fixed shroud, with both a rotating base and top. This particular configuration does not appear to have been treated in the past, previous studies have concerned themselves with the shroud rotating with the same angular velocity as one o f the discs. Again we confine our calculations to the case A = 1, and throughout we use our uniform (fine) grid z~r -- Az = 0.0125. The first set of results so obtained, are shown in Figs. 3(a), (b), and are for the case Re = 0, a = 1/2, i.e. the upper disc is rotating in the same direction as the lower disc, with half the angular velocity o f the lower disc. The lines of constant angular velocity are shown in Fig. 3(a), while radial angular velocity distributions, at selected z stations, are shown in Fig. 3(b). Note that for Re = 0, ~"-~ ff = 0, for all c~. Indeed, analytic solutions for the

Flow between two rotating shrouded discs

189

v/r

z l.¶ ,8

.4 2 I

I

|

I

~2

"4

6

.8

m ~l

r

I

I

I

I

!

2

.4

8

18

t.

0

(a)

r

(b)

Fig.

3. a = 1/2, R e = 0. (a)

Lines of constant angular velocity. (b) Angular velocity distributions.

angular velocity are possible when R = 0 (see, for example, Pao [ 12]). However, the analysis is straightforward (and is not repeated here), and results in Bessel series. In fact, it would appear that from the computational point of view, our finite difference approach is probably the more straightforward to obtain numerical results. The next set of results, shown in Figs. 4(a)-(d), is for the example a = 1/2 (again), with Re = 200. Perhaps the most important feature of these results are the two stagnation points which are seen to occur [Fig. 4(a)] on the upper disc, at around r = 0.23, and on the shroud, near z = 0.63. Associated with these stagnation points is the (quite) large counterrotating flow region in the upper half of the domain. These results also suggest that a boundary layer collision is occurring on the shroud, near the stagnation point, with the subsequent development of a shear layer. Z

Z

I

1. .4 .3

lB

,8

16

.6

14

.4

.2

2

0

I

I

I

I

2

4

.6

8

t

r

0

.

.1

I

I

I

I

2

4

6

.8

(a)

1.

r

(b)

v/t

z 1.

.8 Z

1 l~

Z.2

Z" 8 ~

O

I

I

I

I

2

14

16

,8

(c)

1.

.4 z''5 .2

r

O

I

I

I

I

2

4

.6

8

I.

(d)

Fig. 4. a = 1/2, Re = 200. (a) Streamlines. (b) Lines of constant angular velocity. (c) Angular velocity distributions. (d) Lines of constant vorticity.

190

P.W.

DUCK

The corresponding lines of constant angular velocity are shown in Fig. 4(b), while the radial angular velocity distributions are shown in Fig. 4(c). In the vicinity of the lower disc, these results are much changed from the Re = 0 results [Figs. 3(a), (b)]. The increase of angular velocity close to the join of the fixed shroud and the rotating disc (also observed in the a = 0 examples) appears to be a nonlinear effect [not present in Fig. 3(b)] and this time also occurs near the top (slower rotating) disc. A further interesting feature of the Re = 200 results is the apparent "bunching" of these distributions near the axis of symmetry, indicating an approximately uniform angular velocity along the axis of symmetry. Finally, lines of constant vorticity are shown in Fig. 4(d). The next choice of parameters is ~ = 1/2, Re = 400. From Fig. 5(a) we see that at this higher Reynolds number the stagnation point on the upper disc has moved radially outwards, to around r = 0.37, while the stagnation point on the shroud is little changed in position, at around z = 0.64. Lines of constant angular velocity are shown in Fig. 5(b), and the corresponding radial angular velocity distributions are shown in Fig. 5(c). These distributions again indicate a local peaking of the angular velocity close to r = 1 and z = 0 on z = 1. Although these peaks are approximately the same size as those for Re = 200, they are rather more pronounced, indicating a thinning of the boundary layer. This trend is confirmed by comparing the vorticity lines in Fig. 4(d) with those in Fig. 5(d). We now turn to two examples for which cz = 1, corresponding to the case of the upper and lower discs possessing identical angular velocities. For this choice of c=, regardless of the Reynolds number, the problem possesses a natural symmetry/antisymmetry about the z = 1/2 plane. In particular we have if(r, z -

1/2) = -~b(r, - z

+ 1/2),

co(r, z -

1/2) = - c o ( r , - z

+ 1/2),

(3.5)

z

t, 0

0 .8

6

4

.2

o

.z

.4

.6

.e

~,

r

o

|

I

l

I

.z

.4

.s

.e

(o)

~.

{b)

v/r

I, ---I. .8 z = .t

/

Z = .9

.6

Z" 5

Z=

Z'A

.4

Z=.9

.2

I

!

I

I

I

.2

.4

6

8

{.

(c]

r

0

I

I

I

I

.Z

4

.6

,8

I.

(d}

Fig. 5 . , = 1/2, Re = 400. (a) Streamfines. (b) Lines of constant angular velocity. (c) Angular velocity distributions. (d) Lines of constant vorticity.

Flow between two rotating shrouded discs

191

v/r

I,

Z = .05 Z=.l Z= 2 1=.3

6

,6

.5 .2

I

I

I

I

.2

.4

.6

8

"-

t.

r

0

I

I

I

I

2

4

6

e

(a)

Fig.

"~

r

(b)

6. c~ = 1, R e = 0 . (a)

Lines of constant angular velocity. (b) Angular velocity distributions.

While fl(r, z - 1/2)

=

f l ( r , - z + 1/2)

(3.6)

(implying that the angular velocity is symmetrical about z = 1/2). Consequently, in the figures relating to a = 1 (Figs. 6 and 7), we shall show only half the range o f z . Figure 6(a) shows lines of constant angular velocity for Re = 0, while Fig. 6(b) shows the corresponding radial angular velocity distributions at selected z stations. As in the previous Re = 0 example (for a = 1 / 2 ) , we see these distributions decline monotonically with r. The next computation was performed with a = 1, Re = 400. Figure 7(a) shows the streamlines for this flow, and indicates (as may be anticipated from symmetry considerations) a shear layer lying along z = 1/2. In addition, we see a (small) recirculatory region close to the point of symmetry (i.e., r = 0, z = 1/2), giving four counterrotating regions in all. Lines of constant angular velocity are shown in Fig. 7(b), and the corresponding distributions are in Fig. 7(c). The previously mentioned "peaking" of these distributions again Z

.6

.4

.8

175

.5

.2

,2

0

I 2

I .4

I .6

I ,2

I .8

1.

I .4

r

I .6

I .8

1.

I .6

I ,8

I.

(b)

(a) vlr

.8

z

= .6

.6

.4

.4

.2

.2

.2

.4

.6

8

4.

(c)

Fig.

7.,

r

0

I .2

I .4

r

(d) ( b ) Lines of constant angular velocity. (c) Angular velocity distributions. (d) Lines ofconstam vorticity.

= 1, R e = 4 0 0 . (a) S t r e a m l i n e s .

192

P.W. DUCK

occurs, while again the angular velocity along the axis of rotation is approximately constant. Vorticity contours are shown in Fig. 7(d). This figure suggests that a boundary layer collision occurs on the shroud, close to z = 1/2, although the resulting shear layer appears rather weak, at least for this moderate value of Re. The next choice of a was - 1 / 2 , i.e. the upper disc is rotating in the opposite sense, and with half the rate of that of the lower disc. The first Reynolds n u m b e r taken was Re = 0, and in Fig. 8(a) we see that to quite a good approximation the dividing line between the clockwise and anticlockwise flow (in the azimuthal plane) lies along z - 0.6. The corresponding distributions [Fig. 8(b)] again exhibit a monotonic decrease/increase towards zero. Figures 9 relate to the example a = - 1 / 2 , Re = 400 and the first of these [Fig. 9(a)] shows two regions of counterrotating fluid, with the dividing streamline lying entirely between z - 0.66 and z - 0.73. The line of zero angular velocity [see Fig. 9(b)] has become more deformed compared to the Re = 0 case, and to a first approximation follows the path of the dividing streamline. Once again the angular velocity distributions [Fig. 9(c)] generally exhibit peaks close to r = 1, although, perhaps surprisingly the distribution for z = 0.95 shows only a very weak m a x i m u m . This is presumably because this line z = 0.95, lies well inside the upper boundary layer (confirmed by the vorticity contours Fig. 9(d) in which the angular velocity is fairly independent of r. The final (and extreme) value of a taken was a = - 1 , with the angular velocity of the top disc being equal and opposite to that of the lower disc. Because of the antisymmetry of the problem about z = 1/2, we chose to halve the computational domain by setting ~(r, z = 1/2) = g'(r, z = 1/2) = 9(r, z = 1/2) = 0,

for all r,

(3.7)

and solved in just 0 ~< z ~< 1/2 (but still used Ar = AZ = 0.0125). Consequently in presenting results for this final choice of a, we shall show results for 0 < z ~< 0.5 only. Figure 10(a) shows the lines of constant angular velocity for the case Re = 0, and Fig. 10(b) the corresponding distribution, which shows similar features to the previously discussed Re = 0 cases. We now turn to the example of Re -- 400, a = - 1, results for which are shown in Figs. l l(a)-(d). The flow pattern is presented in Fig. 1 l(a) (remembering that there will exist an equal, contrarotating eddy in the upper half plane z > 1/2). This feature also suggests the existence of a shear layer along z = 1/2, similar to that found in the a = 1 case, caused by some form of boundary layer collision near r = 1, z = 1/2. Figure l l(b) shows angular velocity contours, and Fig. l 1(c) the corresponding radial distributions. The vorticity contours [Fig. 1 l(d)] confirm the presence/development of the shear layer along z = 1/2. v/r

I,

.8

.6

z" .t

Z" .25

.4

-.2S

.2

0 .I

0

I .2 Z ".7S

I .4 y

I

I .

-2 Z s.9 -.4 I .2

I .4

I 6

(a)

I .8

I. (b]

Fig. 8. a = -1/2, Re = 0. (a) Lines of constant angular velocity. (b) Angular velocitydistributions.

193

Flow between two rotating shrouded discs

1.

,g

.6

.4

.2

I

I

I

.2

.6

.8

0

~.

I

I

I

I

2

.4

6

.8

1

.8

1.

(b)

(a) v/r

t.

.8

.6

Z z • .0

-

.4

z • .1

-.

-

-.

t

.2

.8

0

'

.6

-.2

4 z ~ .95

-.4

.2

-.6

o

2

.4

(c)

6

r

(d)

Fig. 9. a = - 1 / 2 , Re = 400. (a) Streamlines. (b) Lines of constant angular velocity. (c) Angular velocity distributions. (d) Lines of constant vorticity.

Figures 12(a)-(d) relate to the example Re = 800. It is interesting to note (and also perhaps a little surprising) that this contra-rotating configuration was a good deal more stable, from the numerical point-of-view, than the other cases tackled, in particular that of the stationary upper disc case (a = 0) considered previously. v/r

1

Z'.05 8

Z " .1

.6

Z'2

.4

.2

0

I

I

I

I

I

2

4

6

8

t

(a)

r

~

Z'.3 Z'4

~

I

I

I

2

4

.6

.8

I.

r

(b)

Fig. 10. a = - 1 , Re = 0. (a) Lines of constant angular velocity. (b) Angular velocity distributions.

P. W. DUCK

194

I

I

I

I

.2

4

.6

B

0

I.

I

I

I

I

2

4

.6

8

(a)

t

(b)

v/r

I

8 z 6

6 Z '--05

4

2

2

I

I

I

I

2

4

6

8

I 2

r

i

I

#,

.6

(c)

Fig.

I 8

I

(d) ( b ) Lines of constant angular velocity. (c) Angular velocity distributions. (d) Lines of constant vorticity.

11. c~ = - 1, R e = 4 0 0 . (a) S t r e a m l i n e s .

Z

.6 .t

0

.4 .3

.2

.2

I

I

I

I

2

.4

6

.8

F 1

0

r

I

I

I

I

.2

.4

6

.8

(a)

t.

r

(b)

v/r

1.

.8

z Z" ,05

.6

.6

4

2

.2

I

I

I

I

I

12

14

.6

I8

~I

(c)

0

.2

.4

.6

8

r

(d)

Fig. 12. a = - 1, Re = 800. (a) Streamlines. (b) Lines of constant angular velocity. (c) Angular velocity distributions. (d) Lines of constant vorticity.

Flow between two rotating shrouded discs

195

z

.6

,4

.2

0

I

#

I

I

.2

.4

.6

.8

1

r

0

I

I

I

I

,2

4

.6

8

1.

6

.8

4.

(o)

r

(b)

v/r 1

z

Z • .05

.6

~°.t .2

0

.2

4

.6

.0

1.

r

0

2

(c)

4

r

(d)

Fig. 13. a = - 1 , Re = 1200. (a) Streamlines. (b) Lines of constant angular velocity. (c) Angular velocity distributions. (d) Lines of constant vorticity.

Figure 12(a) shows a general weakening of the radial and axial flow in the core, when compared to Fig. 11(a) while the flow near the shroud is rather more intense than in the previous, lower Reynolds number. Comparing the angular velocities however [Fig. 12(b)], with the corresponding contours for Re = 400 [Fig. 11(b)], there is a general increase in the angular velocity in the core, particularly close to the axis of the cylinder. This trend is confirmed in Fig. 12(c) showing the angular velocity distributions at selected z locations. Further these distributions exhibit increasingly sharp maxima, close to the shroud. Vorticity contours are shown in Fig. 12(d). Figures 13(a)-(d) are for the example Re = 1200, the highest Reynolds number attempted (the numerical scheme experienced a deterioration in stability as Re was increased beyond this value); interestingly, no underrelaxation was required at this value of Reynolds number. The general trends outlined above were repeated. One further trend also becomes clear, namely the rapid decay of the angular velocity away from the moving disc. For example on the axis of the cylinder, at just a distance of 0.05 above the disc, the angular velocity is just 23% of the angular velocity on the disc itself. Close to the axis of rotation, however, if we then move further from the lower disc, the fluid experiences an increase in angular velocity, before reducing again as the stationary disc is approached. At the largest value of Reynolds number taken, with the grid size used (0.0125), it is probable that some deterioration in accuracy is beginning; as a rough estimate, we may expect the thinnest boundary layer on the walls to be O(Re-'/2), i.e. 0(0.03), and so certainly for larger Reynolds numbers a grid refinement would be necessary. This, combined with the deterioration in stability appears to be a restriction on the iterative and difference schemes used. In the next section we go on to draw several general conclusions arising from this study. 4. C O N C L U D I N G

REMARKS

One point worthy of mention from the onset is that none of our calculations showed any signs of nonuniqueness, for any value of a. In the case of unconfined rotational flows,

196

P.W. DuCK

nonuniqueness appears to be a c o m m o n p h e n o m e n o n , while no nonuniqueness has been reported to date in the case o f confined rotating flows o f this class. It is quite clear that the shroud enclosing the discs has a quite p r o f o u n d effect on the flow. For example, even close to the axis o f rotation, the walls are seen to exert a strong influence; this is evidenced by the various spatial undulations o f the streamlines close to r = 0 [see, for example, Fig. 3(a)]. A second feature illustrating this point is the increase in angular velocity in certain regions close to r = 0, for example, Fig. 13(b). Further (and as a consequence), our numerical results suggest neither the Stewartson [3] or Batchelor [2] similarity solutions are applicable, even close to the axis o f s y m m e t r y , at least for the m o d e r a t e range of Reynolds n u m b e r s tackled in this study. Finally we add a note of caution. A n u m b e r o f the flows calculated possess shear layers [for example, see Figs. 1 l(a), 12(a), 13(a)]. Similarity flows of this general type are discussed by Batchelor [2], who also noted that flows of this type are likely to be very unstable (at least in the limit of large Reynolds numbers), a n d as such are likely to be difficult to realize experimentally. Related to the presence of these shear layers it would appear that unless the position of these is k n o w n a priori and they lie along lines o f constant coordinate, then it would appear that the use of n o n u n i f o r m grids m a y be rather inappropriate, and a uniform grid (such as used in this study) is s o m e w h a t m o r e justifiable. Acknowledgement--The author wishes to thank Dr. S. D. R. Wilson for arousing his interest in problem of

this type.

1. 2. 3. 4. 5. 6. 7. 8.

9. 10. 1I. 12. 13. 14. 15. 16.

REFERENCES T. yon Karman, Z. Angew Math. Mech. 1, 233 (1921). G. K. Batcheior, (2. J. Mech. Appl. Math. 4, 29 (1951). K. Stewartson, Proe. Camb. Phil. Soc. 49, 333 (1953). G. L. Mellor, P. J. Chapple and V. K. Stokes, Z FluidMech. 31, 95 (1968). S. M. Roberts and J. S. Shipman, J. FluidMech. 73, 53 (1976). N. D. Nguyen, J. P. Ribault and P. Florent, J. FluidMeeh. 68, 369 (1975). P.J. Zandbergen, Springer Lecture Notes in Mathematics Vol. 771, p. 563. Springer-Verlag, New York (1979). D. Dijkstra, J. Engng Math. 14, 133 (1980). H. P. Pao, Phys. Fluids 19, 4 (1971). P. F. Tomlan and J. L. Hudson, Chem. EngngSci. 26, 1591 (1971). G. D. Lehmkuhl and J. L. Hudson, Chem. Engng $¢i. 26, 1601 (1971). H. P. Pao, J. Appl. Mech. 37, 480 (1970). N.J. Lugt and H. J. Haussling, Acta. Mech. 10, 255 (1973). D. Dijkstra and G. J. F. van Heijst, J. FluidMech. 12~, 123 (1983). S. Goldstein (Editor). Modern Developments in Fluid Dynamics. Clarendon Press, Oxford (1938). O. R. Bur~raf, J. FluidMech. 24, 113 (1966).