Prediction of Residence Time Distribution of Stirred Reactors

5686

Ind. Eng. Chem. Res. 2001, 40, 5686-5695

Prediction of Residence Time Distribution of Stirred Reactors Ashwin W. Patwardhan* Department of Chemical Technology, University of Mumbai, Matunga, Mumbai 400 019, India

In the present work, the ability of extending the measured velocity and turbulence field to predict the residence time distribution is explored. The model predictions have been compared with experimental measurements of residence time distribution. The predicted residence time distributions have been explained on the basis of the mixing behavior within the stirred reactor. The effect of operating parameters has been studied. 1. Introduction In chemical process industries, stirred vessels are frequently employed to carry out homogeneous as well as heterogeneous reactions, either in a batch or continuous mode. The duties that an impeller has to perform in a stirred reactor range from simple homogenization to the very complex like gas dispersion, solid suspension, heat transfer, or any combination of the above. It may also happen that the operation of the stirred reactor is controlled by the rate of chemical reaction (homogeneous or heterogeneous) taking place in it. When analyzing stirred reactors, it is commonly assumed that they behave like perfectly mixed tanks. Thus, the concentration everywhere within the tank equals the outlet concentration. This assumption implies that the rate of chemical reaction is constant throughout the tank. This allows the estimation of the reactor performance. This situation arises when the mixing time (blending time) is very small as compared to the mean residence time within the reactor. Under these conditions, the exit age distribution is represented as E(τ) ) 1/τ e-t/τ and the mean residence time is just the ratio of the reactor volume to the volumetric flow rate. When the mean residence time is comparable to or lower than the mixing time, the exit age distribution curve shows substantial deviation from the idealized back-mixed behavior.1 The concentration within the reactor is nonuniform. The extent of reaction of a packet of fluid element depends on its residence time in the reactor. The residence time distribution therefore becomes important. When there are two reactants, the situation is even more complex. Besides the residence time of the fluid elements, the scale and intensity of segregation within the reactor also play an important role. If the feed stream is well mixed, then the residence time distribution is the only important thing. Conventionally, the RTD can be described by various models,1 for example, CSTR’s with exchange volumes, CSTR with a dead zone, CSTR with a bypass, and so forth. However, these models contain many parameters, such as the mean residence time, the volume of the dead zone, the exchange flow rate, the bypass flow rate, etc., that can be varied to fit the experimental data. Many times it can also happen that the experimentally measured residence time distribution data can be fitted with more than one type of model. In that case, extracting reliable information becomes difficult. More* E-mail: [email protected]. Phone: 91-22-414 5616. Fax: 91-22-414 5614.

over the models mentioned above do not consider the flow field within the reactor which results in nonideal behavior in the first place. The mixing time in the reactor is a strong function of the impeller speed. Under turbulent flow conditions, the mixing time is inversely proportional to the impeller speed (Nθmix ) constant). Thus, the proportion of the mixing time and mean residence time can be changed simply by a change in the impeller speed. At high impeller speeds, the mixing time would be much smaller than the mean residence time and the reactor would behave essentially as back-mixed. If the impeller speed is reduced, the mixing time increases and may even become comparable or smaller than the mean residence time. This would lead to completely different exit age distribution. Thus, the exit age distribution would be highly dependent upon the relative magnitude of the mixing time and the mean residence time. The presence of dead zones, improper arrangement of inlet and outlet nozzles (short-circuiting), would complicate the situation further. This would have significant implications for the operation of a stirred reactor. Such experimental investigations have been reported in the published literature.2 In view of this, it would be worthwhile to see whether the exit age distribution could be predicted on the basis of flow field and mixing behavior occurring within the reactor. The exit age distribution once predicted can be used for a quick estimation of the level of conversion within the reactor. The effect of various parameters such as impeller type and impeller spacing on the exit age distribution can be used to give a quick idea about the expected benefits in terms of the reactor throughput, product quality, etc. This will be useful for the screening of different stirred reactor configurations, such as impeller location, location of inlet and outlet nozzles, etc., without actually constructing/modifying the fullscale reactor. Therefore, the aim of the present work was to explore the possibility of developing a model for predicting the residence time distribution based on the actual flow and turbulence fields present within the reactor. 2. Present Work The flow field within the reactor can be measured by a nonintrusive method such as laser doppler velocimetry (LDV). This can give information about the mean velocities and turbulence characteristics (turbulent kinetic energy, turbulent energy dissipation rates, etc.) at a point in the reactor. However, measuring the flow

10.1021/ie0103198 CCC: $20.00 © 2001 American Chemical Society Published on Web 10/31/2001

Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001 5687

Figure 1. Details of the experimental setup and the LDV measurement locations. Location A ) 18 mm below the impeller center plane. Location B ) 12 mm above the vessel bottom. Location C ) 18 mm above the impeller center plane. Location D ) 72 mm above the impeller center plane. Location E ) 36 mm below the impeller center plane. Location F ) 36 mm above the impeller center plane.

field throughout the stirred reactor using LDV is extremely time-consuming. Therefore, a computational fluid dynamics (CFD) model was used to extend the LDV measurements so as to characterize the flow field throughout the stirred reactor. The justification for such approach has been provided in the subsequent sections. The flow and turbulence field calculated by a combination of the LDV and CFD is used to predict the residence time distribution within the stirred reactor. 2.1. Measurement of Flow Field. The experimental setup for the LDV measurements consisted of 0.3 m i.d. transparent acrylic vessel equipped with a standard pitched blade downflow impeller (PBTD). The vessel was fitted with four baffles having width 1/10th of the vessel diameter (fully baffled condition). The impeller could be operated at different speeds using a control panel. The impeller diameter was kept equal to 1/3 of the tank diameter and was located at a distance equal to 1/3 from the tank bottom. The impeller consisted of six blades (W/D ) 0.3) at an angle of 45°. The LDV system consisted of laser beam generator, beam splitter, and focusing arrangement. The entire stirred reactor setup was mounted on a traversing system which enabled us to focus the laser beams at any location within the reactor. The Doppler shift produced was measured by a photo multiplier tube in the forward scatter mode. The data were acquired with the help of a PC. All of the LDV measurements were made midway between the two baffles. The experimental setup is shown in Figure 1. The mean velocities in all of the three directions (vr, vz, and vθ) and turbulent kinetic energy (k) were measured at six axial locations (about 15 points in the radial direction at each of the axial locations). 2.2. Prediction of Flow and Turbulence Fields. The conventional CFD simulations employ various impeller modeling approaches, namely, (i) black-box approach, (ii) multiple reference frames, (iii) inner-outer approach, and (iv) sliding mesh simulations. It has been reported that the sliding mesh3,4 simulations do not give good predictions in the near impeller region. On the other hand, the conventional black-box approach5-8 does not show good predictions in regions away from the impeller. Furthermore, with the above approaches, the predicted power number is not always in good agreement with the experimental measurements.9 For the geometry under investigation, it was observed that the

power number predicted with SM, MRF, and the snapshot approach is 1.3, 0.8, and 3.1, respectively,9 whereas the experimentally observed power number for a standard pitched blade turbine is 1.9. This means that these approaches do not give good predictions of turbulence quantities, and as a result, these models would not accurately predict the turbulent transport of the tracer which would be required in the prediction of the residence time distribution. Even the mean velocity fields predicted by these approaches may not be realistic. Initially, only a conventional black-box approach was used. It was observed that the velocity predictions (especially radial and tangential) were poor far away and far below the impeller. To overcome this problem and to force the CFD simulations to give good predictions (both near and away) from the impeller, the CFD simulations were carried out by specifying the experimental data (vr, vz, vθ, and k) at four axial locations (indicated as A, B, C, and D) in the vessel as shown in Figure 1. The velocity and turbulence fields were thus forced to agree with the experimental measurements, far above and far below the impeller. The CFD model was formulated to solve the Reynolds transport equations along with a standard “k - ” model for turbulence. These equations were written in a generalized manner as

1 ∂ ∂ 1 ∂ ∂φ 1 ∂ (vθφ) + (vzφ) ) rΓ + (rvrφ) + r ∂r r ∂θ ∂z r ∂r ∂r 1 ∂ Γ ∂φ ∂ ∂φ + Γ + Sφ (1) r ∂θ r ∂θ ∂z ∂z

(

( )

)

( )

The source terms for all of the flow variables (Sφ) has been given earlier.10 In the above equation, the variables were made dimensionless with the impeller tip velocity (Utip) and the radius of the vessel (R). Thus, the nondimensionalising parameter for mean velocities was Utip, for the coordinate system it was R, for the turbulent kinetic energy it was Utip2, for the turbulent energy dissipation rate it was Utip3/R, and for the eddy diffusivity it was RUtip. The dimensionless grid size in r, z, and θ directions was 0.033, 0.029, and 0.1308, respectively. The above equations were solved using the experimental data measured by LDV as boundary conditions. This was done in order to “bind” the CFD

5688

Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001

solution to experimental values even very far away from the impeller. In other words, CFD model was used only as a tool for extending the experimentally measured flow field. This has been done mainly because LDV measurements throughout the stirred vessel are extremely time-consuming. Such a combination of the LDV data, and the CFD simulations, enable use of realistic flow fields. 2.3. Prediction of the Residence Time Distribution. The flow and turbulence field predicted using the above methodology was used to solve the conservation equation for an inert tracer. The time averaged equation for the transport of the inert tracer is given below:

1 ∂ ∂ 1 ∂ ∂c 1 ∂ (rv′rc′) + + (rvrc) + (v c) + (vzc) + ∂t r ∂r r ∂θ θ ∂z r ∂r 1 ∂ ∂c ∂ 1 ∂ 1 ∂ (v′ c′) + (v′zc′) ) rDm + r ∂θ θ ∂z r ∂r ∂r r ∂θ Dm ∂c ∂c ∂ + D (2) r ∂θ ∂z m ∂z

(

(

)

) ( )

The correlation between fluctuating velocity and fluctuating concentration can be evaluated with the help of eddy diffusivity and the gradient of mean concentration (similar to the eddy viscosity). On further simplification, the resulting equation becomes

1 ∂ ∂ ∂c 1 ∂c + (rvrc) + (v c) + (vzc) ) ∂t r ∂t r ∂θ θ ∂z 1 ∂ Γm ∂c ∂ ∂c ∂c 1 ∂ rΓm + + Γ (3) r ∂r ∂r r ∂θ r ∂θ ∂z m ∂z

(

)

( ) ( )

The above equations are based on the following assumptions: (i) The correlation between concentration and velocity fluctuations can be represented in terms of eddy diffusivity for transport of the tracer and the mean concentration gradient. This is analogous to the approximation of the Reynolds stress terms. The eddy diffusivity for the transport of the tracer is taken equal to that for the transport of momentum divided by the Turbulent Schmidt number (vm ) v/SCT). The turbulent Schmidt number has been assumed to be unity. (ii) The above equation assumes that there are no other forces acting on the fluid elements of the tracer. This implies that the tracer is neutraly buoyant; therefore, there are no buoyancy or gravity effects. In other words the density and the viscosity of the tracer are identical to that of the liquid in the tank. (iii) It is assumed that the concentration gradients do not affect the velocity and turbulence field in any way. This enables computation of the velocity and turbulence fields independent of the concentration fields. (iv) It has been reported in the previous literature that the velocity and turbulence field is essentially at steady state outside the impeller zone. The velocity field within the impeller zone would be unsteady in nature because of successive passage of the impeller blades. Resolving the unsteady velocity fields in this region is not fairly important for the purpose of RTD estimations. This is because the time scale of the velocity fluctuations in the impeller region is quite small and the magnitudes of the velocities are quite large. As a result, the tracer spends very little time in the impeller zone. The overall residence time distribution would therefore be expected to be controlled by the transport of the tracer outside

the impeller zone, where the velocity field is at steady state. Thus, it would be reasonable to assume that the velocity field is at steady state. In fact, eqs 2 and 3 can also be modified to account for the reactions occurring within the reactor. The resulting equation would then represent the conservation equation for a reactive species. It is possible to solve the resulting equation and predict the performance of the reactor entirely through CFD. However, this is not attempted in the present work. When a homogeneous first-order reaction occurs within the reactor, the residence time distribution is sufficient to characterize the system. Further, when a reactive species is considered, one may need to consider the exothermcity/endothermicity of the reaction. This would necessitate a solution of the equation for the conservation of enthalpy. It may also be necessary to solve for the transport equations, turbulence model equations, enthalpy equations, and species conservation equations, simultaneously. When multiple reactions occur, it would be necessary to solve the conservation equation for each species. When the reactions have order higher than unity, or for reactions involving multiple reactants, an additional term appears in eq 2, namely, C′A2, or C′AC′B. The evaluation of these correlations depends on the intensity and scale of segregation and the extent of mixing. The above aspects make analysis from first principles alone difficult as well as time-consuming. Many times in industrial practice the understanding obtained about the flow field and residence time distribution is sufficient to shed light about the operation of the existing reactors and evaluate the effect of various parameters such as impeller type, impeller spacing, number of impellers, locations of the inlet and outlet nozzles, etc. In view of this, the aim of the present work has been limited to the prediction of the residence time distribution only. In the CFD simulations, the outlet nozzle was modeled as the sink of the tracer. The addition of the tracer pulse was simulated by specifying a very high value of concentration at the inlet of the reactor. The concentration at the outlet was calculated by the solution of the eq 3 by an implicit method using the control volume formulation. The variation of the outlet concentration versus time (either experimentally measured or predicted) was used to compute the mean residence time in the following manner:

∫0∞ tCout(t) dt τ) ∞ ∫0 Cout(t) dt

(4)

The outlet concentration as a function of time was used to compute the exit age distribution, E, curve in the following manner:

E(t) )

Cout(t)

∫0 Cout(t) dt ∞

(5)

The above exercise was repeated for different impeller speeds so as to vary the mixing characteristics (mixing time) in the reactor. The prediction of the residence time distribution was also studied at two different sets of the inlet and outlet nozzle geometries. In the first set, the inlet nozzle was kept near the shaft close to the top surface of the liquid, whereas the outlet nozzle was


Figure 2. Comparison of the CFD predictions and the LDV simulations. (s) CFD predictions; (b) LDV measurements; (A) dimensionless radial velocity at E; (B) dimensionless radial velocity at F; (C) dimensionless axial velocity at E; (D) dimensionless axial velocity at F.

provided near the vessel bottom close to the wall. In the second set, the inlet nozzle was kept near the shaft close to the impeller, whereas the outlet nozzle was provided near the top surface of the liquid close to the wall. In the solution of the above equation, the time step plays an important role especially when the concentration profiles are steep (especially at low values of time). Initially, some simulations were carried out so as to study the sensitivity of the time step and to obtain a time-step independent solution. At the start of the solution, the time step was of the order 0.0025 so as to resolve the steepness of the E-curve. When the E curve flattened out, the time step was gradually increased to about 0.25. 2.4. Measurement of Residence Time Distribution. The residence time distributions were also measured experimentally, in a 0.3 m diameter reactor having provisions to maintain constant level. This reactor was operated at a constant inlet flow rate and the speed of agitation was varied from 40 to 220 rpm in order to vary the mixing time within the reactor. The tracer pulse for the residence time distribution studies was a concentrated solution of sodium chloride (common salt). The outlet conductivity was measured right from the instant of giving the pulse till the conductivity fell to a very low value. Using, this conductivity, time data, the mean residence times (eq 4), and the E curve (eq 5) were determined at various impeller speeds. This served as a validation point for the CFD simulations. 3. Results and Discussion 3.1. CFD Model Validation. As a first step, the CFD model was validated by a comparison of the CFD predictions of vr, vz, vθ, and k with the LDV measurements at different locations. The locations at which the comparisons have been made are indicated by E and F in Figure 1. The comparison is shown in Figures 2 and 3. The prediction of radial velocity at location E (Figure 2A) is good only for r/R g 0.4. At lower r/R values, the predicted radial velocity is lower as compared to the experimental measurement. This could be due to the

Figure 3. Comparison of the CFD predictions and the LDV simulations. (s) CFD predictions; (b) LDV measurements; (A) dimensionless tangential velocity at E; (B) dimensionless tangential velocity at F; (C) dimensionless turbulent k.e. at E; (D) dimensionless turbulent k.e. at F.

disturbances induced by the hub, end of the shaft, etc. which may not have been captured completely. At location F (Figure 2B), the predicted radial velocity is lower in magnitude compared with the LDV measurements. The experimental data at this location itself shows considerable scatter, and even then, the CFD model predicts the correct trend. The predicted values of axial velocity at locations E and F (Figure 2 parts C and D, respectively) are in good agreement with the LDV measurements. The magnitude of the radial velocity at these locations is much smaller than the axial velocity. The convective transport of the tracer (which is of interest for the purpose of residence time distribution) will be predominantly in the axial direction. Hence, it is important that the axial velocity predictions are good. The predicted values of tangential velocity at locations E (Figure 3A) shows that the magnitude of the peak is lower by about 30% as compared with the experimental measurements. At location F (Figure 3B), the LDV measurements show scatter, which is partly responsible for the slight disagreement. The predicted values of turbulent kinetic energy at locations E (Figure 3C) agree fairly well with the LDV measurements. At location F (Figure 3D), the predicted values are much higher in magnitude for r/R e 0.4. This could again be due to the proximity of the shaft which can lead to dampening of the turbulence. However, for larger r/R values, the predictions are in excellent agreement with the LDV measurements. It is worth reemphasizing that none of the impeller modeling approaches9 gives good predictions for all of the flow and turbulence variables throughout the tank. Furthermore, the energy balance cannot be established properly with the unsteady approaches such as sliding mesh, multiple reference frames, inner-outer, etc. In view of this, it was necessary to “bind” the CFD solution to the experimental values at four locations (A, B, C, and D) within the reactor. Only when this was done, the overall energy balance was also satisfied; that is, the predicted power number 1.7 was found to be in good agreement with the experimentally observed power

5690


number of 1.9 for this geometry. Therefore, the flow field estimated in this manner with a combination of LDV and CFD can be used further with sufficient confidence for the prediction of the residence time distribution. In the future, it is likely that the advances in computing power together with better impeller and turbulence models; it would be possible to satisfactorily predict all the flow and turbulence characteristics accurately within the entire tank. 3.2. Mixing Characteristics. The previous literature on blending tends to indicate that the blending process is essentially controlled by the convective transport.11,12 The models proposed by these researchers employ impeller pumping number (which is a measure of the bulk flow, that is, convective transport) to correlate the data on blending time. The compartment models13 take into consideration the interstage flow rate as well as the circulation rate (convective transport) within the zone. The experimental investigations14,15 also indicate that the tracer tends to circulate (and simultaneously mix with the surrounding fluid) along the mean flow generated by the impeller. Patwardhan and Joshi16 have carried out detailed investigation of the blending process for a wide variety of axial flow impellers. They were able to correlate the mixing time with the secondary flow number of the impeller. The contribution of the convective transport and the dispersive transport to any transfer process (momentum or heat or mass) is characterized by means of a Peclet number (Pe). It characterizes the ratio of the convective transport to the dispersive transport and can be defined as

Pe )

UL D

Pe )

U ) xvr2 + vz2 + vθ2

(7)

The characteristic length scale was estimated from the turbulent kinetic energy dissipation rate as

(u′)3 k3/2 )A

(8)

The value of the dispersion coefficient was calculated as

D ) νm )

k2 ν ) cµ ScT

(9)

Using these definitions, the Peclet number can be calculated as 3/2

Pe )

xvr2 + vz2 + vθ2Ak k2 cµ

which gives

(6)

here U, L, and D represent the characteristic velocity scale, length scale, and dispersion coefficient, respectively. Therefore, it was thought desirable to compute the values of the Peclet number and its distribution throughout the stirred reactor using the flow field predicted through CFD. The characteristic velocity scale was calculated as

L)A

Figure 4. Predicted distribution of the Peclet numbers.

(10)

Typically the value of the constant A is taken as unity,

xvr2 + vz2 + vθ2 cµk1/2

(11)

This is analogous to the inverse of the turbulent intensity. As mentioned previously, the values of the mean velocities and turbulent kinetic energy were measured at four locations through LDV and supplied as boundary conditions for CFD simulations. The values measured by LDV would only be sufficient to calculate the Peclet number (eq 11) at these four locations. Thus, the Peclet number was calculated using the predicted local values of mean velocities and the turbulent kinetic energy throughout the reactor. The distribution of the Peclet numbers in the plane midway between the two baffles is shown in Figure 4. From the figure it can be seen that the Peclet number at all locations within the reactor is very high. Near the impeller, the values of the Peclet number are very high between 30 and 40. Close to the wall, the values of Pe are as high as 50. This is because of the decay of turbulence as the wall approaches. In a substantial portion near the vessel bottom, the values of Pe are between 10 and 20. Even near the top liquid level (substantial distance away from the impeller), the values of Pe are quite high. This indicates that the contribution of the convective transport is large in comparison to the turbulent (dispersive) transport. 3.3. Residence Time Distribution. The mean residence time predicted with eq 4 is compared with the experimentally observed mean residence time in Figure 5. The figure also shows the values of the mixing time as a function of impeller speed according to the value of Nθmix ) 40.1 determined experimentally.17


Figure 5. Comparison of the residence time and mixing time. (s) Predicted mean residence time; (b) experimentally measured mean residence time; (- - -) mixing time estimated from the experimental data.17

From the figure it can be seen that the mean residence time predicted by CFD simulations is in reasonably good agreement with the experimental measurements over a wide range of impeller speeds. At high impeller speeds, the mean residence time predicted by the CFD simulations is lower by 10% as compared with the experimental measurements. As the impeller speed increases to a very high value, the mean residence time approaches the value V/F. The experimentally measured mean residence time is calculated from the measured outlet concentration with the help of eq 4. The tail of the outlet concentration curve gets longer as the impeller speed reduces. Thus, the experimentally measured mean residence time distribution shows scatter about the ideal value (V/F). The predicted mean residence time is also calculated from eq 4, but with predicted outlet concentration. As the impeller speed reduces, the peak value in the E curve is shifted to right and the peak becomes broader. As a result of this, the numerator of eq 4 shows slightly higher value than the idealized behavior. Thus, it appears that the mean residence time increases marginally. However, this effect is very small. In fact, a severalfold reduction in the impeller speeds causes a very minor change in the mean residence time. The E curves obtained from the CFD simulations are shown in Figure 6 parts A and B for the first set of nozzle locations. From the figure, it can be seen that the magnitude of the peak in the E curve increases and the peak becomes sharper with an increase in the impeller speed. The peak also shifts closer to zero (leftwards) with an increase in the impeller speed. At large times, all of the E curves approach each other and decrease in an exponential manner (Figure 6B) and compare well with those estimated using mean residence time defined as V/F. Thus, the start of the peak is indicative of the time required for the tracer blob to move from the inlet to the outlet nozzle, whereas the magnitude of the peak is indicative of the extent of mixing (intensity of segregation) of the blob with the surrounding during its travel. The effect of nozzle arrangement on the E curves has been depicted in Figure 6 parts C and D. From the figure, it can be seen that, at 450 rpm, the E curve shifts to a lower value and is of smaller magnitude for the second set of nozzle arrangement as compared to the

Figure 6. E curve predicted by CFD simulations at various impellers speeds. (A) For the first 10 s (nozzle locations set 1). (B) For the first 30 s (nozzle locations set 1). Lines: (s) 450 rpm, (-‚-) 112 rpm, (‚‚‚) 80 rpm, (-‚‚‚-) 41 rpm. (b) Residence time distribution corresponding to ϑ ) V/F. (C and D) Comparison of the E curve for nozzle locations, sets 1 and 2 for the first 10 and 30 s, respectively. Lines: (s) 450 rpm nozzle locations set 1, (-‚-) 450 rpm nozzle locations set 2, (‚‚‚) 80 rpm nozzle locations set 1, (-‚‚‚-) 80 rpm nozzle locations set 2.

first set. On the other hand, the peak in the E curve is prolonged and is of lower magnitude at 80 rpm for the second set of nozzle configuration. The contour plots of the predicted concentration profiles at two impeller speeds (450 and 80 rpm for the first set of nozzle geometries) are compared in Figure 7. At time 0.005 s (Figure 7A), the tracer moves very small distance from the impeller. Because the impeller is PBTD type, it generates a downward flow near the inlet nozzle, as a result the tracer is pulled down toward the impeller. Because this graph corresponds to a time just after the pulse is given, the flow pattern in the liquid phase has not had sufficient time to mix the tracer. Hence, the concentration plots for both the speeds show the tracer as a blob close to the inlet nozzle. The shape and maximum value of the concentration plot is practically the same at both the speeds. At a slightly later time, 0.05 s (Figure 7B), the blob of the tracer still persists. Comparing the two figures, it can be seen that, at 450 rpm, the tracer blob has moved closer toward the impeller as compared to that for 80 rpm and also appears to be elongated slightly. The peak value of concentration within the blob is also lower for 450 rpm as compared to 80 rpm. This means that the tracer mixes more at 450 rpm than at 80 rpm, resulting in longer blob and reduced peak concentration. At a still later time, 0.4 s (Figure 7C), dramatic differences can be seen in the two plots. At 450 rpm, the blob of the tracer has moved substantially as compared to that at 80 rpm. Because the Peclet number in this region is very high, of the order 10-30 (Figure 4), the convective transport of the tracer dominates over the dispersive transport. As a result of this, the blob elongates but does not spread much radially. In contrast, the blob of the tracer is still having more or less the same shape and still lies well above the impeller at 80 rpm. The dominant mode of convective transport of the tracer in this region is the axial component of the velocity. That is why, the tracer elongates and does not spread much radially. Because the predicted axial velocities are in good agreement with the experimental measurements, the predicted mean residence time is

5692


Figure 7. Concentration plots at 450 and 80 rpm (nozzle locations shown by arrows). (A) 0.005 s, (B) 0.05 s, (C) 0.4 s, (D) 1.87 s. Horizontal axis represents dimensionless radial coordinate (0 at the center and 1 at the wall). Vertical axis represents dimensionless axial coordinate (0 at the base and 2 at the liquid surface).

also in good agreement with the experimental measurements. Thus, the poor prediction of radial velocity in some regions do not affect the results to a great extent. In other words, the mean residence time is relatively insensitive to the prediction of the radial velocity. At time 1.87 s, (Figure 7D), one can see a clear difference in the two plots. At 450 rpm, the tracer has spread completely in the vessel, in fact some of it has also gone out of the reactor. From these plots, it can be seen that the concentration at the outlet is much higher at 450 rpm as compared to 80 rpm. At this instant, the E curve corresponding to 450 rpm is close to its maximum value. Whereas, at 80 rpm, the tracer blob still persists and the E curve corresponding to this time is still close to zero.

The present model predictions also agree well with the earlier experimental results.2 It has been observed that,2 as the impeller speed decreases, the exit age distribution curve shows significant deviation from ideal back-mixed behavior. The magnitude of the peak in the E curve was also reported to increase with an increase in the impeller speed. This was attributed to the shortcircuiting from inlet to outlet. In fact, Figure 7 parts C and D show that a part of the tracer has gone out even before it mixed with the rest of the reactor at 450 rpm. The effect of nozzle location is compared in Figure 8 parts A-D for an impeller speed of 450 rpm. From the figure, it can be seen that, at small times (0.005 s; Figure 8A), the tracer blob is close to the inlet nozzle. At 0.05 s (Figure 8B), the tracer blob elongates quicker


Figure 8. Concentration plots at 450 rpm (nozzle locations shown by arrows). (A) 0.005 s, (B) 0.05 s, (C) 0.4 s, (D) 1.87 s. Horizontal axis represents dimensionless radial coordinate (0 at the center and 1 at the wall). Vertical axis represents dimensionless axial coordinate (0 at the base and 2 at the liquid surface).

when the nozzle is close to the impeller. At 0.4 s (Figure 8C), the tracer blob spreads and mixes more when the inlet nozzle is closer to the impeller. As a result, the peak in the E curve shifts to lower values but is lower in magnitude (Figure 6 parts C and D) when the inlet nozzle is closer to the impellers. The effect of nozzle locations at 80 rpm is compared in Figure 9 parts A-D. From the figures, it can be seen that, at small times (0.005 s; Figure 9A), the tracer blob is close to the inlet nozzle. At 0.05 and 0.4 s (Figure 9 parts B and C), the tracer blob elongates quicker when the nozzle is close to the impeller. Beyond this point, the motion of the tracer for the second set of nozzle arrangement is slower because it is traveling further away from the impeller as compared to the first set of

nozzle arrangement. In fact, at 1.87 s (Figure 9D), the tracer has reached quite close to the outlet location for the first set of nozzle geometry, whereas it is still at a large distance from the outlet nozzle for the second set of nozzle geometry. This leads to a shift in the E curve to larger times and overall a reduction in the magnitude of the peak for the second set of nozzle arrangement. This is clearly seen in Figure 6 parts C and D. 4. Conclusions In the present work, the possibility of predicting the residence time distribution based on the velocity and turbulence field has been explored. This was done by characterizing the flow produced by a standard pitched

5694


Figure 9. Concentration plots at 80 rpm (nozzle locations shown by arrows). (A) 0.005 s, (B) 0.05 s, (C) 0.4 s, (D) 1.87 s. Horizontal axis represents dimensionless radial coordinate (0 at the center and 1 at the wall). Vertical axis represents dimensionless axial coordinate (0 at the base and 2 at the liquid surface).

blade turbine, with the help of LDV measurements and CFD simulations. The velocity field obtained in this manner has been used to solve the conservation equation for a tracer and to predict the residence time distribution. The observed mean residence time is in good agreement with the experimental measurements. The observed residence time distribution and the exit age distribution have been explained on the basis of the concentration profiles within the reactor. Nomenclature A ) constant in eq 8 c ) tracer concentration at any location, kmol/m3 CFD ) computational fluid dynamics

CSTR ) continuous stirred tank reactor C′A2 ) time averaged square fluctuating concentration of A, kmol2/m6 C′AC′B ) time averaged product of fluctuating concentrations of A and B, kmol2/m6 Cout(t) ) tracer concentration at the outlet at any time t, kmol/m3 Cµ ) parameter used in k - model, eq 9 D ) dispersion coefficient, m2/s Dm ) molecular diffusivity of the tracer m2/s E(t) ) value of the E curve at time t F ) flow rate of inlet to the reactor, m3/s k ) turbulent kinetic energy, m2/s2 L ) characteristic length scale, m

Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001 5695 LDV ) laser doppler velocimetry MRF ) multiple reference frame PBTD ) pitched blade downflow turbine Pe ) Peclet number R ) vessel radius, m RTD ) residence time distribution r ) radial coordinate, m SM ) sliding mesh ScT ) turbulent Schmidt number Sφ ) source term for a generalized flow variable φ t ) time, s U ) characteristic velocity, m/s u′ ) characteristic velocity scale used in eq 8, m/s Utip ) impeller tip velocity, m/s V ) volume of the stirred reactor, m3 vr, vz, and vθ ) mean velocity in the radial, axial, and tangential direction, respectively, m/s vr′, vz′, and vθ′ ) fluctuating velocity in the radial, axial, and tangential direction, respectively, m/s z ) axial coordinate, m Greek Letters Γ ) eddy diffusivity plus diffusivity for momentum, m2/s Γm ) eddy diffusivity plus molecular diffusivity for the tracer, m2/s ) energy dissipated per unit mass, W/kg θ ) tangential coordinate ν ) eddy diffusivity for momentum, m2/s νm ) eddy diffusivity for transport of tracer, m2/s τ ) mean residence time calculated from eq 4, s φ ) generalized flow variable

Literature Cited (1) Fogler, H. S. Elements of Chemical Reaction Engineering. Prentice Hall: Englewood Cliffs, NJ, 1986. (2) Chollete, A.; Cloutier, L. Mixing Efficiency Determinations for Continuous Flow Systems. Can. J. Chem. Eng. 1959, 37, 107. (3) Bakker, A.; Myers, K. J.; Ward, R. W.; Lee, C. K. The Laminar and Turbulent Flow Pattern of A Pitched Blade Turbine. Chem. Eng. Res. Des. 1996, 74, 485-491. (4) Bakker, A.; Laroche, R. D.; Wang, M. H.; Calabris, R. V. Sliding Mesh Simulation of Laminar Flow in Stirred Reactors. Chem. Eng. Res. Des. 1997, 75, 42-44.

(5) Ranade, V. V.; Joshi, J. B.; Marathe, A. G. Flow generated by pitched blade turbines II: Simulation using k - model. Chem. Eng. Commun. 1989, 81, 225-248. (6) Kresta, S. M.; Wood, P. E.; Prediction of the threedimensional turbulent flow in stirred tanks. AIChE J. 1991, 3, 448-460. (7) Armenante, P. M.; Chou, C. Velocity profiles in a baffled vessel with single or double pitched blade turbines. AIChE J. 1996, 42, 42-54. (8) Harris, C. K.; Roekaerts, D.; Rosendal, F. J. J. Computational Fluid Dynamics for Chemical Reactor Engineering. Chem. Eng. Sci. 1996, 51, 1569-1594. (9) Nere, N. K.; Patwardhan, A. W.; Ranade, V. V.; Joshi, J. B. CFD modeling of stirred tank reactors. Chem. Eng. Sci. 2001 Submitted for publication. (10) Sahu, A. K.; Kumar, P.; Joshi, J. B. Simulation of Flow in Stirred Vessels with Axial Flow Impeller: Zonal Modeling and Optimization of Parameters. Ind. Eng. Chem. Res. 1998, 37, 21162130. (11) McManamey, W. J. A. Circulation Model for Batch Mixing in Agitated, Baffled Vessels. Chem. Eng. Res. Des. 1980, 58, 271275. (12) Rzyski, E. Liquid Homogenization in Agitated Tanks. Chem. Eng. Technol. 1993, 16, 229. (13) Vasconcelos, J. M. T.; Barata, J. M.; Alves, S. S. Transitional Mixing in Multiple Turbine Agitated Tanks. Chem. Eng. J. 1996, 63, 53-58. (14) Williams, R. A.; Mann, R.; Dickin, F. J.; Ilyas, O. M.; Ying, P.; Edwards, R. B.; Rushton, A. Application of Electrical Impedance Tomography to Mixing in Stirred Vessels. AIChE Symp. Ser. 1993, 89, 8-15; No. 293. (15) Mann, R.; Pillai, S. K.; El-Hamouz, A. M.; Ying, P.; Togatorop, A.; Edwards R. B. Computational Fluid Mixing for Stirred Vessels: Progress from Seeing to Believing. Chem. Eng. J. 1995, 59, 39-50. (16) Patwardhan, A. W.; Joshi, J. B.; Relation Between Flow Pattern and Blending in Stirred Tanks. Ind. Eng. Chem. Res. 1999, 38, 3131-3143. (17) Rewatkar, V. B.; Joshi, J. B. Effect of Impeller Design on Liquid Phase Mixing in Mechanically Agitated Reactors. Chem. Eng. Commun. 1991, 91, 322-353.

Received for review April 11, 2001 Revised manuscript received July 23, 2001 Accepted August 15, 2001 IE0103198

Prediction of Residence Time Distribution of Stirred Reactors

Recommend Documents