Gas–liquid flow and bubble size distribution in stirred tanks

Chemical Engineering Science 63 (2008) 2107 – 2118 www.elsevier.com/locate/ces

Gas–liquid flow and bubble size distribution in stirred tanks G. Montante ∗ , D. Horn1 , A. Paglianti Dipartimento di Ingegneria Chimica, Mineraria e delle Tecnologie Ambientali, Università di Bologna, Viale Risorgimento 2-40136 Bologna, Italy Received 19 October 2007; received in revised form 20 December 2007; accepted 8 January 2008 Available online 15 January 2008

Abstract This work is aimed at investigating the turbulent two-phase flow and the bubble size distribution (BSD) in aerated stirred tanks by experiments and computational fluid dynamics (CFD) modelling. The experimental data were collected using a two-phase particle image velocimetry technique and a digital image processing method based on a threshold criterion. With the former technique, the liquid and the gas phase ensemble-averaged mean and r.m.s. velocities are measured simultaneously, while with the latter the dimensions of the bubbles dispersed inside the liquid are evaluated. On the modelling side, a CFD approach, based on the solution of Reynolds average Navier–Stokes equations in an Eulerian framework for both phases, is adopted. As for the bubble dimensions modelling, besides the mono-dispersed assumption, a population balance method, named MUSIG, with bubble break-up and coalescence models is considered. The BSD and the axial and radial velocity of the gas and the liquid phase are presented and discussed. The outcomes of the computational work are evaluated on the basis of the experimental results. 䉷 2008 Elsevier Ltd. All rights reserved. Keywords: Gas–liquid stirred tank; Two-phase PIV; Bubble size distribution; CFD

1. Introduction Gas–liquid stirred tanks are widely adopted in several industrial processes and many efforts have been devoted to their investigation by experimental methods and modelling approaches. Nevertheless, the fluid dynamics complexity of the turbulent, three-dimensional, two-phase flow complemented with coalescence and break-up of bubbles has not allowed yet to devise fully comprehensive modelling procedures and well-established design rules. So far, the local hydrodynamic features of stirred gas dispersions have been investigated in many experimental works, which mainly concerned the measurement of the liquid phase flow field in the presence of bubbles by LDV or PIV techniques (e.g. Rousar and van den Akker, 1994; Morud and Hjertager, 1996; Deen and Hjertager, 2002; Khopkar et al., 2003; Aubin et al., 2004), while very limited data are available on bubble ∗ Corresponding author. Tel.: +39 0512090406; fax: +39 0516347788.

E-mail address: [email protected] (G. Montante). 1 Present address: Prague Institute of Chemical Technology, Department

of Chemical Engineering, 16628 Praha, Czech Republic. 0009-2509/$ - see front matter 䉷 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2008.01.005

velocities and on the turbulent characteristics of the two-phase flow (Deen et al., 2002; Montante et al., 2007). Also, bubble size in stirred dispersions has been measured by different techniques (e.g. Barigou and Greaves, 1992; Takahashi and Nienow, 1993; Machon et al., 1997; Shäfer et al., 2000; Alves et al., 2002; Kamiwano et al., 2003; Hu et al., 2005; Laakkonen et al., 2005) and the available bubble size results have allowed to improve the understanding of important aspects of the stirred tanks applications, such as mass transfer; though, a spread of literature data exists, mainly due to limitations of the different measurement methods (Laakkonen et al., 2005). Shortage of reliable experimental bubble size distribution (BSD) for the assessment of computational fluid dynamics (CFD) simulations was recently pointed out by van den Akker (2006). In the applications of CFD simulations to gas–liquid stirred tanks, the BSD has been modelled only in few cases (Bakker and van den Akker, 1994; Venneker et al., 2002; Kerdouss et al., 2006; Laakkonen et al., 2007), rather the mono-disperse assumption has usually been preferred (e.g. Deen et al., 2002; Lane et al., 2004; Khopkar and Ranade, 2006; Scargiali et al., 2007; Montante et al., 2007). Apparently, although useful

2108

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

results have been obtained so far, quantitative evaluations of the prediction accuracy for liquid and gas turbulent flow field and for BSD are still scant in the literature. Moreover, model parameter fitting was often adopted to predict BSD accurately. Overall, many efforts have been devoted to the development of CFD models for multiphase flow, but often though surprisingly, the more advanced the computational techniques, the lesser their experimental validation. A comprehensive review of the CFD methods for turbulent mixing, including the analysis and the discussion of the literature results on gas–liquid systems, can be found elsewhere (van den Akker, 2006). In this work, we aim at gaining insight in both the two-phase turbulent flow field and the BSD in gas–liquid stirred tanks by up-to-date experimental techniques applied at the same time. The study is limited to volumetric gas holdup lower than 1%. A strict assessment of bubble and liquid mean velocities and bubble diameters predicted by adopting an industrially feasible CFD method based on Reynolds averaged Navier–Stokes (RANS) equations is presented. 2. The gas–liquid stirred tank The work was carried out in a gas–liquid stirred tank equal to that adopted in a previous investigation by Montante et al. (2007). It consisted of a fully baffled flat-bottomed and opentopped cylindrical vessel (T = 23.6 cm diameter) stirred by a standard Rushton turbine (D = T /3 diameter) placed at the distance C = T /2 from the vessel base. The experimental vessel was placed inside a square tank filled with the working liquid for minimising the optical errors. The liquid phase was demineralised water and its level inside the ungassed vessel was maintained at the height H equal to T . Air bubbles were sparged at a distance equal to T /4 from the vessel bottom through a porous membrane fixed on the top of a tube of 3.3 mm in diameter placed on the vessel axis. For the bubble size measurements, the impeller rotational speed, N, and the gas flow rate, QG , were varied from 3.3 to 7.5 s−1 and from 0.02 to 0.07 vvm, respectively, while a single condition (N = 7.5 s−1 and QG = 0.02 vvm) was selected for the PIV measurements as well as for the CFD simulations. 3. The experimental techniques The bubble size measurements were performed by capturing and processing several pictures of the gas–liquid stirred tank, which was illuminated adopting a diffused/incoherent light emitted by five tubular neon lamps at 50 kHz. The number, dimensions and position of the lamps were chosen for obtaining an almost uniform illumination of the selected vessel portion. The digital camera adopted for the image acquisitions (1344 × 1024 pixels CCD) was placed at the vessel side opposite with respect to the lamps, thus enhancing the contrast between the image background and bubble shadow. The exposure time was optimised to avoid the distortion effects due to the bubbles’ motion. Only the bubbles contained in the thin slab identified by a focal depth of a few millimetres were measured; in this way a reliable calibration factor could be adopted,

while renouncing quantitative information on the overall bubble holdup. The calibration factor between the CCD pixels and the real images was determined by means of an object of known size placed in the focused region, that was approximately midway between two consecutive baffles. The minimum detectable bubble size was equal to 0.3 mm. The above-mentioned upper limit of the impeller speed was due to the necessity to avoid significant air entrainment from the vessel top, while the constrain on the maximum gas flow rate was forced by the image quality requirements for bubble boundary detection. After its acquisition, each picture was transferred from the camera to the computer for processing by the shadow processing software (Dantec Dynamics). The adopted procedure started with the analysis of image background intensity and its equalisation for maximising the image contrast and reducing possible errors due to imperfect uniformity of illumination. Then, the bubbles’ edges were detected using the Canny method, based on the identification of local maxima of the image gradients—the gradient being calculated using the derivative of a Gaussian filter. Successively, the edges weaker than a fixed sensitivity threshold were discarded and those included in the output were further analysed for retaining only the edges with closed contours. In this way, the out-of-plane bubbles were filtered out since the intensity of their edges was much weaker than that of the bubbles contained in the focal volume. The threshold value was chosen by a trial and error procedure and the value (typically of the order of 0.25) that maximised the detected bubbles’ number was adopted. Finally, the area of each detected bubble was calculated summing up the total number of pixels contained inside the edge. The bubble equivalent diameter, db , was evaluated as the diameter of a spherical bubble having the same area of the measured object. The liquid and gas turbulent ensemble-averaged flow fields were measured using a Dantec Dynamics PIV system based on the adoption a single laser sheet source (pulsed Nd:YAG laser,  = 532 nm light wavelength, f = 15 Hz, frequency) and two cameras provided with appropriate light filters for separately capturing the light scattered by the liquid seeding particles (fluorescent Rhodamine-B) with one camera (HiSense MK II, 1344 × 1024 pixels CCD) and that of the bubbles with the other (PCO, 1280 × 1024 pixels CCD). The time interval between two laser pulses was equal to 150 s and that from a pulses pair to another was equal to a few seconds. Further details on the experimental system and procedures can be found elsewhere (Montante et al., 2007). At least 2000 image pairs were required to obtain sample-independent velocity fluctuation measurements. The interrogation area was set at 64 × 64 pixel and the cross-correlation of the image pairs was performed on a rectangular grid with 50% overlap between adjacent cells. For improving image resolution, the view area of the cameras was always limited to half of the vessel in both the bubble size and the PIV measurements. In the following, the origin of the coordinate system is located at the centre of the vessel bottom; z and r are the axial and radial coordinates, respectively. The mean axial velocity, U , is positive if directed upwards and the mean radial velocity, V , is positive if directed towards the vessel wall.

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

4. The CFD model equations The simulations of the gas–liquid stirred vessel were based on the Eulerian–Eulerian formulation of the continuity and momentum equations for each of the phases, which are considered as interpenetrating continua interacting by the interphase transfer terms. Because of the turbulent fluid flow regime, the RANS equations for each phase have been selected and their mathematical closure was obtained with the simplest extension of the standard k. turbulence model to two phase flows, named “homogeneous”, based on the assumption that the gas and liquid phases share the same values of k and . Both the mono-disperse assumption and the multiple size group model (MUSIG) proposed by Lo (AEA Technology, 1999) were considered for the bubble size. In the former case, our modelling approach coincide with that usually adopted for solid–liquid and gas–liquid stirred vessels, whose model equations can be found elsewhere (e.g. Montante and Magelli, 2005; Khopkar and Ranade, 2006; Scargiali et al., 2007). With the MUSIG model, the BSD is divided into a suitable number of size classes for which the continuity equations, including source and sink terms due to coalescence and break-up, are solved. For each bubble size class, j , the following continuity equation was considered: j (  ) + ∇ · k,j k Uk = Sj , (1) jt k,j k where  is the phase volume fraction,  is the density, U is the velocity vector and S is the source term associated to bubble coalescence and break-up. The subscripts k and j identify the phase (liquid or gas) and the bubble class, respectively. For k = liquid the source term vanishes and Eq. (1) coincides with the liquid continuity equation. It is worth observing that the continuity equation for the j bubble class divided by the phase density coincides with a population balance equation, that is written in terms of bubble number density. For bubble break-up and coalescence, the models of Luo and Svendsen (1996) and Prince and Blanch (1990), respectively, were selected. The former model was developed for predicting fluid–particle break-up rates in a turbulent dispersion and it is based on the theories of isotropic turbulence and probability, assuming binary break-up. The break-up rate, gB , of bubbles of volume vj and size dbj into bubbles of volume vj · fbv and vj · (1 − fbv ) is given by gB (vj : vj fbv ) 

 = 0.923 2 db ⎛ j × exp ⎝−

1/3 

1

(1 + )2

min

11/3

2/3

12[fbv + (1 − fbv )2/3 − 1] 5/3 l 2/3 dbj 11/3

turbulent energy dissipation rate,  is the surface tension of water, which was always assumed equal to 0.076 N m−1 ,  is a model constant. It is worth noting that no unknown or adjustable parameters are introduced in the model. The coalescence model considers the collision rate of two bubbles as the combination of turbulence, buoyancy and laminar shear effects, while the collision efficiency relates the time required for coalescence and the bubbles’ contact time. The birth rate of bubbles of class j due to coalescence of bubbles of the classes i and k is given by 1 Bc = Aik [u2tk + u2ti ]1/2 e(−tik / ik ) nk ni , 2 j

j

(3)

i=1 k=1

where ni and nk are the concentration of bubbles of diameter dbi and dbk , respectively, Aik is the collision cross-sectional area of the bubble, ut is the turbulent velocity, tik is the time required for coalescence, ik is the contact time of two bubbles. Only the turbulence collision is considered in the present model implementation. Further details on the two models can be found in the original papers by Luo and Svendsen (1996) and Prince and Blanch (1990). As for the bubble flow field, a single set of momentum equations was considered, thus obtaining an average flow field for all the bubble classes. A more general approach would consist in considering each bubble size group as a separate phase, but computational resource limitations and huge computational time requirements make it currently inapplicable to gas–liquid stirred tanks. The liquid and bubble velocity fields were computed solving the following momentum equation: j (lam) (turb) + ∇ · k k (  Uk ) + ∇ · k k Uk Uk + ∇ · k k jt k k = −k ∇pk + k k g + MI K . (4) In all the cases, the momentum transfer between the phases, MI K , included the contribution of drag forces, FD , while the contribution of turbulent dispersion, Ftd , virtual mass, Fvm and lift forces, Flift , was evaluated only in one simulation, based on the mono-disperse assumption. The following form for each of the above mentioned forces was considered: 3 CD   |Ug − Ul |(Ug − Ul ), 4 d g l

 DUg DUl − , Fvm = 0.5g l Dt Dt

FD =

Flift = 0.5g l (Ug − Ul )∇Ug ,

⎞ ⎠ d,

2109

(2)

where fbv represents the breakage volume fraction,  is the size ratio between the eddies in the inertial subrange of isotropic turbulence and of the bubble, min is equal to the ratio of the minimum size of eddies in the inertial subrange of isotropic turbulence and of the bubble diameter dbj ,  is the

Ftd = −0.1l kl ∇l .

(5)

In the MUSIG simulation the drag force was evaluated on the basis of the local bubble diameter. For the drag coefficient, a correlation similar to those adopted by Khopkar and Ranade (2006) and Wang et al. (2006) was selected:



 24 2 1/2 , (6) (1 + 0.15Re0.687 ) , E CD = max o b Reb 3

2110

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

where Reb is the bubble Reynolds number and Eo is the Eötvös number; therefore, depending on the local flow regime, the CD was the maximum value between those calculated by the Schiller and Nauman (1933) or the Ishii and Zuber (1979) correlations. 5. Computational domain and solution procedure The general purpose CFX-4.4 code (AEA Technology, 1999) was adopted for the simulations. The computational domain was built and discretised adopting Fortran user routines; it consisted of half vessel in azimuthal direction, included all the geometrical details of the system and was divided into two blocks for dealing with the relative motion of the rotating impeller and the steady baffles (by using the “sliding grid” algorithm). The domain could be limited to in the azimuthal direction due to the flow periodicity. Two block-structured grids were considered, of about 123 000 and 371 000 cells, respectively. In particular, grid refinement was performed along the angular direction (72 instead of 24 cells), while maintaining the same grid size in the other two directions (92×56 cells along the axial and radial coordinates, respectively). Still and ungassed liquid was considered as the initial condition for the simulations. The bubbles fed from the sparger were modelled adopting the inlet boundary conditions on the sparger top surface. The vessel top was modelled as a “degassing boundary”, i.e. a stationary, flat frictionless wall, at which the bubbles disappear due to a mass sink term depending on the outward normal velocity and the local volume fraction, as m ˙ = g g Ac Ug ,

(7)

where Ug is the gas axial velocity across the lower boundary of the cell and Ac is its area. Long time steps corresponding to 60◦ of impeller rotations were chosen in all cases, since time step refinement did not produce different results as discussed previously (Montante et al., 2007; Scargiali et al., 2007). Convergence and steady state attainment were guaranteed by almost constant and very low residuals values (10−6 ), by identical gas inlet and outlet flow rates and by identical flow variables after two subsequent complete impeller revolutions. In no case the model equations parameters were adjusted. The MUSIG simulations were performed adopting the coarser grid and considering the bubble diameter comprised between the minimum and the maximum value of 0 and 15 mm, respectively, with equal diameter division and 16 bubble classes. As a result, the mean equivalent diameter for the smaller bubble class was equal to 0.74 mm and that of the larger class was equal to 14.60 mm. Due to the coalescence phenomena observed at the sparger orifice, the mean diameter of bubbles at the sparger exit was taken as equal to that of the 13th MUSIG class, corresponding to 11.7 mm. For each time step, that consisted of 30 iterations, 25 internal iterations were fixed for the MUSIG equations. As a result, the time required for accomplishing a time step in the MUSIG simulation was six times longer than that of the fixed bubble diameter simulations.

6. Results and discussion 6.1. Bubble size distribution The BSD in the lower part of the vessel, including the impeller, and the upper half of the vessel (from the impeller disk to the liquid free surface) were typically obtained collecting 300 images; in each image about 600–800 bubbles were detected depending on the operating conditions. To obtain a detailed picture of the BSD in the different vessel regions, 13 zones, whose location is depicted in Fig. 1(a), have been identified and analysed: an example of the BSD in a few of these zones is shown in Fig. 1(b)–(d) for a given experimental condition. The mean value of the bubble sizes over the whole dispersion at different gas flow rates and fixed rotational speed is given in terms of different diameter definitions in Fig. 2, namely: number diameter, d10 , Sauter diameter, d32 , and weight mean diameter, d43 . Although the above-mentioned experimental constraints have limited the range of QG variation, the results clearly show that the bubble dimensions increase at increasing gas flow rate. This result suggests that probably the effect of bubble coalescence is not negligible even at the low gassing rates considered in this work. The cumulative distributions of the bubble equivalent diameter obtained at fixed rotational speed (N = 7.5 s−1 ) and two different gas flow rates in the lower part of the vessel are shown in Fig. 3. As can be observed, about 50% of the bubbles are characterised by an equivalent diameter lower than 0.8 mm and 90% of them are smaller than 2 mm. Similar results were obtained in the vessel upper part. A slight increase in the gas flow rate corresponds to a moderate increase in bubble size, that is reasonably due to coalescence. In Fig. 4, the difference produced by rotational speed variation in the upper half vessel can be appreciated: the bubble size is slightly lower for higher rotational speed, as one would expect. The cumulative distributions of bubble equivalent diameter, area and volume are depicted in Fig. 5 for a specific experimental condition. The area, aj , and the volume, Vj , of bubbles in the j -class have been evaluated as abj = nbj

db2j 4

,

vbj = nbj

db3j 6

,

(8)

where nbj and dbj are the number of bubbles and their equivalent diameter in the j -class, respectively. The results shown in Fig. 5 highlight that, while the most numerous bubbles are characterised by a diameter lower than 1 mm, the area and the volume are mainly associated with bubbles larger than 2.5 and 3.5 mm, respectively. This feature is particularly important to consider for interpreting PIV measurements and CFD results, because the experimental PIV data are mainly related to the value of d10 , while the gas holdup is related to d43 and mass transfer rate is related to d32 . Overall, in the 13 different vessel regions mentioned previously d10 ranges from 0.75 to 1.1 mm, the Sauter diameter, d32 , from 1.68 to 2.60 mm and d43 from 2.33 to 3.6 mm. The values of d10 are shown in Fig. 6 and clearly indicate that

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

C

D

B

2111

E

F

G A

B

C

G

F

D

E

B (down)

A

0

1

2 3 db [mm]

4 0

1

G (up)

2 3 db [mm]

4 0

1

2 3 db [mm]

4

Fig. 1. (a) Zone identification, (b)–(d) typical bubble size distributions for N = 7.5 s−1 and QG = 0.02 vvm in different vessel regions.

5.0 Cumulative number [-]

1.0

db [mm]

4.0 3.0 2.0 1.0 0.0 0.00

0.8 0.6 0.4 0.2 0.0

0.01

0.02

0.03

0.04 0.05 QG [vvm]

0.06

0.07

0.08

0.0

2.0

4.0

6.0

db [mm]

Fig. 2. Results of the bubble size measurements at different gas flow rate and fixed impeller speed (N = 7.5 s−1 ). Squares: d10 ; triangles: d32 ; diamonds: d43 .

Fig. 3. Cumulative distribution of bubble equivalent diameter in the lower half of the vessel. N = 7.5 s−1 . Solid symbols: QG = 0.02 vvm; empty symbols: QG = 0.07 vvm.

the smaller bubbles are located in the impeller discharge region (d10 = 0.75 mm). Out of this region, coalescence phenomena prevail in the zones B and C below the impeller, where the bubbles’ residence time is probably longer than in the upper zones: this is most likely due to the downward flow of the liquid phase, whose velocity magnitude in the axial direction close to vessel wall is of the same order of magnitude of the bigger bubble rise velocity (UL is equal up to 30 cm s−1 ). In zone D the same

increase of the bubble dimension due to coalescence is apparent below and above the impeller. In the upper zone E, bubble entrainment from the surface is most likely responsible for the bubble size increase, while in the corresponding E zone at the bottom the bubble diameter decreases probably due to the impact of the jet with the vessel bottom. In the central G region, negligible differences can be observed below and above the impeller. Overall, the two-loop fluid dynamic similarity results

2112

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

in equipment configuration and working conditions, the main trends—namely, spatial distribution of mean bubble diameter, influence of impeller speed and gas flow rate on average bubble size—is generally in good agreement with the previous findings (Alves et al., 2002; Kamiwano et al., 2003; Laakkonen et al., 2005).

Cumulative number [-]

1.0 0.8 0.6 0.4 0.2

6.2. Two-phase fluid flow

0.0 2.0

0.0

6.0

4.0 db [mm]

Fig. 4. Cumulative distribution of bubble equivalent diameter in the upper half of the vessel. QG = 0.07 vvm. Solid symbols: N = 7.5 s−1 ; empty symbols: N = 3.3 s−1 .

Cumulative number [-]

1.0 0.8 0.6 0.4 0.2 0.0 2.0

0.0

6.0

4.0 db [mm]

8.0

Fig. 5. Cumulative distribution of bubble equivalent diameter (), area () and volume (). Upper half of the vessel, QG = 0.02 vvm, N = 7.5 s−1 .

d10 [mm]

1.2

0.8

0.4 A

B

C

D zones

E

F

G

Fig. 6. Bubble number diameter, d10 , in different vessel zones (see Fig. 1(a) for the zone identification). Solid symbols: lower half vessel; empty symbols: upper half vessel.

in a similar BSD. The values of d32 and d43 follow the same qualitative trends of d10 ; therefore, the corresponding figures are not reported to avoid repetitions. Though direct comparison of bubble sizes and BSD with the literature data was unfeasible because of slight differences

The gas and liquid turbulent flow fields were measured on a vertical plane mid-way between two baffles at N = 7.5 s−1 and QG =0.02 vvm. As with the bubble size measurements, the flow field on the whole plane resulted from separate measurements performed on half vessel elevations for improving the image resolution. The single phase velocity field was also measured at the same conditions, for comparison purposes. The mean velocity features of the two-phase flow field and their comparison with that of the single phase at similar operating conditions have been recently discussed by Montante et al. (2007). In this work, we focus on the axial and radial mean velocity of the liquid and bubbles and on the velocity fluctuations of the two phases, at specific locations on the measurement plane. In Fig. 7, the mean radial and axial velocity profiles of the liquid phase and the bubbles are compared at (approximately) the same axial locations (the gas and liquid cameras were not equal, and the measurement coordinates do not coincide perfectly). As can be observed, below (z/T = 0.30), above (z/T = 0.70) and very close to the impeller disk (z/T = 0.49) very small differences are found between the liquid and the bubble mean radial velocity for all the radial positions and the selected axial locations, while remarkable slip velocities are found in the axial direction, apart from the impeller region where both phases are pumped radially out of the impeller and the axial velocities are almost nil. Above and below the impeller, the slip velocity is not constant along the radius: it significantly decreases at some locations towards the vessel wall, a feature that most probably depends on the BSD along the radial direction, which is induced by the local flow field. The r.m.s. radial and axial velocity components of the gas and liquid at the same axial locations of Fig. 7, are shown in Fig. 8. As for the gas phase, below and close to the impeller the radial and axial velocity fluctuations profiles do not differ significantly, while more marked deviations appear above the impeller, where the axial component is about twice the radial one in most radial locations. In any case, the bubble r.m.s. velocities are always bigger than those of the liquid phase; these differences between the two fluctuating components are similar at all the selected elevations. Finally, the effect of the bubbles on the liquid turbulent characteristics can be appreciated by comparing the radial and axial r.m.s. components of the liquid phase in the single and two-phase conditions. For this purpose, the ratio of the liquid velocity fluctuation in ungassed (single phase) and gassed conditions is shown in Fig. 9, thus representing turbulence reinforcement when it is lower than 1 and turbulent dampening in the opposite case. As can be observed, the effect of

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

0.4

0.3

z/T=0.30

z/T=0.30

0.3

0.2

Vr.m.s./Vtip [-]

Vmean /Vtip [-]

0.3

0.1 0.0 -0.1

0.2 0.2 0.1 0.1

-0.2 -0.3 0.0

0.1

0.3

0.2

0.4

0.0

0.5

0.0

0.1

0.2

r/T [-] 0.4

0.4

0.5

0.3

0.4

0.5

0.3

0.4

0.5

0.3

z/T=0.49

z/T=0.49

0.3

0.2

Vr.m.s./Vtip [-]

Vmean /Vtip [-]

0.3 r/T [-]

0.3

0.1 0.0 -0.1 -0.2

0.2 0.2 0.1 0.1

-0.3 0.0

0.1

0.2

0.3

0.4

0.0

0.5

0.0

0.1

0.2

r/T [-] 0.4

r/T [-] 0.3

z/T=0.71

0.3

z/T=0.71

0.3

0.2

Vr.m.s./Vtip [-]

Vmean /Vtip [-]

2113

0.1 0.0 -0.1 -0.2

0.2 0.2 0.1 0.1

-0.3 0.0

0.1

0.2

0.3

0.4

0.5

r/T [-] Fig. 7. Axial (diamonds) and radial (squares) mean velocity components at three elevations. Solid symbols: bubbles; empty symbols: liquid.

bubbles on the liquid turbulent characteristics is generally noticeable, but not particularly marked. The major difference pertains the radial r.m.s. velocity above the impeller, but probably the shape change in the upper loop in the gassed systems (Montante et al., 2007) is responsible for the difference. Overall, the measurements of two-phase turbulent characteristics confirm that beside the well-known limitations of the RANS-k.-based CFD simulations on the turbulent treatment in single-phase flows (the three r.m.s. velocity components are generally different), a further approximation is introduced by the “homogeneous” assumption in the case of two-phase flows (the turbulent velocity fluctuations of liquid and bubbles are

0.0 0.0

0.1

0.2 r/T [-]

Fig. 8. Axial (diamonds) and radial (squares) r.m.s. velocity components at three elevations. Solid symbols: bubbles; empty symbols: liquid.

not coincident). The gas and liquid velocity fluctuations can be adopted as a benchmark for advancement in two-phase turbulence modelling, whose impact would be particularly important for gas–liquid flows, due to their role on the estimation of bubble break-up and coalescence rates. 6.3. CFD results The CFD simulations of the gas–liquid stirred vessel already performed adopting the mono-disperse assumption for bubble size had provided good overall predictions of the liquid and gas mean velocity field, although some discrepancies with

2114

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

0.4

z/T=0.30

Liquid r/T=0.46

0.2 U/Vtip [-]

r.m.s.SP/r.m.s.TP [-]

1.6

1.0

0.0

-0.2 0.4 0.0

0.1

0.2

0.3

0.4

0.5

-0.4

r/T [-]

0.0

0.2

0.4

0.6

0.8

1.0

0.6

0.8

1.0

z/T [-] z/T=0.49

0.4 Gas r/T=0.46

0.2

1.0 U/Vtip [-]

r.m.s.SP/r.m.s.TP [-]

1.6

0.4

0.0

-0.2 0.0

0.1

0.2

0.3

0.4

0.5

r/T [-] -0.4 0.0

r.m.s.SP/r.m.s.TP [-]

1.6

z/T=0.71

0.2

0.4 z/T [-]

Fig. 10. Mean axial velocity profiles of liquid and bubbles at r/T = 0.46. Symbols: PIV data. Lines: CFD simulations (mono-disperse assumption). Thick line: coarse grid and pressure boundary condition on top. Thin line: fine grid and degassing boundary condition on top. Dotted line: fine grid, degassing boundary condition and MI K from Eq. (5).

1.0

0.4 0.0

0.1

0.2

0.3

0.4

0.5

r/T [-] Fig. 9. Ratio of the r.m.s. velocity components of the liquid phase in ungassed and in gassed conditions at three elevations. Axial component: diamonds; radial component: squares.

respect to the PIV data were noticed, particularly in the gas axial velocity above the impeller (Montante et al., 2007). Before performing a comparative analysis of the results obtained with the two simulation approaches adopted in this work (namely, fixed bubble diameter vs. MUSIG model), the influence of the grid size, of the boundary conditions for the open surface at the vessel top and of inter-phase forces other than drag on the mean flow field will be further examined for the monodisperse case. In order to carry out strict evaluation of possible differences in the simulations results, we will focus our attention on the mean

axial velocity profile at a fixed radial position (r/T = 0.46) along the vessel height (Fig. 10). As can be observed, the profiles obtained by either considering the drag force only or including also the added mass, lift and turbulent dispersion forces do not differ significantly, the differences being just noticeable for both phases. Also a grid size refinement of about three times on the tangential direction does not change the results, as is proved by the almost coincident axial velocities below the impeller. The only important difference in the results is due to the open surface treatment, which was changed from “pressure boundary” to “degassing boundary” condition; it is noticeable in both the liquid and gas phase profiles in the upper vessel portion depicted in Figs. 10a and b, respectively. In particular, a more limited discrepancy between experimental and simulated gas velocities was obtained close to the vessel top, accompanied by a reduction in the computational time for achieving the solution. It is important pointing out that the agreement of the predictions with the experimental mean flow field of the two phases is well represented everywhere by the profiles shown in Fig. 10,

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

5.0

2115

0.3 Liquid r/T=0.425

3.0

Vmean /Vtip [-]

db [mm]

4.0

2.0 1.0

0.0

0.0 d10

d32

d43 -0.3 0.0

5.0

0.2

0.4

0.6

0.8

1.0

0.6

0.8

1.0

z/T [-]

4.0 Liquid r/T=0.252

2.0

0.4

1.0 0.0 d10

d32

d43

Vmean /Vtip [-]

db [mm]

0.6 3.0

0.2

0.0

5.0 -0.2

db [mm]

4.0

0.0

0.2

0.4 z/T [-]

3.0

Fig. 12. Comparison of PIV and CFD mean velocity profiles for the liquid. Diamonds: PIV, axial; squares: PIV, radial. Thick lines: MUSIG simulation; thin lines: mono-disperse bubble size simulation. (a) r/T = 0.43; (b) r/T = 0.26.

2.0 1.0 0.0 d10

d32

d43

Fig. 11. Experimental (solid symbols) and MUSIG-calculated (empty symbols) bubble size. (a) Impeller region; (b) below the impeller; (c) above the impeller.

except in the upper region close to the shaft where marked differences were observed between PIV data and CFD results. In particular, the CFD gas velocity vectors are practically vertical in the shaft region, while the PIV data exhibit a significant radial component, probably due to gas entrainment from the open surface. Though, neglecting this discrepancy seems acceptable since modelling the air entrainment observed in the experiments would require a specific treatment for the gas–liquid interface that is beyond the capability of the present simulations. Therefore, we do not concentrate our analysis on that specific region. After these preliminary analyses, the following conditions were selected for effecting the MUSIG simulations: coarser grid, drag force only for the interphase exchange terms, MI K , and degassing boundary condition for the treatment of the open top surface. The MUSIG simulations were performed with the objective of verifying the actual possibility to provide reliable information on BSDs; in fact, they are of crucial importance

for stirred gas dispersion, since the main objective in several industrial applications is to maximise the contact area between the two phases. The values of d10 , d32 and d43 in three significant vessel regions are compared in Fig. 11: in particular those relevant to zone A are shown in Fig. 11a, while those of the zones B–G below and above the impeller are depicted in Figs. 11b and c, respectively. As can be observed, the simulation provides reliable values of d10 in all the three regions, a different level of accuracy for the predicted d43 depending on the zone, while worse agreement is apparent for d32 . Overall theses results are similar and complementary with respect to those of Laakkonen et al. (2007). The influence of the MUSIG simulation approach on the flow field predictions can be observed in Figs. 12 and 13, where the mean axial and radial velocity profiles at selected locations (r/T = 0.26 and 0.43), for the liquid and gas flow field, respectively, are shown. As can be observed, the liquid axial and radial velocity profiles follow the experimental trends quite closely (Fig. 12). With the MUSIG approach, a slight improvement with respect to the mono-disperse assumption simulations can be observed at both the selected radial position, while the differences are more marked for the gas axial velocity shown in Fig. 13, especially above the impeller, where the gas axial

2116

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

simulations in this respect. Indeed, in industrial applications, for accurate mass transfer rate previsions the retention time of the bubbles and the contact area between the phases have to be known. Since the bubble size interval for the simulations was purposely by far wider than that obtained from the experiments, probably better attention should be paid to the influence of the bubble size interval, its discretisation method and the number of bubble classes on the prediction accuracy in order to devise a completely successful modelling approach.

0.3

Vmean /Vtip [-]

Gas r/T=0.427

0.0

7. Conclusions

-0.3 0.0

0.2

0.4

0.6

0.8

1.0

0.6

0.8

1.0

z/T [-] 0.6 Gas z/T=0.26

Vmean /Vtip [-]

0.4

0.2

0.0

-0.2 0.0

0.2

0.4 z/T [-]

Fig. 13. Comparison of PIV and CFD mean velocity profiles for bubbles. Circles: PIV, axial; triangles: PIV, radial. Thick lines: MUSIG simulation; thin lines: mono-disperse bubble size simulation. (a) r/T = 0.43; (b) r/T = 0.26.

velocity is much closer to the experimental profile with respect to the case of constant bubble diameter. Apart from the sparger region, the gas holdup predictions show very low gas concentration below the impeller, bigger local holdup in the upper half of the vessel, significant gas accumulation at the rear of the blades and modest accumulation on the rear of baffles. These features are in qualitative agreement with the visual observation of the experimental conditions, but a quantitative evaluation on the reliability of the calculations with respect to the prevision of the gas holdup will be possible as soon as our experimental investigation will be extended to local gas holdup measurements. To date, accurate prediction of gas holdup distribution in gas–liquid stirred tanks have been obtained by RANS-based CFD simulations, also under the assumption of constant bubble size (Khopkar and Ranade, 2006; Scargiali et al., 2007). Overall, this study shows that reliable predictions of the twophase mean velocity fields and bubble size can be attained. Notably, no parameter was adjusted in the adopted break-up and coalescence models. The systematic underestimation of the Sauter diameter is however apparent in all cases. Further investigation is, therefore, required to improve the quality of the

In this work we have given a combined view of the flow field and the bubble size distribution in a gas–liquid stirred vessel equipped with a Rushton turbine. The experimental part of the investigation has provided insight in the complex fluid dynamics of the systems and benchmark results for validation of CFD models. The RANS-based approach selected for the present simulations was shown to be suitable for obtaining realistic information on the mean flow field of the two phases. The MUSIG model implemented with the Luo and Svendsen (1996) and Prince and Blanch (1990) models for bubble break-up and coalescence is able to provide improved predictions of the gas mean flow field, reliable prevision of the bubble number diameter and of the weight mean diameter, while improvement in the Sauter diameter calculations is expected from further modelling refinements. The MUSIG model is computationally demanding since the time required for achieving the solution is considerably longer than with the two fluid model approach based on a single bubble diameter, but already affordable even with current single processor computers. It is a fully predictive method, since only an initial guess on the bubble size interval has to be provided, with respect to the case of the mono-disperse assumption, where a diameter (or bubble rise velocity) has to be fixed. Therefore, the computational approach presented in this work seems to be promising and feasible as a design tool for industrial gas–liquid stirred tanks, at least for low gas holdup conditions.

Notation a A Ac Bc C CD db d10 d32 d43 D Eo fbv

area of bubble, m2 collision cross-sectional area of bubbles, m2 area of the computational cell, m2 birth rate of bubbles off-bottom impeller clearance, m bubble drag coefficient, dimensionless bubble diameter, m bubble number (arithmetic mean) diameter, m volume-to-surface (Sauter) diameter, m weight mean diameter, m impeller diameter, m Eötvös number, dimensionless breakage volume fraction, dimensionless

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

FD Flift Ftd Fvm g gB H k MI K n nb N p QG r Reb S t T U U ut V vb Vtip z

drag force, N lift force, N turbulent dispersion force, N virtual mass force, N gravity acceleration, m s−2 break-up rate, s−1 ungassed liquid height, m turbulent kinetic energy, m2 s−2 momentum transfer forces, N concentration of bubbles of diameter db , m−3 number of bubbles, dimensionless agitation speed, s−1 pressure, N m−2 gas flow rate, vvm radial coordinate, m bubble Reynolds number, dimensionless source term, kg s−1 time, s vessel diameter, m mean velocity vector, m s−1 mean axial velocity, m s−1 turbulent velocity, m s−1 mean radial velocity, m s−1 volume of bubble, m3 impeller tip speed, m s−1 axial coordinate, m

Greek letters       

model constant, dimensionless turbulent kinetic energy dissipation, m2 s−3 size ratio, dimensionless density, kg m−3 surface tension, N m−1 contact time, s stress tensor, N m−2 volumetric fraction, dimensionless

Subscripts j k

bubble size class liquid or gas phase

Acknowledgements This work was financially supported by the Italian Ministry of University and Research (MIUR) and the University of Bologna under PRIN 2005 programme. The assistance of Dr. Daniele Fajner in carrying out the experimental programme is gratefully acknowledged. References AEA Technology. CFX-4.4 Flow Solver: User Guide. Computational Fluid Dynamics Services, AEA Industrial Technology, Harwell, Oxfordshire, UK, 1999. Alves, S.S., Maia, C.I., Vasconcelos, J.M.T., 2002. Experimental and modelling study of gas dispersion in a double turbine stirred tank. Chemical Engineering Science 57, 487–496.

2117

Aubin, J., Le Sauze, N., Bertrand, J., Fletcher, D.F., Xuereb, C., 2004. PIV measurements of flow in an aerated tank stirred by a down- and uppumping axial flow impeller. Experimental Thermal and Fluid Science 28, 447–456. Bakker, A., van den Akker, H.E.A., 1994. A computational model for the gas–liquid flow in stirred reactors. Chemical Engineering Research & Design 72, 594–606. Barigou, M., Greaves, M., 1992. Bubble-size distribution in a mechanically agitated gas–liquid contactor. Chemical Engineering Science 47, 2009–2025. Deen, N.D., Hjertager, B.H., 2002. Particle image velocimetry measurements in an aerated stirred tank. Chemical Engineering Communications 189, 1208–1221. Deen, N.D., Solberg, T., Hjertager, B.H., 2002. Flow generated by an aerated Rushton impeller: two-phase PIV experiments and numerical simulations. Canadian Journal of Chemical Engineering 80, 638–652. Hu, B., Pacek, A.W., Stitt, E.H., Nienow, A.W., 2005. Bubble sizes in agitated air–alcohol systems with and without particles: turbulent and transitional flow. Chemical Engineering Science 60, 6371–6377. Ishii, M., Zuber, N., 1979. Drag coefficient and relative velocity in bubbly, droptlet or particulate flows. A.I.Ch.E. Journal 25, 843–855. Kamiwano, M., Kaminoyama, M., Nishi, K., Shirota, D., 2003. The measurement of bubble diameter distributions and liquid side mass transfer coefficients in a gas–liquid agitated vessel using a real-time, high-speed image processing system. Chemical Engineering Communications 190, 1096–1114. Kerdouss, F., Bannari, A., Proulx, P., 2006. CFD modelling of gas dispersion and bubble size in a double turbine stirred tank. Chemical Engineering Science 61, 3313–3322. Khopkar, A.R., Ranade, V.V., 2006. CFD simulation of gas–liquid stirred vessel: VC, S33, and L33 flow regimes. A.I.Ch.E. Journal 52, 1654–1672. Khopkar, A.R., Aubin, J., Xuereb, C., Le Sauze, N., Bertrand, J., Ranade, V.V., 2003. Gas–liquid flow generated by a pitched-blade turbine: particle image velocimetry measurements and computational fluid dynamics simulations. Industrial & Engineering Chemical Research 42, 5318–5332. Laakkonen, M., Moilanen, P., Miettinen, T., Saari, K., Honkanen, M., Saarenrinne, P., Aittamaa, J., 2005. Local bubble size distributions in agitated vessel comparison of three experimental techniques. Chemical Engineering Research & Design 83, 50–58. Laakkonen, M., Moilanen, P., Alopaeus, V., 2007. Modelling bubble size distribution in agitated vessels. Chemical Engineering Science 62, 721–740. Lane, G.L., Schwarz, M.P., Evans, G.M., 2004. Numerical modelling of gas–liquid flow in stirred tanks. Chemical Engineering Science 60, 2203 –2214. Luo, H., Svendsen, H.F., 1996. Theoretical model for drop and bubble breakup in turbulent dispersions. A.I.Ch.E. Journal 42, 1225–1233. Machon, V., Pacek, A.W., Nienow, A.W., 1997. Bubble sizes in electrolyte and alcohol solutions in a turbulent stirred vessel. Chemical Engineering Research & Design 75, 339–348. Montante, G., Magelli, F., 2005. Modelling of solids distribution in stirred tanks: analysis of simulation strategies and comparison with experimental data. International Journal of Computational Fluid Dynamics 19, 253–262. Montante, G., Paglianti, A., Magelli, F., 2007. Experimental analysis and computational modelling of gas–liquid stirred vessels. Chemical Engineering Research & Design 85, 647–653. Morud, K.E., Hjertager, B.H., 1996. LDA measurements and CFD modelling of gas–liquid flow in a stirred vessel. Chemical Engineering Science 51, 233–249. Prince, M.J., Blanch, H.W., 1990. Bubble coalescence and break-up in airsparged bubble columns. A.I.Ch.E. Journal 36, 1485–1499. Rousar, I., van den Akker, H.E.A., 1994. LDA measurements of liquid velocities in sparged agitated tanks with single and multiple Rushton turbines. Institution of Chemical Engineers Symposium Series, vol. 136, pp. 89–96. Scargiali, F., D’Orazio, A., Grisafi, F., Brucato, A., 2007. Modelling and simulation of gas–liquid hydrodynamics in mechanically stirred tanks. Chemical Engineering Research & Design 85, 637–646.

2118

G. Montante et al. / Chemical Engineering Science 63 (2008) 2107 – 2118

Schiller, L., Nauman, A., 1933. Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung. Zeitschrift Vereines deutscher Ingenieure 77, 318–320. Shäfer, M., Wächter, P., Durst, F., 2000. Experimental investigation of local bubble size distributions in stirred vessels using phase doppler anemometry. In: van den Akker, H.E.A., Derksen, J.J. (Eds.), 10th European Conference on Mixing. Elsevier, Delft, The Netherlands, pp. 205–212. Takahashi, K., Nienow, A.W., 1993. Bubble sizes and coalescence rates in an aerated vessel agitated by a Rushton turbine. Journal of Chemical Engineering of Japan 26, 536–542.

van den Akker, H.E.A., 2006. The details of turbulent mixing process and their simulation. Advances in Chemical Engineering 31, 151–229. Venneker, B.C.H., Derksen, J.J., van den Akker, H.E.A, 2002. Population balance modelling of aerated stirred vessels based on CFD. A.I.Ch.E. Journal 48, 673–685. Wang, T., Wang, J., Jin, Y., 2006. A CFD–PBM coupled model for gas–liquid flows. A.I.Ch.E. Journal 52, 125–140.