- Email: [email protected]
A finite difference simulation of a continuous flow centrifuge
Computer Methods in Applied Mechanics and Engineering 93 (1991) 401-414 North-Holland
A finite difference simulation of a continuous flow centrifuge J.W. Frederick ~, R.J. Ribando and H.G. Wood Department of Mechanical and Aerospace Engineering, School of Engineering and Applied Science, University of Virginia, Charlottesville, VA 22901, USA Received 8 March 1990
The secondary flow of an incompressible liquid in a continuous process, vertical tubular centrifuge is studied using finite difference solutions to the governing equations. In particular the secondary flow associated with the addition and removal of mass or momentum from the container is computed and analyzed. For parameters typical of the liquid centrifuges studied here, the governing equations are linear. Thus various feed and withdrawl configurations may be analyzed using superposition of a few solutions to the governing equations. Even aside from the secondary flow induced by the density difference between the species being separated, the addition and withdrawl of fluid is found to have a major impact on the secondary flow pattern.
1. Introduction
In this paper a hypothetical tubular centrifuge is analyzed using finite difference solutions to the governing equations for a homogeneous fluid. While a real centrifuge obviously entails a density difference either between two fluids or between a fluid and solid particles dispersed in it, this study is restricted to species of negligible property difference so that the dynamical effects of perturbations due to sources and sinks of mass and momentum can be studied. For most cases of interest it can be shown that the governing equations are linear to good approximation (as measured by the Rossby number, a non-dimensional parameter to be discussed later), and thus superposition of separate simple effects can be used as a tool in calculating and analyzing more complicated flows. This approach is used extensively in the design and analysis of gas centrifuges for uranium enrichment [1]. In gas centrifuge analysis the effects of withdrawl scoops and mass injection have been simulated using point or distributed sources of mass and momentum [2]. There have been very few studies other than in the gas centrifuge literature involving strongly rotating flows with continuous throughput, but many studies of incompressible, single-phase flow in closed, rotating cylinders [3-13]. Of particular interest to this work are the studies concerning axial countercurrent flow and those few pertaining to mass inflow and outflow. The problem of axial countercurrent flow driven by an independently rotating endcap has been studied analytically, numerically and experimentally by many researchers, including Now at Exxon Research and Engineering, Florham Park, NJ, USA. 0045-7825/91/$03.50 © 1991 Elsevier Science Publishers B.V. All rights reserved
402
J.W. Frederick et at., A finite difference simulation of a centrifuge
Shadday et al. [12] and Ribando and Shadday [13]. These papers report flow fields predicted from finite difference simulation of the governing equations and have direct application to the present research. Both considered an incompressible fluid in a partially filled, rotating cylinder with an independently rotating endcap and a free surface, reported results obtained from a numerical code and compared them with experimental laser Doppler velocimeter measurements. The latter paper [13] used a code similar to that used in [12], but modified with an analytical matching condition at the endcaps in lieu of grid refinement to study the same types of flow at lower Ekman numbers. Some of the computational methods used in these two papers have been applied in the present research and are discussed more fully later. A schematic diagram of the experimental apparatus used by Krauss [14] in the flowvisualization studies which motivated the numerical studies reported here, is show in Fig. 1.
tI
Feed Removal
t
__~ _ _t_ 'l Baffle Plate
Free Surface
I I
Inviscid Region ',4'
I I I
I,I "1IIE''~ I
Feed Addition
Feed Disk
(Rotating)
I I I I I
I, I, I /3
Baffle Plate
E~3Jl'- -~11
(Independently
IF J:
/2
Rotating) r.-
Feed Removal
Fig. 1. Schematic of experimental apparatus used by Krauss [14]. Dotted lines indicate extent of various boundary layers.
J.W. Frederick et al., A finite difference simulation of a centrifuge
403
The dotted lines indicate the extent of boundary layers, which will be discussed in the next section. This apparatus was similar to that used by Popp [15] and Shadday et al. [12, 13], but included provisions for continuous fluid addition and removal. In the Krauss experiments a countercurrent secondary flow was driven either by over-rotation of the top lid or by insertion of a drag probe into the fluid near the bottom of the container. The entire process was made continuous by addition of feed at the inner free surface with fluid removal through the top and bottom lids. The feed was introduced through a perforated disk rotating at the same angular velocity as the cylinder. Thus the feed was accelerated to approximately the solid body velocity before contacting the main body of fluid. To reduce disruption cf flow, the fluid was removed from both the top and bottom lids through scoops located behind baffles.
2.
Theory
The unique features of rotating flows must be taken into account in the selection and implementation of any numerical scheme. Therefore, in this section a few of the basics of rotating flows are reviewed [16]. The governing equations for an incompressible, viscous, constant property fluid include conservation of mass:
V.q=O
(1)
and momentum: dq + ½V(q. q) + (V x q) x q + 2 / ] x q + K~ x ( ~ x r) at
1Vp+F
vVx(Vxq)
P
(2)
Here q is the velocity in the rotating reference frame, ~ is the angular velocity of the cylinder, and p and v are the density and kinematic viscosity of the fluid, respectively. The body force is assumed to be conservative, i.e., F = - V V , so that it and the centrifugal acceleration 1p(/~ x r). ( ~ x r) may be combined with the pressure P to form the reduced pressure [16]: p = P +
pV- _p(a x r). ( a x r).
(3)
These equations may be expressed in dimensionless form by scaling the length, time and velocity using the characteristic values L, ~-~ and U, respectively. In the system discussed first, L is the length of the cylinder, and U is the peripheral velocity of the independently rotating endcap measured in the rotating reference frame. Making the appropriate substitutions, the governing equations are V.q=0
and
O._._qq+tq.Vq+2,~xq=_Vp_EVxVxq,
(4),(5)
Ot
where q is now dimensionless, and the hat (^) denotes a unit vector. These equations are
J.W. Frederick et al., A finite difference simulation of a centrifuge
404
characterized by two dimensionless groups; the Ekman number E = v/J2L 2 ,
(6)
and the Rossby number e=U/OL.
(7)
The Ekman number compares the typical viscous force to the typical Coriolis force and is analogous to an inverse Reynolds number. The Rossby number measures the relative size of the convective acceleration and the Coriolis force, giving an estimate of the importance of non-linear terms. Examination of (5) reveals that the Ekman number multiplies the most highly differentiated terms. For large Ekman numbers (approaching unity), viscous forces are significant throughout the entire flow field. Conversely for small Ekman numbers, viscous effects are limited to thin boundary layers and inertial forces dominate the interior of the flow. In the present study, the Ekman number ranged from 10 -3 to 10 -1°, and so both regimes are of interest. Viscous layers may also occur in the flow interior to counteract any effects that tend to produce sharp or discontinuous velocity profiles [16]. The structure of the boundary layers depends upon the underlying force balance. For example, the flow in a cylinder with an independently rotating endcap has been examined extensively and will be referred to throughout this discussion. For this situation, the boundary layers develop as shown by the dotted lines in Fig. 1. The solution of the linearized equations, i.e., with the Rossby number approaching zero, as well as experimental evidence, reveals that momentum from the independently rotating endcap is transported to fluid through a viscous boundary layer of thickness of order E ~/2. Within this Ekman boundary layer the Coriolis force is balanced by viscous shear and the pressure is a constant of order O(E). In steady flow, fluid near the endcap is spun up to a greater angular velocity by direct viscous action. The increased Coriolis force overcomes the pressure gradient along the endcap and the fluid is propelled radially outward. To compensate for the mass flow in the boundary layer, a small normal flux of order E t/2 is required from the inviscid interior. At the rotor sidewall the flow is expelled into two vertical shear layers. The flow proceeds axially between the horizontal Ekman layers in a Stewartson boundary layer of thickness of order E 1/3. In the azimuthal plane, another layer of thickness of order E 1/4 matches the azimuthal fluid velocity generated by the differential disk rotation to the no-slip condition required at the rotor side wail. These boundary layers are all shown schematically in Fig. 1. These processes produce a secondary countercurrent flow imposed on the main solid body rotation, and it is flows of this type which are studied in the remainder of this paper. The horizontal Ekman boundary layers account for most of the radial transport of fluid while the vertical Stewartson boundary layers account for most of the axial transport. Between these boundary layers is an inviscid core characterized by a balance between the Coriolis force and the pressure gradient. For flows in which the deviation from solid body rotation is slight, i.e., corresponding to e ~ 0 ; steady ( O q / O t = 0); and inviscid (E = 0), the momentum equation reduces to × q = -vP.
(8)
J.W. Frederick et al., A finite difference simulation of a centrifuge
405
The curl of this expression combined with (4) gives the Taylor-Proudman Theorem: (~-V)q=O.
(9)
Under the restrictions given, the fluid velocity becomes two-dimensional, i.e., independent of the z-coordinate along the rotation axis. This requires every fluid particle in a vertical column to have the same velocity; hence each column moves as a single elongated vertical fluid element [16]. The above assumptions are very restrictive, however, and in actual experiments the Taylor-Proudman Theorem may not always strictly apply, although there will still be the tendency toward two-dimensionality, which grows stronger with decreasing Ekman number.
3. Numerical scheme
Most aspects of the numerical scheme employed in this study have been reported elsewhere [12, 13, 17]. Therefore, only a few highlights of the scheme and changes made to allow economical calculations at very low Ekman numbers (order 10-8-10 -9) will be discussed. For all the cases to be reported in the next section, the fluid was considered to be isothermal, incompressible and constant property. Under these restrictions the full, non-linear axisymmetric governing equations including continuity and conservation of radial, azimuthal and axial momentum are 1 Oru + Ow
r a--7
Tz =°'
O u + 1 Oru 2 + ~ _ O'-'-[ r O-----rOz
a-7
r
b'7
(10) 20+
o .... r
+ 20+r
+v Or
u=
--+ O'r r
Or
0z 2 / '
or +oz2/'
Ow + 1 a r u w + aw- = ap + v /' 1 a r aw + a2w] 0--7 r Or Oz - O"-z \ r -~r Or ~ / "
(11) (12) (13)
As before, these equations are written relative to a coordinate system rotating with the container, and the pressure indicated is the reduced pressure divided by the (constant) density. These equations may be non-dimensionalized as before to reveal the Ekman and Rossby numbers as parameters. ThE computational region is a vertical cross-section of the fluid-filled annulus (see Fig. 1) with azimuthal velocities perpendicular to this plane. A variable mesh is used so that grid points may be concentrated in areas having large velocity gradients. A typical computational mesh uses 20 radial grid points and 31 axial grid points. A staggered grid arrangement defines the axial and radial velocities at the appropriate edges of the pressure control volume, and the azimuthal velocity component is placed coincident with the radial velocity component. Previous algorithms [12] differed little from the original primitive variable algorithms of Harlow and Welch [18] and Hirt et al. [19] except in the treatment of the Coriolis terms. For instance, with the superscripts n and n + 1 indicating the present and advanced time levels,
J.W. Frederick et el,, A finite difference simulation of a centrifuge
406
respectively, the radial momentum equation (eq. (11)) is differenced as n -t- I
n
ui'J - ui'J At
n+
t! + 1
--Pi+l,j
(convection)','j + zatui.j + •
Ar
.
n
+ (vlscous)~ j.
•
(14)
With the azimuthal and radial velocities defined at the same point, a simple substitution allows fully-implicit treatment of the Coriolis terms, avoids a time step limitation of the form A t < 1 / 2 0 [20], and precludes excitation of inertial waves in the solution by the finite differencing [17]. The resulting form for the radial velocity and a similar form for the axial velocity are substituted into the continuity equation at the new time level to form a Poisson equation for pressure. With the pressure thus determined, velocities are updated and the process is repeated until convergence. For low Ekman numbers (10-5-10 -6) the explicit treatment of the viscous terms becomes a problem at the ends of the cylinder, when the grid is refined to resolve the Ekman layers. To alleviate the resulting severe time step restriction, equivalent matching conditions which have the same effect on the interior flow as would resolution of the boundary layers are used [16]. These take the form E Ij2 0
wr-
2r
Or r ( v r -
E I/2 0
Vr)'
wB-
2r
Or r ( ° B - VB)"
(15), (16)
Here v n and u r are the fluid velocities at the outer edge of the boundary layers; V8 and V r are the velocities of the differentially-rotating endcaps. In the case studied here, V r = 0. The implementation of these matching conditions into a finite difference algorithm is discussed in [13] and the theoretical basis is covered in [21]. For still lower Ekman numbers (order 10-7-10-9), and assuming the endwall boundary layer matching has been implemented as just discussed, the explicit treatment of radial diffusion terms at the sidewalls becomes the limiting factor in the timestep. Under these conditions and assuming that the non-linear terms are negligible (corresponding to low Rossby number), implementation of an Alternating Direction Implicit (ADI) scheme [22] for the diffusion terms becomes desirable. With the ADI scheme, allowable timesteps were found to be an order of magnitude greater than otherwise. The implementation of boundary conditions is straightforward. Slip and no-slip boundary conditions are accomplished by reflection with the appropriate sign change. Mass input or output through a vertical boundary involves specifying that prescribed value as the appropriate one of the four fluxes in the continuity statement at the adjacent pressure control volume. A prescribed mass flow through an endcap must be added to that determined from the matching conditions (eqs. (15) and (16)). The free surface was treated as free slip, impermeable and perfectly vertical.
4. Results - C l o s e d cylinder
The algorithm using ADI with the analytic matching condition at the end caps was used to simulate the flow in continuous flow centrifuges for several sample problems. For these linear
J.W. Frederick et al., A finite difference simulation of a centrifuge
407
cases, it was important to know the conditions under which they were applicable and the accuracy of the results as determined by comparison with experimental data. In all of the following cases involviag mass inflow and outflow, only the linear solution was obtained. It was therefore important to determine the point where non-linear effects become important and the linear code is no longer strictly applicable. Specifically, the Rossby number where the flow began to deviate significantly from the linear case was determined. This was accomplished by comparison of numerical and experimental results for the problem of a closed, rotating cylinder with an over-rotated endcap. Figure 2 shows results from the non-linear, finite difference code with resolved Ekman boundary layers. Figure 2(a) shows the computed streamlines when the top lid is over-rotated by 5%, and Figs. 2(b) and 2(c) show the measured [15] and computed radial distributions of the axial and azimuthal velocities. The case of an over-rotated lid with e = 0.01 is taken as the standard with which the non-linearity of the simulations at higher Rossby numbers will be compared. For purposes of comparison, the case of e = 0.01 (which is not shown) is assumed 0-61
I
|
I
I
I
Z/L=0.55
/ -I
• 0.o
o.5
1.o
(b)
R'
0.4t-
~
V,
1.0
0.032 Z/L=0.55
0.016 i
i W' o.o0o
Z,
0 0 0 0 0
-0.016
0,0 D.O
0.5
(a)
R'
1.0
-0.0320 •0
(c)
I
0.51
I
1.0
R'
Fig. 2. Streamlines (a) and measured and computed radial distributions of the azimuthal (b) and axial (c) velocity for a 5% over-rotated (e = 0.05) top lid. E = 2.2 x 10 -6. Streamline increment A~, = 0.75 x 10 -~ mZ/s. In this and all remaining figures the horizontal coordinate has been expanded by 4.75, velocities are normalized by the disk overspeed velocity (U = e ~ R ) , and the radial position is normalized by film thickness. Also the cylinder height is 18.86 cm and the film thickness is 1.98 cm.
408
J.W. Frederick et al., A finite difference simulation of a centrifuge
to be perfectly linear. Therefore the streamline increment in Fig. 2(a) has been set at half the difference between the maximum and minimum streamline values obtained from the case where e = 0.01. This case of e = 0.05 is nearly linear so there are only 10 streamlines. Figures 2(b) and 2(c) show that except near the free surface, there was qualitatively good agreement between the velocity results using the finite difference model and Popp's experimental results. The free-surface discrepancy was due to the limitations of the finite difference code in handling the free surface. As noted earlier, the free surface was assumed to be free slip, impermeable and perfectly vertical in the computer model. In the experiments there was some sag in the free surface (the centrifugal force was about 84 g's) and some additional radial displacement due to the effects of the over-rotated lid. The combination of these effects resulted in experimental measured velocities at the free surface that were slightly lower than those predicted by the numerical code. Simulation results for the case of a 10% over-rotated lid were only slightly different and are not presented. Figure 3(a) shows the computed streamlines for the case of e = 0.20 using the non-linear, resolved boundary layer, finite difference code. The radial distributions of the axial and azimuthal velocity profiles at a point slightly above the midplane for this same case and for
o,, i 0.6
V,
I
I
I
] I
0.2 0.0
I
3.0
(b) 1 0 ,~
I
I
0.5
1,0
R'
0,04
I
,,,
I
I
Z / L = 0,55
0.02
W ' 0.00
-0.02
_
......
0,20
----
0.10 Linear
.0
0,5
(a)
R'
1.0
-0.04 0 0 •
I
t 0.5
(c)
R'
I
.0
Fig. 3. Streamlines for a 20% (e = 0.20) over-rotated top lid (a) and computed radial distribution of the azimuthal (b) and axial (c) velocity for e--0.20 and 0.10 using non-linear code and from linear code. E = 2 . 2 x 10 -6, A~b = 3.0 X 10-" m~'/s.
409
J.W. Frederick et al., A finite difference simulation of a centrifuge
e = 0.10 computed by the non-linear code are shown in Figs. 3(b) and 3(c). Also shown are the results computed from the linear code that uses the analytic matching condition for the Ekman boundary layers. It is clear from these figures that the non-linear effects become important between e = 0.1 and 0.2. The streamlines in Fig. 3(a) show a significantly different flow pattern when compared with the streamlines shown in Fig. 2(a). There is some short-circuiting of the countercurrent flow near the top outer wall. In addition, the axial and azimuthal velocity profiles, when normalized by the disk peripheral velocity, show an increasing deviation from the linear case.
5. Results- Continuous flow centrifuge Having established in the previous section that non-linear effects become important somewhere between e = 0.1 and e = 0.2 for the case of an over-rotated endcap, this section looks at a number of cases involving inflow and outflow. These cases were chosen to illustrate the general characteristics of typical continuous flow centrifuges (Fig. 1). The centrifuge being modeled has the same dimensions as the cylinder modeled in the previous section. The first case involves mass inflow (with no axial or azimuthal shear stress) at the inner free surface and mass outflow through the top endcap. The streamlines and radial distribution of the axial and azimuthal velocity profiles computed for this case are shown in Fig. 4. Feed enters through a single grid point 1.0 cm in height midway between top and bottom endcaps and is removed through four radial grid points totaling 0.5 cm radius midway along the top endcap. The inlet
1.C
I
I
I
0
0.11 W'
.
I
0
V' -0.4
0.03
I
~
I
-.
-0.8-
-0,05 0.11" W'
Z,
V.
0.03
"0.
-0.8~
-0.05
-
0.11
"
1i
ZIL=0.15
W, 0.03
V,
-o.sh
-0.05 0.( u.u
~.0
"O-
0.0
I
I 0.5
I
1.0
0.0
I
I 0.5
I
1.0
(a) R' (b) R' (c) R' Fig. 4. Computed streamlines (a) and radial distributions of the axial (b) and azimuthal (c) velocity for mass inflow at the inner surface and mass outflow out the top. AO = 0.73 x 10-~ m2/s, e = 0.05.
J.W. Frederick et al., A finite difference simulation of a centrifuge
410
flow was 2.0 liter/min, corresponding to a residence time of 60 sec. The axial and azimuthal velocities were normalized in the same manner as in Fig. 2 for e = 0.05 (5% over-rotated lid). The streamline increment in Fig. 4(a) is the same as that used for the linear case of a 5% over-rotated lid. Figure 4(a) shows 12 streamlines indicating that the strength of this mass-driven flow is somewhat greater than that of a 5% over-rotated lid. The same figure also shows that most of the inflow is transported vertically through shear layers along the free surface to the Ekman boundary layers and out through the top endcap. The axial flow occurred in shear layers of order E 1/3 and most of the radial transport occurred in the E ~/~ Ekman boundary layers. Inlet flow diverted to the bottom portion of the centrifuge made its way to the top outlet through the middle of the inviscid core in a vertical column directly below the outlet. This flow pattern is very much different from that expected in a non-rotating system of similar geometry. It is, however, in agreement with the Taylor-Proudman theorem in that the flow tends to be two-dimensional in the fluid interior. The azimuthal velocity profiles in Fig. 4(c) shows that in the absence of an applied azimuthal shear stress supplied with the mass inflow, the main flow in the centrifuge is less than that corresponding to solid body rotation. Fluid exits the cylinder at solid body velocity at a larger radius than that of the fluid entering the cylinder at somewhat less than solid body velocity, resulting in the centrifuge acting as an angular momentum sink, i.e., angular momentum must be added to the fluid. Most of the angular momentum required by this action is transported to the fluid through viscous forces in the Ekman boundary layers, and thus requires negative azimuthal velocities. Because of the method used to add mass to the centrifuge, the absence of an applied azimuthal shear stress means that the fluid entered the cylinder at less than the solid body velocity in this simulation. 1.0
I
I I
0.03 ZIL= 0.85 W' -0.01
V'
I
I
°.O.I Z/L=O.85
-0.05
0.( OJ,
-o.olI.J °'°3l F~,~_.z,,L = O.S5
W'
Z,
O.
\
-0,0
Z/L=0.55
0.1
0.03
0.8
~ , . ~ Z/L=O.15 W' -0.01
~
(a)
'
'
R'
1.0
0,0
(b)
0.4
Z/L=0.15
-0,05 0.0!,0
V'
0.0 I
I 0.5
R'
I
I 1.0
0.0
(c)
I 0.5
I 1.0
R'
Fig. 5. Computed streamlines (a) and radial distributions of the axial (b) and azJ.muthal (c) velocity for an applied azimuthal shear stress. E = 2.2 × 10 -~, AqJ = 0.73 x 10 -6 m 2 / s , e = 0.05.
J.W. Frederick et al., A finite difference simulation of a centri[uge
411
The results displayed in Fig. 5 show the effects of an azimuthal shear stress applied to the inner free surface at the same point where feed was introduced in Fig. 4. This azimuthal shear stress is applied by prescribing an azimuthal velocity 5% faster than the solid body velocity at the boundary point grid cell along the inner free surface. In Fig. 5(a) the streamline increment is again 0.73 x 10 -6 m 2 / s and the flow produces only 5 streamlines indicating that the strength of this flow is much less than that of a 5% over-rotated lid. The azimuthal and axial velocities in Figs. 5(a) and 5(b) are normalized by the disk overspeed velocity U--= el-lR with e = 0.05. The flow is confined to a region near the free surface and does not penetrate deeply into the fluid layer. The flow in the meridionai plane is counterclockwise in the upper region and clockwise in the bottom half of the centrifuge. The direction of circulation would be reversed if the cylinder were spinning in the opposite direction or if a negative shear stress were applied. The additional angular momentum imposed on the system is dissipated through viscous forces in the Ekman boundary layers on both endcaps. The final set of plots in this series illustrates a combination of sources and sinks of mass and momentum. Figure 6 shows the computed results for a mass inflow of 2.01 [ min at solid body velocity at the inner free surface and mass outflow of 1.0 l/rain through the center of each endcap at solid body velocity. The bottom lid is under-~'otated by 5% to simulate the effect of a drag probe inserted into the fluid near the bottom of the container. This solution was obtained by a linear combination of four solutions: an over-rotated lid (Fig. 2), mass flow out the endcap (Fig. 4), azimuthal shear stress (Fig. 5) and axial shear stress. Mass flow through the bottom endcap was obtained by inverting the solution to mass outflow through the top lid
I
0.12 W' 0.04
I
I
=
I
O'il V, -0.5
~
-0.04-
o
v
-L
o.1~;
-0.5 -1. 0.0
~-
W, 0.04
Z/L=0.15
-
-0.04
u.u
,.O R'
ZI L = 0 . S J
-T
-o.o4.t
(a)
'
Z/L= 0"15~'~
V, -0.5 -1.0
0.0 (b)
I
"1.
O.l Tw
I
I
I
0.5 R'
I
1.0
0.0 (c)
I
0.5 R'
1.0
Fig. 6. Computed streamlines (a) and radial distributions of the axial (b) and azimuthal (c) velocity for mass inflow at the inner surface, mass outflow out the top and bottom with a 5% under-rotated (e = 0.05) bottom endcap. E = 2.2 x 1,0-~, AO = 0.73 x 10 -~ m2/s.
J.W. Frederick et al., A finite difference simulation of a centrifuge
412
using the coordinate transformation z' = 1 - z / L . The flow for an under-rotated lower lid was obtained by transforming the solution for the over-rotated top lid and changing the signs of all azimuthal velocities. Fluid entering at solid body velocity was obtained by adding the appropriate amount of azimuthal and axial shear stress at the free surface. Examination of Fig. 6 shows that the most of the feed flows upward in the Stewartson boundary layer to the top endcap resulting in good countercurrent flow in the upper half of the centrifuge. Flow on the bottom endcap is fed from the Stewartson E 1/3 boundary layer along the outer wall. Consequently a nearly stagnant region appears in the lower half of the centrifuge near the free surface. It should also be noted that the azimuthal flow is always slower than wheel flow away from the outer wall, which may indicate a decrease in the driving force for separation. The previous case studies were characterized by a relatively large Ekman number (2.2 x 10 -6 ) that allowed good illustration of the vertical boundary layer regions. Centrifuges of more practical interest are characterized by much smaller Ekman numbers, on the order of 10 -8 to 10 -9. This is caused both by higher speeds and longer rotors than used in the bench scale experiments referred to earlier. The main effect of the lower Ekman numbers is to compress the main axial and radial flows into thinner boundary layers. The computed streamlines and radial distributions of the axial and azimuthal velocity for the case of E = 2.2 x 10 -8 with a 5% over-rotated top lid are shown in Fig. 7. Thi~ case, like the others in this section, used the ADI with matching scheme. Full resolution of the Ekman layers at this small Ekman number would have been prohibitively expensive. Comparison of Fig. 7 with Fig. 2 shows the effect of a smaller Ekman number on the flow in identical cylinders. In
1.C
0.81-
0.04
W'
ZlL=0.85
V'
Z/L=0.55
I V'
0.00
-0,04 O.Oa~
Z'
W'
0.00
-0.04
V'
Z/L=0.15
%.
0.0(
"0.0~
(a)
R'
Z/L=0.85
0.4[-
01
Z/L=0.55
0.4
Z/L=O,15 0.4
0.0 I
1.0
I
0.8
0.¢ u.u
I
O.
o.o4~ W'
I
0.0
(b)
I 0.5
R'
I
I 1,0
0,0
(C)
I 0.5
I 1.0
R'
Fig. 7. Computed streamlines (a) and radial distributions of the axial (b) and azimuthal (c) velocity for a 5% over-rotated (e = 0.05) top lid. E = 2.2 x 10 -8, A~k = 0.83 x 10 -7 m2/s.
J.W. Frederick et al., A finite difference simulation of a centrifuge
413
addition to the flow being confined to thinner boundary layers, the amount of countercurrent circulation is much less. A streamline i n c r e m e n t of 0.83 x 10 -7 m2/s p r o d u c e d 10 streamlines indicating that the a m o u n t of flow was r e d u c e d by an o r d e r of magnitude.
6. Conclusion Five individual cases and one composite study have b e e n shown in an attempt to understand the basic flow characteristics of an incompressible fluid in a centrifuge. A linear version of the governing equations was m o d e l e d to allow study of very small E k m a n n u m b e r flows and to allow easy combination of the resulting solutions. A small n u m b e r of distinct problems were t h e r e f o r e simulated to enable a multitude of general solutions to be obtained easily without using large amounts of c o m p u t e r resources. From these studies the basic nature of strongly rotating flows have b e e n illustrated using the streamline presentations. T h e thin boundary layer structure along the endcaps and sidewalls, and the t e n d e n c y towards two-dimensionality have b e e n clearly shown. In addition, these two characteristics have been shown to play a very i m p o r t a n t role in understanding the structure of the flow in the mass throughput cases.
References [1] H.G. Wood and J.B. Morton, Onsager's pancake approximation for the fluid dynamics of a gas centrifuge, J. Fluid Mech. 101 (1980) 1-31. [2] H.G. Wood and G. Sanders,. Rotating compressible flows with internal sources and sinks, J. Fluid Mech. 127 (1983) 299-313. [3] F. Bein and S.S. Penner, Velocity profiles in steady and unsteady rotating flows for a finite cylindrical geometry, Phys. Fluids 13 (1970) 1661-1665. [4] H.P. Pao, Numerical computation of a confined rotating flow, ASME Trans. 6 (1970) 480-487. [5] H.P. Pao, Numerical solution of the Navier-Stokes equations for flows in the disk-cylinder system, Phys. Fluids 15 (1972) 4-11. [6] W.B. Watkins and R.G. Hussey, Spin-up from rest in a cylinder, Phys. Fluids 20 (1977) 1596-1604. [7] E.R. Benton, Vorticity dynamics in spin-up from rest, Phys. Fluids 22 (1979) 1250-1251. [8] C.W. Kitchens Jr., Navier-Stokes solutions for spin-up in a filled cylinder, AIAA J. 18 (1980) 929-934. [9] M. Bertela and F. Gori, Laminar flow in a cylindrical container with a rotating cover, J. Fluids Engrg. 104 (1982) 31-39. [lfi] J.M. Hyun, Spin-up from rest in a differentially rotating cylinder, AIAA J. 21 (1982) 1278-1282. [11] J.M. Hyun, F. Leslie, W.M. Fowlis and A. Warn-Varaas, Numerical solutions for spin-up from rest in a cylinder, J. Fluid Mech. 127 (1983) 263-281. [12] M.A. Shadday, R.J. Ribando and J.J. Kauzlarich, Flow of an incompressible fluid in a partially filled, rapidly rotating cylinder with a differentially rotating endcap, J. Fluid Mech. 130 (1983) 203-218. [13] R.J. Ribando and M.A. Shadday, The Ekman matching condition in a partially filled, rapidly rotating cylinder, J. Comp. Phys. 53 (1984) 266-288. [14] R. Krauss, University of Virginia, Written communication (1986). [15] C.G. Popp, Countercurrent velocity measurements from the rotating reference frame, MS Thesis, University of Virginia, 1987. [16] H.E Greenspa- - r ~ Theory of Rotating Fluids (Cambridge Univ. Press, New York, NY, 1968). [17] G.E. Heuser, RJ. Ribando and H.G. Wood, III, A numerical simulation of inertial waves in a rotating fluid, Comput. Methods Appl. Mech. Engrg. 57 (1986) 207-222.
414
J.W. Frederick et al., A finite difference simulation of a centrifuge
[18] F.H. Harlow and J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (1965) 2182-2189. [19] C.W. Hirt, B.D. Nichols and N.C. Romero, SOLA- A numerical solution algorithm for transient fluid flows, LA-5852 (1975) Los Alamos National Laboratory. [20] R.D. Beardsley, K.D. Saunders, A.C. Warn-Varnas and J.M. Harding, An experimental and numerical study of the secular spin-up of a thermal stratified rotating fluid, 3. Fluid Mech. 93 (1979) 161-164. [21] H.P. Greenspan, Simulation of countercurrent flow in a gas centrifuge, UVA-ER-746-82U, University of Virginia, 1982. [22] P.3. Roache, Computational Fluid Dynamics (Hermosa, Albuquerque, 1976).
A finite difference simulation of a continuous flow centrifuge J.W. Frederick ~, R.J. Ribando and H.G. Wood Department of Mechanical and Aerospace Engineering, School of Engineering and Applied Science, University of Virginia, Charlottesville, VA 22901, USA Received 8 March 1990
The secondary flow of an incompressible liquid in a continuous process, vertical tubular centrifuge is studied using finite difference solutions to the governing equations. In particular the secondary flow associated with the addition and removal of mass or momentum from the container is computed and analyzed. For parameters typical of the liquid centrifuges studied here, the governing equations are linear. Thus various feed and withdrawl configurations may be analyzed using superposition of a few solutions to the governing equations. Even aside from the secondary flow induced by the density difference between the species being separated, the addition and withdrawl of fluid is found to have a major impact on the secondary flow pattern.
1. Introduction
In this paper a hypothetical tubular centrifuge is analyzed using finite difference solutions to the governing equations for a homogeneous fluid. While a real centrifuge obviously entails a density difference either between two fluids or between a fluid and solid particles dispersed in it, this study is restricted to species of negligible property difference so that the dynamical effects of perturbations due to sources and sinks of mass and momentum can be studied. For most cases of interest it can be shown that the governing equations are linear to good approximation (as measured by the Rossby number, a non-dimensional parameter to be discussed later), and thus superposition of separate simple effects can be used as a tool in calculating and analyzing more complicated flows. This approach is used extensively in the design and analysis of gas centrifuges for uranium enrichment [1]. In gas centrifuge analysis the effects of withdrawl scoops and mass injection have been simulated using point or distributed sources of mass and momentum [2]. There have been very few studies other than in the gas centrifuge literature involving strongly rotating flows with continuous throughput, but many studies of incompressible, single-phase flow in closed, rotating cylinders [3-13]. Of particular interest to this work are the studies concerning axial countercurrent flow and those few pertaining to mass inflow and outflow. The problem of axial countercurrent flow driven by an independently rotating endcap has been studied analytically, numerically and experimentally by many researchers, including Now at Exxon Research and Engineering, Florham Park, NJ, USA. 0045-7825/91/$03.50 © 1991 Elsevier Science Publishers B.V. All rights reserved
402
J.W. Frederick et at., A finite difference simulation of a centrifuge
Shadday et al. [12] and Ribando and Shadday [13]. These papers report flow fields predicted from finite difference simulation of the governing equations and have direct application to the present research. Both considered an incompressible fluid in a partially filled, rotating cylinder with an independently rotating endcap and a free surface, reported results obtained from a numerical code and compared them with experimental laser Doppler velocimeter measurements. The latter paper [13] used a code similar to that used in [12], but modified with an analytical matching condition at the endcaps in lieu of grid refinement to study the same types of flow at lower Ekman numbers. Some of the computational methods used in these two papers have been applied in the present research and are discussed more fully later. A schematic diagram of the experimental apparatus used by Krauss [14] in the flowvisualization studies which motivated the numerical studies reported here, is show in Fig. 1.
tI
Feed Removal
t
__~ _ _t_ 'l Baffle Plate
Free Surface
I I
Inviscid Region ',4'
I I I
I,I "1IIE''~ I
Feed Addition
Feed Disk
(Rotating)
I I I I I
I, I, I /3
Baffle Plate
E~3Jl'- -~11
(Independently
IF J:
/2
Rotating) r.-
Feed Removal
Fig. 1. Schematic of experimental apparatus used by Krauss [14]. Dotted lines indicate extent of various boundary layers.
J.W. Frederick et al., A finite difference simulation of a centrifuge
403
The dotted lines indicate the extent of boundary layers, which will be discussed in the next section. This apparatus was similar to that used by Popp [15] and Shadday et al. [12, 13], but included provisions for continuous fluid addition and removal. In the Krauss experiments a countercurrent secondary flow was driven either by over-rotation of the top lid or by insertion of a drag probe into the fluid near the bottom of the container. The entire process was made continuous by addition of feed at the inner free surface with fluid removal through the top and bottom lids. The feed was introduced through a perforated disk rotating at the same angular velocity as the cylinder. Thus the feed was accelerated to approximately the solid body velocity before contacting the main body of fluid. To reduce disruption cf flow, the fluid was removed from both the top and bottom lids through scoops located behind baffles.
2.
Theory
The unique features of rotating flows must be taken into account in the selection and implementation of any numerical scheme. Therefore, in this section a few of the basics of rotating flows are reviewed [16]. The governing equations for an incompressible, viscous, constant property fluid include conservation of mass:
V.q=O
(1)
and momentum: dq + ½V(q. q) + (V x q) x q + 2 / ] x q + K~ x ( ~ x r) at
1Vp+F
vVx(Vxq)
P
(2)
Here q is the velocity in the rotating reference frame, ~ is the angular velocity of the cylinder, and p and v are the density and kinematic viscosity of the fluid, respectively. The body force is assumed to be conservative, i.e., F = - V V , so that it and the centrifugal acceleration 1p(/~ x r). ( ~ x r) may be combined with the pressure P to form the reduced pressure [16]: p = P +
pV- _p(a x r). ( a x r).
(3)
These equations may be expressed in dimensionless form by scaling the length, time and velocity using the characteristic values L, ~-~ and U, respectively. In the system discussed first, L is the length of the cylinder, and U is the peripheral velocity of the independently rotating endcap measured in the rotating reference frame. Making the appropriate substitutions, the governing equations are V.q=0
and
O._._qq+tq.Vq+2,~xq=_Vp_EVxVxq,
(4),(5)
Ot
where q is now dimensionless, and the hat (^) denotes a unit vector. These equations are
J.W. Frederick et al., A finite difference simulation of a centrifuge
404
characterized by two dimensionless groups; the Ekman number E = v/J2L 2 ,
(6)
and the Rossby number e=U/OL.
(7)
The Ekman number compares the typical viscous force to the typical Coriolis force and is analogous to an inverse Reynolds number. The Rossby number measures the relative size of the convective acceleration and the Coriolis force, giving an estimate of the importance of non-linear terms. Examination of (5) reveals that the Ekman number multiplies the most highly differentiated terms. For large Ekman numbers (approaching unity), viscous forces are significant throughout the entire flow field. Conversely for small Ekman numbers, viscous effects are limited to thin boundary layers and inertial forces dominate the interior of the flow. In the present study, the Ekman number ranged from 10 -3 to 10 -1°, and so both regimes are of interest. Viscous layers may also occur in the flow interior to counteract any effects that tend to produce sharp or discontinuous velocity profiles [16]. The structure of the boundary layers depends upon the underlying force balance. For example, the flow in a cylinder with an independently rotating endcap has been examined extensively and will be referred to throughout this discussion. For this situation, the boundary layers develop as shown by the dotted lines in Fig. 1. The solution of the linearized equations, i.e., with the Rossby number approaching zero, as well as experimental evidence, reveals that momentum from the independently rotating endcap is transported to fluid through a viscous boundary layer of thickness of order E ~/2. Within this Ekman boundary layer the Coriolis force is balanced by viscous shear and the pressure is a constant of order O(E). In steady flow, fluid near the endcap is spun up to a greater angular velocity by direct viscous action. The increased Coriolis force overcomes the pressure gradient along the endcap and the fluid is propelled radially outward. To compensate for the mass flow in the boundary layer, a small normal flux of order E t/2 is required from the inviscid interior. At the rotor sidewall the flow is expelled into two vertical shear layers. The flow proceeds axially between the horizontal Ekman layers in a Stewartson boundary layer of thickness of order E 1/3. In the azimuthal plane, another layer of thickness of order E 1/4 matches the azimuthal fluid velocity generated by the differential disk rotation to the no-slip condition required at the rotor side wail. These boundary layers are all shown schematically in Fig. 1. These processes produce a secondary countercurrent flow imposed on the main solid body rotation, and it is flows of this type which are studied in the remainder of this paper. The horizontal Ekman boundary layers account for most of the radial transport of fluid while the vertical Stewartson boundary layers account for most of the axial transport. Between these boundary layers is an inviscid core characterized by a balance between the Coriolis force and the pressure gradient. For flows in which the deviation from solid body rotation is slight, i.e., corresponding to e ~ 0 ; steady ( O q / O t = 0); and inviscid (E = 0), the momentum equation reduces to × q = -vP.
(8)
J.W. Frederick et al., A finite difference simulation of a centrifuge
405
The curl of this expression combined with (4) gives the Taylor-Proudman Theorem: (~-V)q=O.
(9)
Under the restrictions given, the fluid velocity becomes two-dimensional, i.e., independent of the z-coordinate along the rotation axis. This requires every fluid particle in a vertical column to have the same velocity; hence each column moves as a single elongated vertical fluid element [16]. The above assumptions are very restrictive, however, and in actual experiments the Taylor-Proudman Theorem may not always strictly apply, although there will still be the tendency toward two-dimensionality, which grows stronger with decreasing Ekman number.
3. Numerical scheme
Most aspects of the numerical scheme employed in this study have been reported elsewhere [12, 13, 17]. Therefore, only a few highlights of the scheme and changes made to allow economical calculations at very low Ekman numbers (order 10-8-10 -9) will be discussed. For all the cases to be reported in the next section, the fluid was considered to be isothermal, incompressible and constant property. Under these restrictions the full, non-linear axisymmetric governing equations including continuity and conservation of radial, azimuthal and axial momentum are 1 Oru + Ow
r a--7
Tz =°'
O u + 1 Oru 2 + ~ _ O'-'-[ r O-----rOz
a-7
r
b'7
(10) 20+
o .... r
+ 20+r
+v Or
u=
--+ O'r r
Or
0z 2 / '
or +oz2/'
Ow + 1 a r u w + aw- = ap + v /' 1 a r aw + a2w] 0--7 r Or Oz - O"-z \ r -~r Or ~ / "
(11) (12) (13)
As before, these equations are written relative to a coordinate system rotating with the container, and the pressure indicated is the reduced pressure divided by the (constant) density. These equations may be non-dimensionalized as before to reveal the Ekman and Rossby numbers as parameters. ThE computational region is a vertical cross-section of the fluid-filled annulus (see Fig. 1) with azimuthal velocities perpendicular to this plane. A variable mesh is used so that grid points may be concentrated in areas having large velocity gradients. A typical computational mesh uses 20 radial grid points and 31 axial grid points. A staggered grid arrangement defines the axial and radial velocities at the appropriate edges of the pressure control volume, and the azimuthal velocity component is placed coincident with the radial velocity component. Previous algorithms [12] differed little from the original primitive variable algorithms of Harlow and Welch [18] and Hirt et al. [19] except in the treatment of the Coriolis terms. For instance, with the superscripts n and n + 1 indicating the present and advanced time levels,
J.W. Frederick et el,, A finite difference simulation of a centrifuge
406
respectively, the radial momentum equation (eq. (11)) is differenced as n -t- I
n
ui'J - ui'J At
n+
t! + 1
--Pi+l,j
(convection)','j + zatui.j + •
Ar
.
n
+ (vlscous)~ j.
•
(14)
With the azimuthal and radial velocities defined at the same point, a simple substitution allows fully-implicit treatment of the Coriolis terms, avoids a time step limitation of the form A t < 1 / 2 0 [20], and precludes excitation of inertial waves in the solution by the finite differencing [17]. The resulting form for the radial velocity and a similar form for the axial velocity are substituted into the continuity equation at the new time level to form a Poisson equation for pressure. With the pressure thus determined, velocities are updated and the process is repeated until convergence. For low Ekman numbers (10-5-10 -6) the explicit treatment of the viscous terms becomes a problem at the ends of the cylinder, when the grid is refined to resolve the Ekman layers. To alleviate the resulting severe time step restriction, equivalent matching conditions which have the same effect on the interior flow as would resolution of the boundary layers are used [16]. These take the form E Ij2 0
wr-
2r
Or r ( v r -
E I/2 0
Vr)'
wB-
2r
Or r ( ° B - VB)"
(15), (16)
Here v n and u r are the fluid velocities at the outer edge of the boundary layers; V8 and V r are the velocities of the differentially-rotating endcaps. In the case studied here, V r = 0. The implementation of these matching conditions into a finite difference algorithm is discussed in [13] and the theoretical basis is covered in [21]. For still lower Ekman numbers (order 10-7-10-9), and assuming the endwall boundary layer matching has been implemented as just discussed, the explicit treatment of radial diffusion terms at the sidewalls becomes the limiting factor in the timestep. Under these conditions and assuming that the non-linear terms are negligible (corresponding to low Rossby number), implementation of an Alternating Direction Implicit (ADI) scheme [22] for the diffusion terms becomes desirable. With the ADI scheme, allowable timesteps were found to be an order of magnitude greater than otherwise. The implementation of boundary conditions is straightforward. Slip and no-slip boundary conditions are accomplished by reflection with the appropriate sign change. Mass input or output through a vertical boundary involves specifying that prescribed value as the appropriate one of the four fluxes in the continuity statement at the adjacent pressure control volume. A prescribed mass flow through an endcap must be added to that determined from the matching conditions (eqs. (15) and (16)). The free surface was treated as free slip, impermeable and perfectly vertical.
4. Results - C l o s e d cylinder
The algorithm using ADI with the analytic matching condition at the end caps was used to simulate the flow in continuous flow centrifuges for several sample problems. For these linear
J.W. Frederick et al., A finite difference simulation of a centrifuge
407
cases, it was important to know the conditions under which they were applicable and the accuracy of the results as determined by comparison with experimental data. In all of the following cases involviag mass inflow and outflow, only the linear solution was obtained. It was therefore important to determine the point where non-linear effects become important and the linear code is no longer strictly applicable. Specifically, the Rossby number where the flow began to deviate significantly from the linear case was determined. This was accomplished by comparison of numerical and experimental results for the problem of a closed, rotating cylinder with an over-rotated endcap. Figure 2 shows results from the non-linear, finite difference code with resolved Ekman boundary layers. Figure 2(a) shows the computed streamlines when the top lid is over-rotated by 5%, and Figs. 2(b) and 2(c) show the measured [15] and computed radial distributions of the axial and azimuthal velocities. The case of an over-rotated lid with e = 0.01 is taken as the standard with which the non-linearity of the simulations at higher Rossby numbers will be compared. For purposes of comparison, the case of e = 0.01 (which is not shown) is assumed 0-61
I
|
I
I
I
Z/L=0.55
/ -I
• 0.o
o.5
1.o
(b)
R'
0.4t-
~
V,
1.0
0.032 Z/L=0.55
0.016 i
i W' o.o0o
Z,
0 0 0 0 0
-0.016
0,0 D.O
0.5
(a)
R'
1.0
-0.0320 •0
(c)
I
0.51
I
1.0
R'
Fig. 2. Streamlines (a) and measured and computed radial distributions of the azimuthal (b) and axial (c) velocity for a 5% over-rotated (e = 0.05) top lid. E = 2.2 x 10 -6. Streamline increment A~, = 0.75 x 10 -~ mZ/s. In this and all remaining figures the horizontal coordinate has been expanded by 4.75, velocities are normalized by the disk overspeed velocity (U = e ~ R ) , and the radial position is normalized by film thickness. Also the cylinder height is 18.86 cm and the film thickness is 1.98 cm.
408
J.W. Frederick et al., A finite difference simulation of a centrifuge
to be perfectly linear. Therefore the streamline increment in Fig. 2(a) has been set at half the difference between the maximum and minimum streamline values obtained from the case where e = 0.01. This case of e = 0.05 is nearly linear so there are only 10 streamlines. Figures 2(b) and 2(c) show that except near the free surface, there was qualitatively good agreement between the velocity results using the finite difference model and Popp's experimental results. The free-surface discrepancy was due to the limitations of the finite difference code in handling the free surface. As noted earlier, the free surface was assumed to be free slip, impermeable and perfectly vertical in the computer model. In the experiments there was some sag in the free surface (the centrifugal force was about 84 g's) and some additional radial displacement due to the effects of the over-rotated lid. The combination of these effects resulted in experimental measured velocities at the free surface that were slightly lower than those predicted by the numerical code. Simulation results for the case of a 10% over-rotated lid were only slightly different and are not presented. Figure 3(a) shows the computed streamlines for the case of e = 0.20 using the non-linear, resolved boundary layer, finite difference code. The radial distributions of the axial and azimuthal velocity profiles at a point slightly above the midplane for this same case and for
o,, i 0.6
V,
I
I
I
] I
0.2 0.0
I
3.0
(b) 1 0 ,~
I
I
0.5
1,0
R'
0,04
I
,,,
I
I
Z / L = 0,55
0.02
W ' 0.00
-0.02
_
......
0,20
----
0.10 Linear
.0
0,5
(a)
R'
1.0
-0.04 0 0 •
I
t 0.5
(c)
R'
I
.0
Fig. 3. Streamlines for a 20% (e = 0.20) over-rotated top lid (a) and computed radial distribution of the azimuthal (b) and axial (c) velocity for e--0.20 and 0.10 using non-linear code and from linear code. E = 2 . 2 x 10 -6, A~b = 3.0 X 10-" m~'/s.
409
J.W. Frederick et al., A finite difference simulation of a centrifuge
e = 0.10 computed by the non-linear code are shown in Figs. 3(b) and 3(c). Also shown are the results computed from the linear code that uses the analytic matching condition for the Ekman boundary layers. It is clear from these figures that the non-linear effects become important between e = 0.1 and 0.2. The streamlines in Fig. 3(a) show a significantly different flow pattern when compared with the streamlines shown in Fig. 2(a). There is some short-circuiting of the countercurrent flow near the top outer wall. In addition, the axial and azimuthal velocity profiles, when normalized by the disk peripheral velocity, show an increasing deviation from the linear case.
5. Results- Continuous flow centrifuge Having established in the previous section that non-linear effects become important somewhere between e = 0.1 and e = 0.2 for the case of an over-rotated endcap, this section looks at a number of cases involving inflow and outflow. These cases were chosen to illustrate the general characteristics of typical continuous flow centrifuges (Fig. 1). The centrifuge being modeled has the same dimensions as the cylinder modeled in the previous section. The first case involves mass inflow (with no axial or azimuthal shear stress) at the inner free surface and mass outflow through the top endcap. The streamlines and radial distribution of the axial and azimuthal velocity profiles computed for this case are shown in Fig. 4. Feed enters through a single grid point 1.0 cm in height midway between top and bottom endcaps and is removed through four radial grid points totaling 0.5 cm radius midway along the top endcap. The inlet
1.C
I
I
I
0
0.11 W'
.
I
0
V' -0.4
0.03
I
~
I
-.
-0.8-
-0,05 0.11" W'
Z,
V.
0.03
"0.
-0.8~
-0.05
-
0.11
"
1i
ZIL=0.15
W, 0.03
V,
-o.sh
-0.05 0.( u.u
~.0
"O-
0.0
I
I 0.5
I
1.0
0.0
I
I 0.5
I
1.0
(a) R' (b) R' (c) R' Fig. 4. Computed streamlines (a) and radial distributions of the axial (b) and azimuthal (c) velocity for mass inflow at the inner surface and mass outflow out the top. AO = 0.73 x 10-~ m2/s, e = 0.05.
J.W. Frederick et al., A finite difference simulation of a centrifuge
410
flow was 2.0 liter/min, corresponding to a residence time of 60 sec. The axial and azimuthal velocities were normalized in the same manner as in Fig. 2 for e = 0.05 (5% over-rotated lid). The streamline increment in Fig. 4(a) is the same as that used for the linear case of a 5% over-rotated lid. Figure 4(a) shows 12 streamlines indicating that the strength of this mass-driven flow is somewhat greater than that of a 5% over-rotated lid. The same figure also shows that most of the inflow is transported vertically through shear layers along the free surface to the Ekman boundary layers and out through the top endcap. The axial flow occurred in shear layers of order E 1/3 and most of the radial transport occurred in the E ~/~ Ekman boundary layers. Inlet flow diverted to the bottom portion of the centrifuge made its way to the top outlet through the middle of the inviscid core in a vertical column directly below the outlet. This flow pattern is very much different from that expected in a non-rotating system of similar geometry. It is, however, in agreement with the Taylor-Proudman theorem in that the flow tends to be two-dimensional in the fluid interior. The azimuthal velocity profiles in Fig. 4(c) shows that in the absence of an applied azimuthal shear stress supplied with the mass inflow, the main flow in the centrifuge is less than that corresponding to solid body rotation. Fluid exits the cylinder at solid body velocity at a larger radius than that of the fluid entering the cylinder at somewhat less than solid body velocity, resulting in the centrifuge acting as an angular momentum sink, i.e., angular momentum must be added to the fluid. Most of the angular momentum required by this action is transported to the fluid through viscous forces in the Ekman boundary layers, and thus requires negative azimuthal velocities. Because of the method used to add mass to the centrifuge, the absence of an applied azimuthal shear stress means that the fluid entered the cylinder at less than the solid body velocity in this simulation. 1.0
I
I I
0.03 ZIL= 0.85 W' -0.01
V'
I
I
°.O.I Z/L=O.85
-0.05
0.( OJ,
-o.olI.J °'°3l F~,~_.z,,L = O.S5
W'
Z,
O.
\
-0,0
Z/L=0.55
0.1
0.03
0.8
~ , . ~ Z/L=O.15 W' -0.01
~
(a)
'
'
R'
1.0
0,0
(b)
0.4
Z/L=0.15
-0,05 0.0!,0
V'
0.0 I
I 0.5
R'
I
I 1.0
0.0
(c)
I 0.5
I 1.0
R'
Fig. 5. Computed streamlines (a) and radial distributions of the axial (b) and azJ.muthal (c) velocity for an applied azimuthal shear stress. E = 2.2 × 10 -~, AqJ = 0.73 x 10 -6 m 2 / s , e = 0.05.
J.W. Frederick et al., A finite difference simulation of a centri[uge
411
The results displayed in Fig. 5 show the effects of an azimuthal shear stress applied to the inner free surface at the same point where feed was introduced in Fig. 4. This azimuthal shear stress is applied by prescribing an azimuthal velocity 5% faster than the solid body velocity at the boundary point grid cell along the inner free surface. In Fig. 5(a) the streamline increment is again 0.73 x 10 -6 m 2 / s and the flow produces only 5 streamlines indicating that the strength of this flow is much less than that of a 5% over-rotated lid. The azimuthal and axial velocities in Figs. 5(a) and 5(b) are normalized by the disk overspeed velocity U--= el-lR with e = 0.05. The flow is confined to a region near the free surface and does not penetrate deeply into the fluid layer. The flow in the meridionai plane is counterclockwise in the upper region and clockwise in the bottom half of the centrifuge. The direction of circulation would be reversed if the cylinder were spinning in the opposite direction or if a negative shear stress were applied. The additional angular momentum imposed on the system is dissipated through viscous forces in the Ekman boundary layers on both endcaps. The final set of plots in this series illustrates a combination of sources and sinks of mass and momentum. Figure 6 shows the computed results for a mass inflow of 2.01 [ min at solid body velocity at the inner free surface and mass outflow of 1.0 l/rain through the center of each endcap at solid body velocity. The bottom lid is under-~'otated by 5% to simulate the effect of a drag probe inserted into the fluid near the bottom of the container. This solution was obtained by a linear combination of four solutions: an over-rotated lid (Fig. 2), mass flow out the endcap (Fig. 4), azimuthal shear stress (Fig. 5) and axial shear stress. Mass flow through the bottom endcap was obtained by inverting the solution to mass outflow through the top lid
I
0.12 W' 0.04
I
I
=
I
O'il V, -0.5
~
-0.04-
o
v
-L
o.1~;
-0.5 -1. 0.0
~-
W, 0.04
Z/L=0.15
-
-0.04
u.u
,.O R'
ZI L = 0 . S J
-T
-o.o4.t
(a)
'
Z/L= 0"15~'~
V, -0.5 -1.0
0.0 (b)
I
"1.
O.l Tw
I
I
I
0.5 R'
I
1.0
0.0 (c)
I
0.5 R'
1.0
Fig. 6. Computed streamlines (a) and radial distributions of the axial (b) and azimuthal (c) velocity for mass inflow at the inner surface, mass outflow out the top and bottom with a 5% under-rotated (e = 0.05) bottom endcap. E = 2.2 x 1,0-~, AO = 0.73 x 10 -~ m2/s.
J.W. Frederick et al., A finite difference simulation of a centrifuge
412
using the coordinate transformation z' = 1 - z / L . The flow for an under-rotated lower lid was obtained by transforming the solution for the over-rotated top lid and changing the signs of all azimuthal velocities. Fluid entering at solid body velocity was obtained by adding the appropriate amount of azimuthal and axial shear stress at the free surface. Examination of Fig. 6 shows that the most of the feed flows upward in the Stewartson boundary layer to the top endcap resulting in good countercurrent flow in the upper half of the centrifuge. Flow on the bottom endcap is fed from the Stewartson E 1/3 boundary layer along the outer wall. Consequently a nearly stagnant region appears in the lower half of the centrifuge near the free surface. It should also be noted that the azimuthal flow is always slower than wheel flow away from the outer wall, which may indicate a decrease in the driving force for separation. The previous case studies were characterized by a relatively large Ekman number (2.2 x 10 -6 ) that allowed good illustration of the vertical boundary layer regions. Centrifuges of more practical interest are characterized by much smaller Ekman numbers, on the order of 10 -8 to 10 -9. This is caused both by higher speeds and longer rotors than used in the bench scale experiments referred to earlier. The main effect of the lower Ekman numbers is to compress the main axial and radial flows into thinner boundary layers. The computed streamlines and radial distributions of the axial and azimuthal velocity for the case of E = 2.2 x 10 -8 with a 5% over-rotated top lid are shown in Fig. 7. Thi~ case, like the others in this section, used the ADI with matching scheme. Full resolution of the Ekman layers at this small Ekman number would have been prohibitively expensive. Comparison of Fig. 7 with Fig. 2 shows the effect of a smaller Ekman number on the flow in identical cylinders. In
1.C
0.81-
0.04
W'
ZlL=0.85
V'
Z/L=0.55
I V'
0.00
-0,04 O.Oa~
Z'
W'
0.00
-0.04
V'
Z/L=0.15
%.
0.0(
"0.0~
(a)
R'
Z/L=0.85
0.4[-
01
Z/L=0.55
0.4
Z/L=O,15 0.4
0.0 I
1.0
I
0.8
0.¢ u.u
I
O.
o.o4~ W'
I
0.0
(b)
I 0.5
R'
I
I 1,0
0,0
(C)
I 0.5
I 1.0
R'
Fig. 7. Computed streamlines (a) and radial distributions of the axial (b) and azimuthal (c) velocity for a 5% over-rotated (e = 0.05) top lid. E = 2.2 x 10 -8, A~k = 0.83 x 10 -7 m2/s.
J.W. Frederick et al., A finite difference simulation of a centrifuge
413
addition to the flow being confined to thinner boundary layers, the amount of countercurrent circulation is much less. A streamline i n c r e m e n t of 0.83 x 10 -7 m2/s p r o d u c e d 10 streamlines indicating that the a m o u n t of flow was r e d u c e d by an o r d e r of magnitude.
6. Conclusion Five individual cases and one composite study have b e e n shown in an attempt to understand the basic flow characteristics of an incompressible fluid in a centrifuge. A linear version of the governing equations was m o d e l e d to allow study of very small E k m a n n u m b e r flows and to allow easy combination of the resulting solutions. A small n u m b e r of distinct problems were t h e r e f o r e simulated to enable a multitude of general solutions to be obtained easily without using large amounts of c o m p u t e r resources. From these studies the basic nature of strongly rotating flows have b e e n illustrated using the streamline presentations. T h e thin boundary layer structure along the endcaps and sidewalls, and the t e n d e n c y towards two-dimensionality have b e e n clearly shown. In addition, these two characteristics have been shown to play a very i m p o r t a n t role in understanding the structure of the flow in the mass throughput cases.
References [1] H.G. Wood and J.B. Morton, Onsager's pancake approximation for the fluid dynamics of a gas centrifuge, J. Fluid Mech. 101 (1980) 1-31. [2] H.G. Wood and G. Sanders,. Rotating compressible flows with internal sources and sinks, J. Fluid Mech. 127 (1983) 299-313. [3] F. Bein and S.S. Penner, Velocity profiles in steady and unsteady rotating flows for a finite cylindrical geometry, Phys. Fluids 13 (1970) 1661-1665. [4] H.P. Pao, Numerical computation of a confined rotating flow, ASME Trans. 6 (1970) 480-487. [5] H.P. Pao, Numerical solution of the Navier-Stokes equations for flows in the disk-cylinder system, Phys. Fluids 15 (1972) 4-11. [6] W.B. Watkins and R.G. Hussey, Spin-up from rest in a cylinder, Phys. Fluids 20 (1977) 1596-1604. [7] E.R. Benton, Vorticity dynamics in spin-up from rest, Phys. Fluids 22 (1979) 1250-1251. [8] C.W. Kitchens Jr., Navier-Stokes solutions for spin-up in a filled cylinder, AIAA J. 18 (1980) 929-934. [9] M. Bertela and F. Gori, Laminar flow in a cylindrical container with a rotating cover, J. Fluids Engrg. 104 (1982) 31-39. [lfi] J.M. Hyun, Spin-up from rest in a differentially rotating cylinder, AIAA J. 21 (1982) 1278-1282. [11] J.M. Hyun, F. Leslie, W.M. Fowlis and A. Warn-Varaas, Numerical solutions for spin-up from rest in a cylinder, J. Fluid Mech. 127 (1983) 263-281. [12] M.A. Shadday, R.J. Ribando and J.J. Kauzlarich, Flow of an incompressible fluid in a partially filled, rapidly rotating cylinder with a differentially rotating endcap, J. Fluid Mech. 130 (1983) 203-218. [13] R.J. Ribando and M.A. Shadday, The Ekman matching condition in a partially filled, rapidly rotating cylinder, J. Comp. Phys. 53 (1984) 266-288. [14] R. Krauss, University of Virginia, Written communication (1986). [15] C.G. Popp, Countercurrent velocity measurements from the rotating reference frame, MS Thesis, University of Virginia, 1987. [16] H.E Greenspa- - r ~ Theory of Rotating Fluids (Cambridge Univ. Press, New York, NY, 1968). [17] G.E. Heuser, RJ. Ribando and H.G. Wood, III, A numerical simulation of inertial waves in a rotating fluid, Comput. Methods Appl. Mech. Engrg. 57 (1986) 207-222.
414
J.W. Frederick et al., A finite difference simulation of a centrifuge
[18] F.H. Harlow and J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (1965) 2182-2189. [19] C.W. Hirt, B.D. Nichols and N.C. Romero, SOLA- A numerical solution algorithm for transient fluid flows, LA-5852 (1975) Los Alamos National Laboratory. [20] R.D. Beardsley, K.D. Saunders, A.C. Warn-Varnas and J.M. Harding, An experimental and numerical study of the secular spin-up of a thermal stratified rotating fluid, 3. Fluid Mech. 93 (1979) 161-164. [21] H.P. Greenspan, Simulation of countercurrent flow in a gas centrifuge, UVA-ER-746-82U, University of Virginia, 1982. [22] P.3. Roache, Computational Fluid Dynamics (Hermosa, Albuquerque, 1976).