Computational Flow Modeling and Visualization in the Annular Region

Ind. Eng. Chem. Res. 2007, 46, 8343-8354

8343

Computational Flow Modeling and Visualization in the Annular Region of Annular Centrifugal Extractor Sandesh S. Deshmukh, Sreepriya Vedantam, and Jyeshtharaj B. Joshi* Institute of Chemical Technology, UniVersity of Mumbai, Matunga, Mumbai-400 019, India

Sudhir B. Koganti Indira Gandhi Centre for Atomic Research, Kalpakkam, TN-603 102, India

Flow between two concentric cylinders, with either or both of them rotating, has potential advantages over the conventional process equipment. This flow, which is also termed Taylor-Couette flow and exhibits a variety of flow regimes, has been studied using computational fluid dynamics (CFD) simulations. The onset of centrifugal instability, the various cell patterns, and the velocity profiles have been predicted and compared with the available experimental data. To extend the base of experimental information, new measurements have been made using laser Doppler velocimetry (LDV) and particle image velocimetry (PIV). For the entire range of experimental data, the Reynolds stress model (RSM) was determined to exhibit good predictive ability for the mean components of velocity and the turbulent kinetic energy. The wavelengths of the vortices in various regimes have also been determined using CFD and were observed to be in good conformity with all the available experimental results. Furthermore, very good agreement between the predicted energy dissipation rate with the energy input rate was observed over a wide range of speeds and annular gaps. 1. Introduction A system that involves two concentric cylinders, with either or both rotating (which is also known as a Taylor-Couette contactor), has attracted attention as efficient process equipment. The duct through which the fluid flows is the annulus between the two cylindrical surfaces, which is small, compared to the radii of the coaxial cylinders. The flow is circumferential due to the shearing action of the rotating cylinders. Furthermore, when such equipment is operated in a continuous mode, an axial component is introduced and a pressure differential is setup between the ends of the annulus. In addition to the axial and circumferential flow, three-dimensional secondary flows are generated when the rotational speed exceeds a certain value. Typical vortex patterns, in the absence of an axial flow, are shown in Figure 1. Such a flow pattern is of great interest, because it provides high values of the heat- and mass-transfer coefficients. Some important industrial operations that can be performed advantageously in Taylor-vortex equipment include liquid-liquid extraction,1,2 emulsion polymerization,3 synthesis of silica particles,4 heterogeneous catalytic reactions,5 and also bio-reactions.6 In all these cases, for the rational design, it is desirable to develop design procedures for the prediction of pressure drop, flow regimes, axial mixing, and the heat- and mass-transfer coefficients. The reliability of such procedures is dependent on the extent of knowledge of the flow pattern in Taylor-vortex contactors. The work described here forms part of an investigation designed to develop an improved understanding of the Taylorvortex single-phase and multiphase flows, particularly in an annular centrifugal extractor (ACE), which operates as a mixer and settler with a very low residence time. The importance of an accurate description of the flow pattern has been highlighted through the use of computational fluid dynamics (CFD). In the * To whom all correspondence should be addressed: Tel.: 00-9122-2414 0865. Fax: 00-91-22-2414 5614. E-mail address: [email protected].

present work, single-phase simulations within the annulus of a Taylor-vortex flow system have been conducted. Both horizontal and vertical geometries have been considered. 2. Previous Work The Taylor-Couette flow has been extensively investigated since the end of the nineteenth century, and the state-of-the-art reviews have been presented by Koschmieder7 and Vedantam and Joshi.8 A brief account has been presented below. The first set of experiments in Taylor-Couette flow systems were reported by Couette9 and Mallock.10 During their drag measurement experiments, they reported instability at a certain rotational speed of the cylinders. However, this instability was first analyzed by Rayleigh11 for inviscid flows. For this case, he concluded that the flow is unstable if the two cylinders rotate in opposite directions. If they rotate in the same direction, then the motion is stable if Ωoro2 is greater than Ωiri2. For real fluids, Taylor12 theoretically analyzed the stability problem using the Bessel Fourier series approach and showed good agreement with his experimental measurements. However, the analysis was restricted to narrow gaps (ri/ro ) 0.94). The case of narrow gaps has also been addressed by Jeffreys,13 and the principal objective was to understand the analogy between the TaylorCouette flow and the Rayliegh-Benard convection. He used a linear theory of stability, and the exercise of analogy led to the introduction of new dimensionless number, which was later named the Taylor number (Ta) by Chandrasekhar.14 The value of Ta at which the first transition occurs from Couette flow (CF) to Taylor-vortex flow (TVF) was designated as the critical Taylor number (TaCr). Chandrasekhar14 obtained a more-general solution for instability for the entire range of rotational ratio of the cylinders for both narrow as well as wide gaps. The validity of Chandrasekhar’s analysis was confirmed by Donnelly15 by measuring torque for the two cases of radius ratio, η ) 0.95 and 0.5 with an aspect ratio of 50 and 5, respectively.

10.1021/ie061587e CCC: $37.00 © 2007 American Chemical Society Published on Web 06/29/2007

8344

Ind. Eng. Chem. Res., Vol. 46, No. 25, 2007

Figure 1. System of concentric cylinders and various vortex patterns in the absence of axial flow.

Figure 3. (A) Onset of centrifugal instability, shown through a comparison of the CFD predictions with the experimental data of Taylor,12 with H ) 0.084 m and d ) 2.35 mm. (B) Effect of viscosity on the centrifugal instability (legend for panel B: ([) experimental (from Taylor12), (1) µ ) 10-3 kg/(m s), (2) µ ) 10-4 kg/(m s), (3) µ ) 10-5 kg/(m s), (4) µ ) 10-6 kg/(m s), (5) µ ) 10-7 kg/(m s), and (6) Rayleigh line).

Figure 2. Details of the geometry used in the experimental measurements and simulations. (All dimensions in figure are given in millimeters.)

The case of wide gaps was also addressed by Roberts16 and he found the critical Taylor number TaCr to vary with the gap, according to the following relationship:

TaCr ) 1589.2η-1.0964

(1)

For the above case, the outer cylinder was stationary. In eq 1, when applied to the Taylor experiment (η ) 0.94), the value of TaCr becomes 1701, which is in fairly close agreement with the value reported by Taylor (1703).12 The effect of net axial flow on stability was analyzed by Chandrasekhar,17 and he showed that the introduction of axial flow delays the occurrence of instability, according to

TaCr ) 1708 + 27.15Rez2

(2)

As Ta continues to increase beyond TaCr, for the case of a stationary outer cylinder, the flow structure changes to wavy vortex flow (WVF) (TaCr < Ta < 100TaCr), and then to chaotic vortex flow (CVF) (100TaCr < Ta < 1000TaCr), and finally to turbulent Taylor vortex flow (TTVF) (Ta > 1000TaCr). These flow patterns are known as supercritical Taylor vortices, having a characteristic feature of asymmetric disturbances. Therefore, the linear theory of instability is not appropriate, and, hence, Stuart,18 Davey,19 and DiPrima and Eagles20 have used the nonlinear theory of instability to understand the supercritical vortices. Coles21 visualized wavy vortices by suspending aluminum particles in silicon oil and has reported interesting observations. The wavy nature of the Taylor vortex was determined to be dependent on the way in which the rotational speed of the cylinders is varied. The numbers of vortices were also determined to be dependent on the methodology used to increase or decrease speed. Burkhalter and Koschmieder22 performed experiments with a stationary outer cylinder, with radius ratios of η ) 0.727 and 0.505 and aspect ratios of either Γ ) 51 or 28. They provided

Ind. Eng. Chem. Res., Vol. 46, No. 25, 2007 8345

rotating plates at the gas-liquid free surface. The important objective of investigation was to bring out the role of bounded ends where the zero slip was assumed. The range of Taylor number was that of WVF. The authors have reported interesting observations: (i) The size of the two end vortices increases with an increase in the value of Ta, (ii) The size of the bulk vortices increases with an increase in the overall length of annulus; however, it reaches an asymptotic value when H/d is >51, and (iii) The number of vortices reduce with an increase in the gap. Koschmieder23 measured the wavelength of turbulent Taylor vortices up to 40 000TaCr, using visualization experiments. Although he did not theoretically analyze the flow systems, his studies reported the effect of start-up procedure on the wavelength of vortices. Andereck et al.24 used flow visualization (laser light scattering) to construct flow maps that show the various regimes from Couette flow to turbulent Taylor vortex flow. The time dependences of the flows have been studied by measuring the intensity of laser light scattered by polymeric flakes. Andereck et al.24 used a batch mode of operation. Later, Lueptow et al.25 presented flow maps in the presence of axial flow. They also measured the wavelengths of vortices for different flow regimes. For highly turbulent regimes, Parker and Merati26 used laser Doppler velocimetry (LDV) to measure three components of mean velocity and turbulent intensity at various circumferential planes. For the aspect ratio of 4 and 20, they studied the end effects on the vortices. To predict the magnitude of vortex velocities and behavior beyond the critical point, Baier27 used a CFD approach. Haut et al.28 also performed CFD simulations in the annular region in horizontal co-rotating cylinders, for the wavy vortex flow and the turbulent Taylor vortex flow. They used the standard k- model to account for the turbulence parameters. The velocity profiles obtained from the simulations for both regimes were in good agreement with their experimental results. The maximum angular velocity they used was 5 rad/s (Ta ) 5.35 × 105). They have investigated instantaneous and time-averaged velocity fields experimentally, using particle image velocimetry (PIV). From the foregoing discussion, it is clear that substantial experimental and theoretical work has been published for the regime transitions. However, only few studies have been reported on the quantitative flow pattern. Therefore, it was thought desirable to use CFD simulation to predict flow patterns over a wide range of Reynolds and Taylor numbers. It was also thought desirable to undertake measurement using LDV and PIV. In the present work, an attempt has been made to present an extensive comparison between all the experimental measurements reported in the published literature (including the present measurements) and the CFD predictions. These include the three velocity components, the turbulent kinetic energy, and the dissipation rate. Extensive comparison has been presented for the first transition Taylor number (TaCr) over a wide range of rotor diameter, rotation ratio, annular gap, and liquid viscosity values. Furthermore, a complete energy balance has also been established. It may be emphasized that all the previous investigations were focused only on the annular region. For the real annular centrifugal extractor, flow must also be analyzed in the region below the rotating cylinder. This aspect has been addressed in this investigation by undertaking LDA and PIV measurements, as well as CFD simulations.

3. Mathematical Modeling 3.1. Model Formulation. When the problem is axisymmetric, with respect to geometry and flow conditions, including swirl or rotation, the flow can be modeled in two dimensions (i.e., solve the axisymmetric problem) and include the prediction of the circumferential (or swirl) velocity. This analysis can be verified by comparing the tangential (swirling) velocity profiles with the measured profiles. Hence, in the present case, twodimensional (2D) simulations have been performed for an axisymmetric swirling flow between two concentric cylinders. The Navier-Stokes equations for an incompressible, constant viscosity liquid can be written in cylindrical coordinates as follows. Continuity:

∂F ∂(Fuz) ∂(Fur) Fur + + + )0 ∂t ∂z ∂r r

(3)

Axial Component:

∂(Fuz) 1 ∂(rFuzuz) 1 ∂(rFuruz) + + ) ∂t r ∂z r ∂r ∂uz 2 ∂P 1 ∂ + b) rµ 2 - (∇‚V ∂r r ∂r ∂z 3

[(

)]

(4)

Radial Component:

∂(Fur) 1 ∂(rFuzur) 1 ∂(rFurur) ∂P + + )+ ∂t r ∂z r ∂r ∂r ∂ur 2(∇‚V ∂ur ∂uz b) 1 ∂ 1 ∂ rµ + rµ 2 + r ∂z ∂z ∂r r ∂r ∂r 3 2 ur 2µ(∇‚V b) Fuθ 2µ 2 + + (5) 3r r r

[(

)]

[(

)]

wherein

∇‚V b)

∂uz ∂ur ur + + ∂z ∂r r

(6)

Tangential Component:

∂(Fuθ) 1 ∂(rFuzuθ) 1 ∂(rFuruθ) + + ) ∂t r ∂z r ∂r ∂uθ 1 ∂ ∂ uθ 1 ∂ rµ + 2 r3µ r ∂z ∂z ∂r ∂r r r

( )

[ ( )]

-

Furuθ (7) r

In swirling flows, the conservation of angular momentum has a tendency to create a free vortex flow in which the angular velocity increases as the radius decreases. In an ideal vortex flow, the centrifugal forces created by a circumferential motion are in equilibrium with the radial pressure gradient, i.e.,

∂P FΩ2 ) ∂r r

(8)

As the distribution of angular momentum in nonideal vortex evolves, the form of this radial pressure gradient also changes, driving radial and axial flows, in response to highly nonuniform pressures that result therein. 3.2. Turbulence Modeling. For turbulence modeling, two models have been used, viz, the standard k- model and Reynolds stress model (RSM). A few of the details are given below.

8346


and the rate of dissipation of turbulent kinetic energy, , is given as

(( ) ) (( ) ) ( ) ( )

µt ∂ ∂ ∂ 1 ∂ 1 ∂ r µ+ + (F) + (Fur) + (Fuz) ) ∂t r ∂r ∂z r ∂r σ ∂r µt ∂ 2 ∂ µ+ + C1 Gk - C2F (10) ∂z σ ∂r k k In these equations, Gk represents the generation of turbulence kinetic energy because of the mean velocity gradients, and the turbulent viscosity (µt) is computed as

()

µt ) Fcµ

Figure 4. Parity plot for comparison of the power dissipation predictions between correlation by Arafat et al.31 and CFD for experimental geometry using water.

k2

(11)

and the model constants are C1 ) 1.44, C2 ) 1.92, and Cµ ) 0.09. 3.2.2. The Reynolds Stress Model (RSM). In this model, individual Reynolds stresses (u′iu′j) are computed via a differential transport equation. Thus, the RSM model solves six Reynolds stress transport equations. Along with these, an equation for dissipation rate is also solved. The exact form of Reynolds stress transport equations is derived by taking moments of an exact momentum equation. This is a process wherein the exact momentum equations are multiplied by a fluctuating property, with the product then being Reynolds averaged. Using Cartesian coordinates for conciseness, the transport equations for Reynolds stresses, Fu′iu′j, are given by

∂[ Fu′iu′ju′k + P(δkju′i + δiku′j) ] ∂(Fu′iu′j) ∂(Fuku′iu′j) )+ + ∂t ∂xk ∂xk

[ ] (

)

∂(u′iu′j) ∂uj ∂ui ∂ µ - F u′iu′k + u′ju′k + ∂xk ∂xk ∂xk ∂xk P

Figure 5. Onset of centrifugal instability: Comparison of CFD predictions with the experimental data of Chandrasekhar,14 with H ) 0.084 m and d ) 2.35 mm. Table 1. Comparison of Energy Dissipation a Dissipated Energy (kW/tonne) angular velocity (rad/s)

Arafat et al.31 correlation

RSM

k-

62.8 125.6

0.56 3.97

0.6 3.85

0.62 3.90

Determined using the expression 2πF ∫rrio ∫z0 r dr dz.

The governing equations for the geometry under consideration are given by the following models. 3.2.1. The Standard k- Model. The turbulence kinetic energy, k, is given as

∂ 1 ∂ ∂ (Fk) + (rFkur) + (Fkuz) ) ∂t r ∂r ∂z µt ∂k µt ∂k ∂ 1 ∂ r µ+ + µ+ + Gk - F (9) r ∂r σk ∂r ∂z σk ∂z

( ( ) ) (( ) )

∂u′i ∂u′j ∂u′i ∂u′j - 2µ (12) ∂xj ∂xi ∂xk ∂xk

The second term on the right-hand side represents turbulent diffusive transport. It has been simplified in FLUENT software32 to use a scalar turbulent diffusivity as ∂/∂xk [µt/σk ∂(u′iu′j) /∂xk] (from Launder et al.29) with σk ) 0.82 and µt is computed using eq 11. The turbulence kinetic energy was obtained by taking the trace of Reynolds stress tensor,

1 k ) (u′iu′i) 2

CFD

a

( )

(13)

The fourth term in eq 12 represents the pressure strain term, modeled using the linear pressure strain model explained in FLUENT software,32 whereas the last term (dissipation rate of turbulent stresses) is modeled as (2/3)δijF. Here, is computed with the same transport equation (eq 10) as that used by the k- model. It may be noted, at this stage, that, in the FLUENT software, the equations are written and solved in the Cartesian system and an appropriate transformation of coordinates is included during solution of the equations and presentation of the results. 3.3. Boundary Conditions. The angular velocities of the walls have been specified (in the range of 0-125.6 rad/s). The rotation ratio of the cylinders was varied from -2.24 to 1. The negative value indicates the counter rotation of the cylinders, and the zero value indicates the stationary outer cylinder. Across


Figure 6. Velocity contours (υ ) 0.552, η ) 0.826): (a) axial velocity (expressed in units of m/s), (b) azimuthal velocity (expressed in units of m/s) (c) radial velocity (expressed in units of m/s), and (d) stream function (expressed in units of kg/s).

are obtained for the sudden start of the experiment; henceforth, the velocities have been specified for a rotating boundary and transient simulations were performed. In the case of a gradual increase of speed, simulations were started at a lower rotational speed initially. After attaining steady state, the speed was increased again. This stepwise procedure was repeated until the final angular velocity was obtained. 3.4. Boundary Conditions for Reynolds Stresses. The standard wall functions are based on those of Launder and Spalding.30 At the walls, the near-wall values of the Reynolds stresses and are computed from wall functions. The explicit wall boundary conditions are applied for Reynolds stresses, using log-law and the assumption of equilibrium, thus disregarding convection and diffusion in the transport equations for stresses. 4. Numerical Framework: Method of Solution Figure 7. Vortex spacing as a function of the rotation ratios: comparison of CFD predictions with the experimental data of Taylor,12 with H ) 0.084 m and d ) 0.0235 m.

the flow regime, no variation of velocities along the tangential direction was considered (∂/∂θ ) 0). In the case of stability studies, axial boundary conditions were given as being transitionally periodic. It has been specified that no-slip conditions hold at all the walls. The different lengths of the annuli covered are 0.01, 0.305, 0.9, and 1 m. The range of gap width used in these simulations has been 2-15 mm. The Ta range used for pure circumferential flow with no axial flow, in simulations, was 0-2 × 109. In the case of the present LDV and PIV measurements for turbulent Taylor-vortex flow investigation, the two angular velocities (62.8 and 125.6 rad/s) have been chosen. The results

With the finite-volume formulation, all the simulations were performed using axisymmetric 2D grids. The commercial FLUENT software (Versions 4.5 and 6) have been used in all the studies. While a uniform scheme that consisted of 200 × 100 grids in axial and radial directions were used in the case of the predictions of previous experimental data, nonuniform grids are used for the experimental geometry. A segregated implicit solver method with implicit linearization was used for the momentum equations. The momentum equations have been discretized with the quick upwind scheme, and for the pressure velocity coupling, the PISO scheme has been used. For the pressure equation, the PRESTO scheme was used. Eccentricity (ratio of the offset distance of the cylinder axis to the average gap width) was assumed to be zero. Data were collected at specified points to track the development of the flow and confirm that the asymptotic solution was attained.

8348


Figure 8. Formation of the secondary set of Taylor vortices at a high rotation ratio, with H ) 0.084 m, d ) 0.0235 m, and υ ) -2.38.

Figure 9. Flow patterns obtained by PIV and CFD at (A) 62.8 rad/s and (B) 125.6 rad/s.

5. Experimental Section The details of geometry investigated are shown in Figure 2. The upper end was kept stationary, while the bottom end was modified by keeping the rotor 8 mm above the bottom of the stationary outer cylinder. The shaded portion shows the measurement window for flow visualization and anemometry experiments. The rotor and outer vessel was composed of acrylic, which has the refractive index of 1.49. The difficulties arising in the measurements near the wall that were due to curvature effects have been resolved using a refractive index matching technique that was previously used by Parker and Merati.26 For this issue, a solution of sodium iodide in water was used. The viscosity and the density of the solution were measured and found to be 1.5233 cS and 1.837 g/cm3 at 24 °C,

respectively, and the refractive index of solution was measured to be 1.48 ( 0.05. The rotor speeds at which measurements have been conducted were 62.8 and 125.6 rad/s with a sudden start at the respective speeds. 5.1. PIV Measurements. The PIV measurements were performed using a TSI PIV system. The laser source adopted was a pulsed Nd:YAG laser that had a pulse duration of 6 ns and was synchronized with the camera using a synchronizer. The optics included a combination of cylindrical and spherical lenses, attached in view of creating a thin laser sheet ∼1 mm thick. The images were captured using a high-resolution 4M CCD camera (with a freqeuncy of 15 Hz), placed perpendicular to the laser sheet, with 2048 × 2048 pixels. The liquid was seeded with silver-coated hollow glass particles with a mean


Figure 10. Comparison of the axial and radial velocity predictions (A) 62.8 rad/s and (B) 125.6 rad/s with experimental measurements at the axial location z/H ) 0.60 ((s) RSM model, (- - -) k- model, and (0) experimental results).

diameter equal to 20 µm. The interrogation area was set at 64 × 64 pixels with a 50% overlap resulting in ∼1500 vectors for the images. The time difference between the two laser pulses was optimized based on Nyquest criterion. The personal computer (PC) functioned also as a central data acquisition, as well as a processing unit. Image processing of vectors extracted from PIV images is done by passing the signal-to-noise ratio using threshold value of 4 so that most of the bad vectors are eliminated. The post-processing is then done where vectors are compared with neighboring vectors. The vectors that vary by more than the validation tolerance from neighborhood average are removed. After the vectors have been validated, the missing points are filled in and the properties of the flow are computed. 5.2. LDV Measurements. A three-beam LDV system was used to measure the two velocity components (axial-radial and axial-tangential) simultaneously. The LDV that was set up was comprised of Dantec 55 X modular series optics, along with electronic instrumentation and a PC. A 5-W Ar-ion laser (Spectra Physics) was used as a source. To identify the flow direction, a frequency shift of 40 MHz is given to one laser beam. All the optics are from Dantec Dynamics, and the system includes the cover and retarder plates (to adjust the beam polarization), beam splitters, Bragg cell, beam expander, and 600-mm front lens for focusing. The measurement volume is formed by three mutually perpendicular laser beams, viz, blue

(588 nm), green (614 nm), and cyan (combination of blue and green) beams. The scattered laser light from the measurement volume is captured by photomultiplier tubes in a forward scatter mode. Because the tap water naturally contains seeding for the scattering of light, the flow was not artificially seeded. Data validation and signal processing (on-line fast Fourier transform) were performed by Burst Spectrum Analyzer (62N40 BSA F60 processor) from Dantec Dynamics, which consists of two velocity channels. The average data rate was ∼700 Hz, with a validity >90%. The entire operation of data acquisition, including the high voltage to photomultiplier tubes and the record length selection for the burst detection were controlled by a personal computer using BSA Flow software version 4.0. 6. Results and Discussions Initially, three-dimensional (3D) as well as 2D axisymmetric simulations were performed to understand the onset of instability. Figure 3A shows that the 3D and 2D simulations give practically the same results. Furthermore, to provide a good basis for the validity of 2D-CFD, it was thought desirable to establish a complete energy balance. The power consumption in an annular centrifugal extractor (ACE) has been analyzed by Arafat et al.31 For this ACE case, where the inner cylinder rotates and the outer one is stationary,

8350


Figure 12. Comparison of CFD predictions of axial velocity with the experimental data of Haut et al.:28 η ) 0.826, Ωi ) 2.15 rad/s, Γ ) 10.5 ((s) CFD and (9) experimental results).

Figure 11. Comparison of the tangential velocity predictions:(A) 62.8 rad/s and (B) 125.6 rad/s with experimental measurements at an axial location of z/H ) 0.60 ((s) RSM model, (- - -) k- model, and (0) experimental results).

the power consumption was correlated by the following equation:

Pw )

(µd)

0.0261Hri3.75(jΩ)2.75F0.75 m

0.25

(14)

where

j ) 0.0554(log Re) + 1.368 Re )

2jFΩrid µ

and

3 × 103 < Re < 1 × 106 For the purpose of prediction, the power consumption for the 2D case was computed from the moment about a specified center. The FLUENT software32 calculates this moment vector by summing the product of the force vectors for each face with the moment vector, i.e., summing the forces on each face about the moment center. Figure 4 shows an excellent agreement between the CFD predictions and the correlation of Arafat et

al.31 Furthermore, the energy input rate (by the rotation of any one or both the cylinders) must equal the energy dissipation rate (by both the viscous and turbulent modes of dissipation). The comparison is shown in Table 1. CFD shows a very good match at lower angular velocity (62.8 rad/s). When the speed is doubled (125.6 rad/s), the predicted values for the RSM and k- models increase by a factor of 6.4 and 6.3, respectively. Furthermore, the experimental value can be increase correspondingly by a factor of 7.1. This can be considered to be fairly close agreement. In view of the successful energy balance, as well as the success of the prediction of the onset of instability in two dimensions, all the subsequent simulations have been performed using the 2D framework. 6.1. Centrifugal Instability. The model validation has been presented in terms of comparison of the onset of instability with the experimental data of Taylor,12 as well as the predictions of Rayleigh11 and Chandrasekhar.14 Figure 3A shows very good agreement between the CFD predictions and the experimental data of Taylor.12 Furthermore, it can be seen that, as the rotation rate of outer cylinder increases, viscous forces become less important, as compared to the centrifugal forces, and the Rayleigh criterion is approached as an asymptote. To understand the effect of viscosity over a wide range of rotational speeds, simulations were performed over a wide range of viscosities ranging from 1 × 10-6 to 2 × 10-3 Pa s. The gradual shift of the stability curve from the Taylor fluid (µ ) 0.001 Pa s) toward the Rayleigh line has been captured in Figure 3B. Because the Taylor methodology developed a criterion for the case of narrow gaps, the cases of wide gaps have been analyzed by Chandrasekhar.14 Figure 5 shows the comparison of the predictions of Chandrasekhar14 and those of CFD and agreement can be seen to be reasonably good. With the onset of instability, as the azimuthal flow is superimposed over the circular Couette-flow, doughnut-shaped vortices form along the length of the annulus. A typical flow pattern arising due to the onset of instability has also been captured from CFD simulations. The stream function contours depict the Taylor vortex pattern as alternate cells, where consecutive cells move in the same direction at their meeting


Figure 13. Comparison of CFD predictions of axial velocity with the experimental data of Parker and Merati26 at z/d ) 0.66, η ) 0.672, Ωi ) 125.6 rad/s, Γ ) ((s) CFD and (9) experimental results).

Figure 15. (A) Velocity time series data obtained by LDV at 62.8 rad/s [location: (r - ri)/d ) 0.6 and z/H ) 0.55] (B) Ratio of eddy viscosity to molecular viscosity at various axial locations: ([) z/H ) 0.22, (4) z/H ) 0.55, and (0) z/H ) 0.88).

Figure 14. Comparison of CFD predictions of radial velocity with the experimental data of Parker and Merati26 at z/d ) 0.66, η ) 0.672, Ωi ) 125.6 rad/s, and Γ ) 20 ((s) CFD and (9) experimental results).

point. For the case of a rotation ratio of 0.552, the radial velocity contours, axial velocity contours, and the swirl velocity contours are shown in Figure 6. This pattern formation is in qualitative agreement with the predictions reported in the literature.27 For quantitative comparison, it was thought desirable to compare the vortex sizes over a wide range of rotational speeds. The results in the form of vortex spacing (the distance between the vortex centers) have been shown in Figure 7. The CFD results are in very good agreement with the experimental data of Taylor.12 In the case of counter-rotating cylinders, as the Ta was increased from TaCr, at a very high rotational speed, a secondary set of vortices started to occur, thus spanning in the gap width between the cylinders. This observation is once again in agreement with those observed in the past.12 This phenomenon has been termed as a secondary bifurcation in the mathematical analysis of Taylor.12 The formation of this secondary set of vortices is shown in Figure 8.

6.2. Flow Patterns. After an occurrence of Taylor instability (Ta ) TaCr), counter-rotating toroidal vortices are observed in Taylor-Couette flow. As the rotational speed increased, the number of vortices decreased. The experimental geometry is similar to that of the ACE, and the present simulations have been performed for the single-phase case. For a batch-mode operation and for the case of 62.8 rad/s, the number of vortices was determined to be 8, which decreased to 6 when the speed was doubled to 125.6 rad/s. Furthermore, elongated vortices at the ends (top and bottom) and compressed vortices at the center were observed, which are in agreement with the observations of Parker and Merati.26 However, the vortex at the bottom is observed to be stretched more than the vortex that is at the top. Although Taylor vortices are observed in the annulus, the flow below the rotating cylinder is a tangential rotating flow. Figure 9 shows velocity vectors obtained in CFD simulations using the RSM model, which was also confirmed by PIV measurements. If end effects are neglected, with η ) 0.75 and a 6-mm annular gap for an aspect ratio of Γ ) 12, the number of vortices would be 12. Flow patterns for the experiments available in the literature show the importance of end effects. For the same operating conditions, as well as ri and ro (Ta ) 2 × 109), as used by Parker and Merati,26 if the top and bottom effects are neglected, the wavelength was determined to be 2d, so that 4 and 20 vortices can be observed with Γ ) 4 and 20, respectively. When

8352


Figure 16. Comparison of turbulent properties obtained by CFD and experimental measurements at (A) 62.8 rad/s and (B) 125.6 rad/s at axial location z/H ) 0.60 ((s) RSM model, (- - -) k- model, and (9) experimental results).

the end effects were included in the simulation, only 2 vortices were observed for the case of Γ ) 4, which is in complete agreement with their experimental findings. For the case of Γ ) 20, CFD simulations showed 16 vortices, versus 18 as reported by Parker and Merati.26 In the context of the aforementioned role of end effects, the comparison with the measurements of Haut et al.28 is noteworthy. For their case of Γ ) 10.5 and an angular velocity of 2.15 rad/s (Ta ) 1.02 × 105), the CFD simulations were in exact agreement with the experimental values of 10 vortices, as observed by Haut et al.28 For a detailed understanding of end effects, the geometry has been further investigated for various angular velocities, up to 70 rad/s (Ta ) 1.05 × 108). Note that, under the same conditions, for an infinitely long column, each vortex is of equal size (9 mm). When the speed was gradually increased, the vortices at the ends started to elongate and middle vortices started to compress. The number of vortices remained 10 until the angular velocity increased to 65 rad/s (Ta ) 9 × 107); however, at 70 rad/s, the number of vortices suddenly decreased to 6. The wavelength measurements by Burkhalter and Koschmieder22 were limited, up to a Taylor number of TaCr ) 80. The end effects were then limited to the vortices subsequent to end plates (elongation of end vortices), whereas all the other vortices were approximately the same size. However, as the Taylor number was gradually increased further,

it was observed that the effect was transferred to the middle vortices. The dominant alternate elongation-compression is observed and this continued until two pairs of vortices disappeared. In the case of a sudden start, at these speeds, the number remained 10 until 30 rad/s (Ta ) 1.93 × 107) and 8 vortices were observed for all other speeds of >30 rad/s until a speed of 70 rad/s was attained. The dominant elongation-compression was observed only at a speed of 30 rad/s, after which the transition occurred. Furthermore, the numbers of vortices was determined to be strongly dependent on the initial conditions. 6.3. Mean Velocities. Mean axial and radial velocities have been compared for both angular velocities of 62.8 and 125.6 rad/s. The mean velocities have been made nondimensional, with respect to the surface velocity of the rotating cylinder. Figure 10 shows the comparison between the experimental measurements using PIV and the CFD predictions. Both k- and RSM models were used. Figure 10A shows that the predictions of both models agree very well with the experimental data when the angular velocity is 62.8 rad/s. However, at 125.6 rad/s, the RSM model can be observed to perform very well (see Figure 10B). The k- model does not provide even qualitative predictions. Furthermore, the validity for axisymmetric swirl modeling was checked by comparing the tangential velocities measured by LDV. In Taylor-Couette-type flow, the tangential velocity is dominant (10 times of axial and radial velocities),


which is maximum at the inner rotating wall and zero at the stationary outer wall. The comparison has been shown in Figure 11. Note that, for the tangential velocity, predictions of both the turbulence models agree very well with the experimental data. Such an agreement validates the assumption of axisymmetric flow made in this work. It was also thought desirable to confirm the validity of the present CFD simulations by comparing them with the experimental results published in the literature. For this purpose, simulations were performed for both horizontal and vertical contactors. In the case of horizontal system, simulations have been performed for an annular gap of 8.5 mm and compared with the experimental data of Haut et al.28 (Ta ) 1.02 × 105) and shown in Figure 12, where very good agreement is observed when the RSM model was used. The RSM model hence has also been used for comparison with the experimental results of Parker and Merati,26 who used vertical geometry. They have measured the velocity components using LDV. The axial and radial velocity components at z/d of 0.66 are shown in Figures 13 and 14, respectively. Again, the agreement can be seen to be good for both the velocity component. The choice of turbulent model in the CFD simulations is observed to be very important to predict the flow accurately. The advantage of the k- model lies in the requirement of lower computational power than that requires for the RSM model. In the case of Taylor-Couette flow, predictions of the k- model were determined to be fairly good for an angular velocity of 62.8 rad/s (see Figure 10A), even though the flow is in the turbulent Taylor vortex regime. For weakly turbulent and wavy vortex regimes, Haut et al.28 could predict the mean axial velocities accurately using the k- model. However, at high velocities (125.6 rad/s), the end effects dominate and influence the flow pattern. The assumption of isotropy in the k- model fails when the flow is highly turbulent. In this case, the RSM model was determined to be more suitable for the turbulent modeling. 6.4. Turbulent Characteristics. Note that the operational range covered in this work falls within the turbulent Taylor vortex regime (Ta > 1000TaCr). For the selected geometry (ri ) 19.5 mm, d ) 6 mm, µ ) 2.75 mPa s, F ) 1837 kg/m3), the TaCr value works out to be 2178.16 At the two angular velocities of 62.8 and 125.6 rad/s, the values of Taylor number are 7.44 × 106 and 3 × 107, respectively, which are at least 3,500 times greater than TaCr. This indicates that the flow is turbulent. The CFD simulation shows that most of the turbulent energy dissipation occurs near the wall boundaries. Even then the flow away from the wall was also found to be turbulent. For instance, (a) the velocity time series data obtained by LDV measurements at 62.8 rad/s [at (r - ri)/d ) 0.6 and z/H ) 0.55] is shown in Figure 15A. It can be seen that the rms velocity is 80.5 mm/s and the turbulent intensity is 32.8%. (b) The ratio of eddy viscosity (µt) to molecular viscosity (µ) has been shown in Figure 15B. It can be observed that the minimum value of µt is at least 2 times higher than µ. The aforementioned two situations clearly show that the flow is turbulent. It should be interesting to note that, for the case of turbulent pipe flow with Re ) 22 000, the maximum turbulent kinetic energy was determined to be 0.014 m2/s2 near the wall, whereas, at the center, the values were as low as 0.003 m2/s2.33 In the case of the ACE, the corresponding values were determined to be 0.02 and 0.004 m2/s2 at an angular velocity of 62.8 rad/s, which are comparable to those of turbulent pipe flow. Note that,

in both cases, the dissipation is highest near the wall and decreases with an increase in the distance from the wall. The prediction of turbulent kinetic energy and its dissipation rates are compared with the measured values by LDV in Figure 16. Very good agreement has been observed between the measured values and those predicted by CFD at lower speed (62.8 rad/s) by both the k- and RSM models. Simulations conducted using RSM shows better predictions. The CFD show high levels of k and in the near-wall region. Furthermore, Figure 16 shows that the comparison for is very good at both velocities, whereas k is overpredicted at the higher speed. It has already been mentioned previously that the volume integral of (or the total predicted energy dissipation) equals the energy input rates, according to the correlation given by Arafat et al.31 7. Conclusions (1) Computational fluid dynamics (CFD) simulations have been performed in the annulus for a wide range of Taylor numbers (Ta) and Reynolds numbers (Re). For turbulence, the Reynolds stress model (RSM) model was determined to be suitable over the entire range of Ta and Re values, whereas the k- model was observed to have restricted applicability for the flow simulation in annular centrifugal extractors (ACEs). (2) The CFD model has been validated by establishing a complete energy balance. Thus, the volume integral of (predicted by CFD) was determined to be equal to the energy input rate. Excellent agreement was also found when the predictions were based on the torque experienced by the inner rotating cylinder. Furthermore, the energy balance has been established over a wide range of Ta values and annular gaps. (3) The transition from purely circumferential flow to Taylor vortex flow has been successfully predicted by CFD simulations. These predictions were determined to be in excellent agreement with (i) the Rayleigh criterion for inviscid flow, (ii) Taylor’s criterion and experimental results for narrow gaps,and (iii) Chandrasekhar’s criterion for wide gaps. Furthermore, the CFD simulations correctly predict the number of cells over a wide range of Ta values. (4) The CFD was also determined to simulate successfully the existence of two rows of Taylor vortices after a certain critical Taylor number (TaCr), in the case of counter-rotating cylinders. (5) Laser Doppler velocimetry (LDV) and particle image velocimetry (PIV) measurements have been performed for the actual ACE geometry. (6) In case of the ACE, elongated vortices were observed at both ends, which is in agreement with previous observations; however, the vortex at the lower end shows more elongation, because of the change in geometry, and the flow is almost tangential in the region below the rotating cylinder. (7) An excellent agreement has been observed between the predicted and the experimental profiles of mean velocity components, turbulent kinetic energy, and the energy dissipation rate. Acknowledgment This work has been part of the project program supported by Board of Research in Nuclear Sciences (BRNS), Sanction No. 2002/34/7-BRNS/140. The authors also acknowledge Dr. L. M. Gantayet for the discussions provided during this study. Nomenclature C1, C2, Cµ ) turbulent parameters in k- model d ) gap width (m)

8354


Gk ) generation of turbulent kinetic energy due to mean velocity gradients (kg m-1 s3) H ) annulus height/length (m) j ) factor defined in eq 16 k ) turbulent kinetic energy (m2/s2) m ) mass of the liquid (kg) P ) pressure (Pa) Pw ) power dissipated (kW/tonne) Re ) Reynolds number (defined in eq 14) Rez ) axial Reynolds number; Rez ) ri2Ω/ν ri ) inner cylinder radius (m) ro ) outer cylinder radius (m) Ta ) Taylor number; Ta ) 2η2/[(1 - η2)d4(Ωi/ν)2] TaCr ) critical Taylor number u ) velocity (m/s) u′ ) fluctuating velocity (m/s) uI ) instantaneous velocity (m/s) Us ) velocity at the surface of rotating cylinder (m/s) Greek Symbols δ ) Kronecker delta Ω ) angular velocity (rad/s) F ) density (kg/m3) ν ) kinematic viscosity (m2/s) µ ) molecular viscosity (Pa s) µt ) eddy or turbulent viscosity (kg m/s) η ) radius ratio; η ) ri/ro υ ) rotation ratio; υ ) Ωo/Ωi ) turbulent energy dissipation rate (m2/s3) λ ) wavelength (m) Γ ) aspect ratio; Γ ) H/d θ ) azimuthal coordinate σ ) quantity which determines stability σk ) turbulent parameter Subscripts i, j, k ) Cartesian coordinate components o ) Outer cylinder r ) radial direction z ) axial direction θ ) azimuthal direction Literature Cited (1) Davis, M. W.; Weber, E. J. Liquid-liquid extraction between rotating concentric cylinders. Ind. Eng. Chem. Res. 1960, 52, 929. (2) Bernstein, G. J.; Grodsvenor, D. E.; Lenc, J. F.; Levitz, N. M. Development and performance of a high-speed annular centrifugal contactor. Report No. ANL-7968, Argonne National Laboratory, Argonne, IL, 1973. (3) Imamura, T.; Saito, K.; Ishikura, S. A new approach to continuous emulsion polymerization. Polym. Int. 1993, 30, 203. (4) Ogihara, T.; Matsuda, G.; Yanagawa, T.; Ogata, N.; Fujita, Nomura M. Continuous synthesis of monodispersed silica particles using CouetteTaylor vortex flow. J. Ceram. Soc. Jpn. Int. 1995, 103, 151. (5) Sczechowski, J. G.; Koval, C. A.; Noble, R. D. A Taylor vortex reactor for heterogeneous photocatalysis. Chem. Eng. Sci. 1995, 50, 3163. (6) Tsao, Y. M. D.; Boyd, E.; Spaulding, G. Fluid dynamics within a rotating bioreactor in space and earth environments. J. Spacecr. Rockets. 1994, 31, 937.

(7) Koschmieder, E. L. Benard Cells and Taylor Vortices; Cambridge University Press: New York, 1993. (8) Vedantam, S.; Joshi, J. B. Design of annular centrifugal extractors: A review. Chem. Eng. Res. Des. 2006, 84 (A5), 522. (9) Couette, M. Etudes sur le frottement des liquides. Ann. Chim. Phys. 1890, 6, 433. (10) Mallock, A. Experiments on fluid viscosity. Philos. Trans. R. Soc. London, A 1896, 187A, 41. (11) Rayleigh, L. On the dynamics of revolving fluids. Proc. R. Soc. London, A 1916, 93, 148. (12) Taylor, G. I. Stability of a viscous liquid contained between two rotating cylinders. Phil. Trans. Soc. London, A 1923, 223, 289. (13) Jeffreys, H. Some Cases of instability in fluid motion. Proc. R. Soc. London, A 1928, 118, 195. (14) Chandrasekhar, S. The stability of viscous flow between rotating cylinders. Mathematica 1954, 1, 5. (15) Donnelly, R. J. Experiments on the stability of viscous flow between rotating cylinders. I. Torque measurements. Proc. R. Soc. London, A 1958, 246, 312. (16) Roberts, P. H. The solution of characteristic value problem. Proc. R. Soc. London, A 1965, 283, 550. (17) Chandrasekhar, S. The stability of spiral flow between rotating cylinders. Proc. R. Soc. London, A 1962, 265, 188. (18) Stuart, J. T. On the nonlinear mechanics of hydrodynamics stability. J. Fluid Mech. 1958, 4, 1. (19) Davey, A. The growth of Taylor vortices in flow between rotating cylinders. J. Fluid Mech. 1962, 14, 336. (20) DiPrima, R. C.; Eagles, P. M. Amplification rates and torques for Taylor vortex flows between rotating cylinders. Phys. Fluids 1977, 20, 171. (21) Coles, D. Transition in circular Couette Flow. J. Fluid Mech. 1965, 21, 385. (22) Burkhalter, J. E.; Koschmieder, E. L. Steady supercritical Taylor vortex flow. J. Fluid Mech. 1973, 58, 547. (23) Koschmieder, E. L. Turbulent Taylor vortex flow, J. Fluid Mech. 1979, 93, 515. (24) Andereck, C. D.; Liu, S. S.; Swinney, H. L. Flow regimes in a circular Couette system with independently rotating cylinders. J. Fluid Mech. 1986, 164, 155. (25) Lueptow, R. M.; Docter, A.; Min, K. Stability of axial flow in an annulus with a rotating inner cylinder. Phys. Fluids A 1992, 4, 2446. (26) Parker, J.; Merati, P. An investigation of turbulent Taylor-Couette flow using Laser Doppler Velocimetry in a refractive index matched facility. Trans. ASME 1996, 118, 810. (27) Baier, G. Liquid-liquid extraction based on a new flow pattern: two-fluid Taylor-Couette flow, Ph.D. Thesis, University of Wisconsin, Madison, WI, 2000. (28) Haut, B.; Ben, A. H.; Coulon, L.; Jacquet, A.; Halloin, V. Hydrodynamics and mass transfer in a Couette-Taylor bioreactor for the culture of animal cells. Chem. Eng. Sci. 2003, 58, 774. (29) Launder, B. E.; Reece, G. J.; Rodi, W. Progress in the development of a Reynolds-stress turbulence closure. J. Fluid Mech. 1975, 68 (3), 537. (30) Launder, B. E.; Spalding, D. B. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 1974, 3, 269. (31) Arafat, H. A.; Hash, M. C.; Hebden, A. S.; Leonard, R. A. Characterization and recovery of solvent entrained during the use of centrifugal contactor, Report No. ANL-02/08, Argonne National Laboratory, Argonne, IL, 2001. (32) FLUENT 6.0 User Guide, FLUENT, Lebanon, NH, 2002. (33) Schildknecht, M., Miller, J. A., Meier, G. E. A The influence of suction on the structure of turbulence in fully turbulent pipe flow. J. Fluid Mech. 1979, 90, 67.

ReceiVed for reView December 11, 2006 ReVised manuscript receiVed March 18, 2007 Accepted May 11, 2007 IE061587E

Computational Flow Modeling and Visualization in the Annular Region

Recommend Documents