Available online at www.sciencedirect.com

Engineering

Procedía

ELSEVIER

Procedía Engineering 42 (2012) 974 - 991

www.elsevier.com/locate/procedia

20th International Congress of Chemical and Process Engineering CHISA 2012 25 - 29 August 2012, Prague, Czech Republic

Hydrodynamics of airlift reactor with internal circulation loop: Experiment vs. CFD simulation

P. Lestinsky a*9 P. Vayrynenb, M. Vecera, K. Wichterlea

aDept. of Chemistry, VSB-Technical University of Ostrava 17. listopadu 15, Ostrava 70833, Czech Republic bDept. ofMaterial Science andEngineering School of Chemical Technology, Aalto University, Vuorimiehentie 2, Espoo, Finland

Effect of geometrical parameters on two phase hydrodynamics of airlift reactor is the main topic of present paper. Laboratory scale apparatus with internal circulation loop consists of concentric draft tube, in which a gas bubbles rising. Setup of draft tube inside of reactor is important geometry parameter and has big influence on two phase hydrodynamics. In experiment was studied influence of changes diameter of draft tube to hydrodynamics in airlift reactor. Results of experiments (liquid velocity and gas hold-up) were compared with the simple CFD simulations performed in COMSOL Multiphysics 3.5a. For each point of gas volumetric flow in simulation, were determined conditions of bubble diameter and bubble drag coefficient. Although bubble break-up and coalescence were neglected, the results of numerical simulation are in pretty good agreement with experimental data.

© 2012 Published by Elsevier Ltd. Selection under responsibility of the Congress Scientific Committee (Petr Kluson)

Keywords: Bubbles; airlift reactor; draft tube; CFD simulation

Nomenclature

Ag cross-section area to gas input (m2) CR turbulence modeling constant, CR = 0.09(1)

a* Corresponding author. Tel.: +420 597 324 230; E-mail address: pavel.lestinsky@vsb.cz .

Abstract

1877-7058 © 2012 Published by Elsevier Ltd. doi:10.1016/j.proeng.2012.07.482

C9i turbulence modeling constant, C9i = 1.44 (1)

C92 turbulence modeling constant, Cv2= 1.92(1)

ck bubble induced turbulence modeling constant, Ck = 0.01-1 (1)

r bubble induced turbulence modeling constant, C9 =1-1.9(1)

Co drag coefficient (1)

db bubble diameter (m)

F volume force (interaction force) (N/m3)

g gravity acceleration, g = 9.81 (m2/s)

Ahr,d difference of level in reverse U-tube manometer in riser and downcomer (m)

Kb,T loss coefficient at the bottom and at the top (1)

M molecular weight of gas, Man = 0.029 (kg/mol)

NP gas mass flux (kg/(m2 s))

Ap pressure difference (Pa)

P hydrostatic pressure (Pa)

Qg input gas volumetric flow (m3)

R universal gas constant, R = 8.314 (J/(mol.K))

sk additional source term (kg/(m3.s))

T temperature of liquid, T = 288 (K)

Ui liquid velocity (m/s)

Ui superficial liquid velocity (m/s)

Ug gas velocity (m/s)

Ug superficial gas velocity (m/s)

Uslip slip velocity (m/s)

Azr.d difference of height between delivery point from reverse U-tube manometer in riser and

downcomer (m)

9 dissipation rate of turbulent energy (m/s3)

eg gas hold-up (1)

egr,gd gas hold-up in riser and downcomer (1)

El liquid hold-up (1)

ni viscosity of liquid (Pa.s)

turbulentviscosity (Pa.s)

Pg density of gas (kg/m3)

Pd density of gas-liquid dispersion (kg/m3)

pi density of liquid (kg/m3)

k kinetic energy (m2/s2)

Gk turbulence modeling constant, ck =1(1)

c9 turbulence modeling constant, c9 =1.3(1)

4 loss coefficient (1)

1. Introduction

The Airlift reactors are finding increasing applications in chemical industry, biochemical fermentation and biological wastewater treatment processes [1,2]. Generally there are two types of airlift reactors: reactor with internal and with external circulation loop.

For design of an airlift reactor with internal circulation loop, it is necessary to have accurate estimates of the phase hold-ups and velocity in the riser and in the downcomer. Several literature studies have focused on the estimation of these hydrodynamics parameters [3-5]. The velocities of the liquid in the downcomer and in the riser are dependent on the frictional losses, which are determined by geometry of the airlift reactor and operating condition.

Several recent publications have established the potential of CFD simulation for describing the hydrodynamics of airlift reactors [6-9].

l.l.Airlift reactor

The upward flowing liquid in the center of the column can be stabilized over the whole height by inserting a draft tube. This type of column is called a bubble column loop reactor. Bubble columns loops are commonly known as Airlift reactors. Air is injected through a sparger at or near the bottom of the up-flow section, the so called riser. The reduced density of the two phase mixture in the riser causes a circulation in the loop in the same way as with an air-lift pump. High gas flow rates are possible. Airlift reactors are widely used in process industries in recent years because of the following advantages: simple construction without moving parts, good mass heat transfer characteristics, efficient mixing with low energy consumption, and simple operation with low cost [10].

1.2. Bubble/low model

Essentially two basic approaches to simulation of flow dynamics of two-phase gas-liquid systems are known. The first is an approach where both the liquid phase and the gas-phase motion are considered as a homogeneous. These two-fluid approximations are presented in Eulerian representation and thus referred to as Euler-Euler simulation. The second approach treats only the liquid-phase motion in an Eulerian representation and computes the motion of the dispersed gas-phase elements in a Lagrangian way by individually tracking them on their way through the liquid body. This approach has been termed Euler-Lagrange representation.

Two fluids Euler-Euler model is basic macroscopic model for two-phase fluid flow. The derivation of the model equations for the two-phase bubbly flow starts with the assumption that both phases can be described as continua, governed by the partial differential equations of continuum mechanism. The phases are separated by an interface, which is assumed to be a surface. At the interface, jump condition for the conservation of mass and momentum can be formulated. Bubble flow is simplification of two-fluid flow model with following assumptions:

density of gas is negligible compared with density of liquid

motion of gas bubble relative to liquid is determined by balance between viscous shear and pressure

pressure field are equivalent for both phases

The sum of the momentum equation for both phases is result from these assumptions. The sum of equation content the momentum equation for liquid velocity, the continuity equation and the transport equation for gas hold-up.

The interaction between the various phases follows from the distribution of the pressure and shear stress at the interface. However, the details of these are not known and cannot be resolved in the Euler-Euler approach, as individual interfaces are not represented in this theoretical framework. Instead, the interaction forces have to be formulated separately and fed back to the momentum equation [11].

1.3. Interfacialforces

The magnitude of F describes the interaction force between the continuous and the dispersed phase. The plus marker is given for gas phase and the minus marker is given for liquid phase. Sokolichin [12,13] considered, that only the pressure and the gravity forces has been effect on stationary gas bubble in stationary liquid in two-phase flow. Since there is usually a relative motion between the bubble and the liquid, the liquid flow around individual bubbles leads to local variations in pressure and shear stress. Usually three different contributions for the interaction force term are taken into account, a friction force term -Ff, an inertia force term - Ft and a lift force term - Ft, leading to the following approximation for the force of interaction:

where Cw is friction coeffcient. Schwarz [14] determined values of Cw = 5.104 kg/(m3.s) for usUp = 0,2 The inertia force is given equation:

F = Ff+Fi + Fi The friction force is calculated by:

F/= -sg Cw(ug- ui)

Fi = -eg Cipi d(ug-ui)/dt

where Ct is coefficient corresponding of amount of liquid, which is moved by gas bubble. Cook and Harlow [15] provided from experiment, that Ct = 0,25.

Last term in equation (15) is the lift force, which Thomas [16] express by relation:

F, = -eg Ct pi (ug - ui)^(vxui) (4)

where Ct is coefficient ofbuoyancy. Delnoij [17] provided this coefficient as Ct = 0,53. Boisson [18] suggested compound of friction force and inertia force together and replace them by the drag force. The drag force is calculated by relation:

Fd = 3/4 Co/db sg si pi usup (ug- ui) (5)

This approach (i.e. the drag force and the lift force were used for simulation) used another authors [1921] and they has been good agreement numerical results with experimental data. On the other hand Joshi [22] used relation:

Fo = (Sg (Pg -pi)g(Ug- u,))/(ug - u,) (6)

On the basic of these examinations, Sokolichin [23] used for determined interfacial force only the lift force and the drag force. Without the lift force, the stationary bubbles hold in the stagnant liquid and without drag force, the bubbles has infinite acceleration.

The COMSOL Multiphysics 3.5a make it possible input of interfacial force fas a constant. The magnitude of the interfacial force has a great influence to the resulting liquid velocity and gas hold-up. The equation (6) was used for our simulation, because it does not contain a constant as drag coefficient and bubble diameter. These constants obtain equation of the slip velocity and we obtained the two independent variables (slip velocity and interfacial force). In simulation we used the simplification of continuity equation; therefore we simplified equation (6) to term:

FD = (ps-pi)g (7)

The interfacial force after substitution of constants is:

F = -9778N/m3 (8)

Azzopardi [11] published a sum of interfacial forces occur between both phases. In our simulations we calculating with the simplification of continuity equation and many variables at interfacial forces is not necessary for us.

1.4. Slip 'velocity and drag coefficient

The slip velocity is name for relative velocity between gas and liquid phase. It is a one of the main parameters of multiphase flow. Together with interfacial force belong to variables, which could affect the result of numerical simulations. The slip velocity can be calculated as:

us,ip = -((4 db Ap)/(3 cD pi)/5 (9)

The most important parameters are the bubble diameter - db and the drag coefficient - cD. Another parameters is the hydrostatic pressure - but his values cannot be affected. The bubble diameter is assumed as a constant (equivalent sphere diameter), over all reactor. This assumption has been used also by other authors [5,24]. In some cases of numerical simulations the bubble diameter was represented by a bubble size distribution. Khrisna [25] uses such approach for three-phase flow reactor with Fisher-Tropsh synthesis, when attention was focused on the kinetics and the mass transfer rate.

Bubble diameter in present numerical simulation follows paper by Johansen [26], who uses eq. (10) as bubble size diameter for case when gas is injected into the bottom part through orifice.

db = 0.35 (Qg2/g)°

Drag coefficient is another important parameter for slip velocity determination. There are many approaches how to determine the drag coefficient which has been published in literature. Clift [27] published beside standard drag curve for solid spherical particle also empirical lines denotes behaviour of drops and bubbles in contaminated and pure liquid systems. Graphical dependence of the drag coefficient vs. the Reynolds number can be seen in Fig. 1.

Fig. 1. Drag coefficient as a function of Reynolds number for water drops in air and air bubbles in water, compared with standard drag curve for rigid spheres

Schiller [28] made a correlation of the drag coefficient vs. the Reynolds number for spherical particles (11), (12).

0<Re < 1000 cD = 24/Re (1 + 0.15Re0687) (11)

Re > 1000

cD = 0.44

Delnoij [29] used this consideration for numerical simulation of bubble column reactor. Application of these equations is safe for bubbles with diameter smaller than 1 mm (spherical shape).

Empiric formula covering effect of surface tension was introduced by Johansen [26]. They used the Eôtvôs number for prediction of the drag coefficient for larger bubbles with non-spherical shapes:

cD = 0.622/(1/Ed + 0.235)

Ranade [30] and Bhanu [31] used this from for calculation of bubble drag coefficient in ultrapure water. The values of drag coefficient corresponding to the curve (13) are denoted as "bubble in pure water" (dashed line in Fig. 1) inthe region of Re > 1000.

1.5. Numerical model

Nowadays, when supercomputer exists, the ways for simulation multiphase flow, counting individual movement and behavior of each bubble (break-up and coalescence, collision dynamics, oscillation of shape) are open. Unfortunately, no commercial software covers all of those processes and specific in-house codes are needed [32].

Anyway, using software such FLUENT or COMSOL for simulation of multiphase flows seems to be not much effective, but it can bring useful information. Simulation of internal loop airlift reactor under conditions similar with experiments was carried out using COMSOL. The governing equations calculated in this case were as follows:

The momentum equation for liquid phase:

s, pi (du)/dt + si pi u,Vu, = -Vp + V[s, fa, + (Vu, + Vu,T)] + s,p,g±F (14)

where is turbulent viscosity defined as:

nT=p,C,1{!/V (15)

where k is the turbulent kinetic energy, y is the dissipation rate of turbulent energy and their values are calculated by closure equations:

p, dk/dt - V[fa + Vk] +p,u,Vk = 1/2 (Vu, + Vuf /-p,<p+Sk (16)

p, dq/dt - V [fa + qT/a9) Vcp] +p, u, Vcp = 1/2 C9l y/k (Vu, + Vuf)2 - p, Cf2 ?2/k + y/k C9 Sk (17)

where Sk is the additional source term comprehensive of turbulence induced by bubble:

sk = -Ck eg Vp usiip (18)

The continuity equation for two phase systems is:

d/dt (s,p, + sgpg) + V(s,p, u, + sgpgUg) =0 (19)

The momentum equation for gas phase:

(5egpg)/at+V(egpgug) = 0 (20)

Invoking the equation of state for each phase:

p, = const.,pg=(Mp)/(R T) (21)

and following closure relation:

s, + sg=l (22)

We obtain a closed system of ten differential and algebraic equations, which describes the dynamics of the two-phase flow and can be solved numerically when specific conditions are reached [33]. Solving of such system on relatively large computing domain needs nonstandard hardware or specific computing conditions (clusters for parallel computing).

Other option is provide initial calculation with simplified model. When low gas fraction occurs in the system, the value of s, is close to one. If we neglected the variation of s, due to the presence of gas in the continuity equation for the liquid phase eq. (19) becomes:

(ds,)/dt+ V(s,u) =0 (23)

It will be transformed into (with assumption e~l):

Vu, = 0 (24)

The gas velocity is then calculated as a sum of liquid velocity, slip velocity and drift velocity:

ug = ui + usup + u ¿rift (25)

where the slip velocity usiip is calculated by equation (9), and drift velocity udrift is the additional condition for employed turbulent model, calculated as:

Udrift = -nr/Qi Veg/ei (26)

The continuity equation can be used for calculation local distribution of gas hold-up.

1.5.1. Boundarycondition

2D axial symmetry was used. Gas input condition at the bottom is given as gas mass flux calculated as product of superficial gas velocity and effective density.

Np = -n (P; Ug) (27) where

Pg=Pgeg (28)

Ug = Q/Ag (29)

Gas outlet boundary condition is used in top, where the gas phase flows outward with the gas velocity - ug. Condition for liquid at same position is treated as outlet with defined pressure. Boundary condition at left hand wall is axial symmetry.

Rests of boundaries are adjusted as walls with a no slip condition. There are generally two approaches for adjusting boundary conditions at solid walls implemented into commercial software. Both are using empirical models to calculate the thin laminar boundary layer adjacent to the solid walls. The first approach is "Wall offset", the second is "Wall offset in viscous units". More details can be found in user manual COMSOL. Wall offset setup has been used here.

Bubble diameters calculated by equation (10) and drag coefficients for contaminated systems (see Fig. 1) related for inlet gas volumetric flow are reported in Tab. 1.

Table 1. Bubble diameter and drag coefficient as function inlet gas volumetric flow

Q db cD

[m3/s] [mm] [ - ]

1.11 * 10-04 5.8 2.30

1.67 * 10-04 6.8 2.46

2.22 * 10-04 7.7 2.58

2.78 * 10-04 8.4 2.67

3.33 * 10-04 9.0 2.75

3.89 * 10-04 9.6 2.81

4.44* 10-04 10.1 2.87

5.00 * 10-04 10.6 2.91

5.56 * 10-04 11.1 2.96

1.5.2.Meshgrid

Triangular mesh elements were used for creating computational domain. Distance between the elements was 0.007 m. Distance between mesh elements at problematic places (both ends of draft tube) was 0.001 m.

1.5.3. Numerical Stabilization

Isotropic diffusion and streamline diffusion, or both can be added for the momentum equations as well as the gas transport equation. The momentum equation is stabilized using streamline diffusion by default. Gas transport equation is stabilized by default using both terms. Note, a suitable scale for the effective gas density and the bubble number density must be specified. More information about numerical stabilization is in the manual COMSOL.

1.5.4. Solving the model Default element type: Analysis type: Solver:

Time stepping:

Segregated groups:

Method:

Lagrange - P2Pi Transient

Time dependent segregated time 0-60 second Relative tolerance 0.001 Absolute tolerance 0.0001

1 - ul, vl (liquid velocity, radial and vertical component), p (pressure), rhogeff (effective gas density)

2 - logk (logarithm of turbulent kinetic energy), logd (logarithm of turbulent dissipation rate)

BDF backward differentiation formulas (it uses variable-order variable-step-size)

2. Experimental

Schematic view of experimental apparatus is in Fig. 2. Pachuca tank is the name used for industry scale clone of such apparatus. It is typical equipment in hydrometallurgy for leaching of key component from ore. Our apparatus consist of glass vessel and has diameter of 15 cm and height of 150 cm. The draft tube is inside. Gas hold-up is measured by pressure difference in reverse U-tube manometer. A Pitot tube has been used for liquid velocity measurements. Geometry of draft tube is varied in length and diameter of draft tube (length of 0.8 and 0.6 m and diameter of 0.1, 0.08 and 0.06 m). The next geometrical variable has been distance of draft tube from bottom, which ranges from 15 mm to 65 mm.

In experiments has been used two-phase system air-water.

Fig. 2. Experimental apparatus; 1-Laboratory airlift reactor, 2-Draft tube, 3-Gas inlet and flow-meters, 4-Reverse U-tube manometer for measure pressure difference (to interpretation of gas hold-up), 5-Reverse U-tube manometer for measure liquid velocity by Pitot tube, 6-Gas bubbles distributor, 7-Reading distance of drafttube from bottom, 8-Escape valve;

Three types of experimental set-ups are reported here. The list of set-up parameters is given in Tab. 2. The distance between draft tube and bottom was calculated by:

hB = 0.4 dinternal (30)

Table 2. Geometry and position of draft tube in the experimental reactor

hB ^internal ^external Adown / Arise hdraft tube

(m) (m) (m) (i) (m)

type 1 0.02 0.052 0.06 7.05 0.8

type 2 0.03 0.072 0.08 3.07 0.8

type 3 0.035 0.092 0.10 1.49 0.8

Main measured quantity is pressure difference resulting to gas hold-up. The gas hold-up is defined as volume fraction of gas phase in the gas-liquid dispersion in the reactor.

egr,d=pi/(pi-pg) hr,/ zr,d (31)

The value of liquid circulation velocity Ulr is one of the most important design parameter for airlifts or Pachuca reactors. Liquid circulation has influence on to the gas hold-up in the vessel, to heat and mass transfer coefficient and to the extent of mixing in the reactors [2]. The liquid velocity has been measured by Pitot tube. Pressure difference has been read from the inverse U-tube manometer:

Ulr = (2gpL Ah/pDr)a5

Gas holdup and liquid circulation velocity have been calculated from the experimental data and results are discussed in the section Results.

3. Results and discussion

The model was constructed in the COMSOL Multiphysics 3.5a, it means draw geometry, set-up boundaries and variables, mesh grid and solver set-up. This model was exported to the Matlab and there was generated loop simulation for all gas input flux, bubbles diameters and drag coefficients. Postprocessing of results and comparison between experimental data and numerical solutions has been carried out also in the Matlab. The values of liquid velocity and gas hold-up were calculated as integrated variables per volume of riser and downcomer:

U,= (i2xru, dU,)/(idV) ,eg=(i2xreg dsg)/(idV) (33)

The main comparison between numerical solution and experimental data is figured out in following figures. The figures are dividing by results of liquid velocity in riser and gas hold-up in riser and in downcomer. These parameters have been plotted for gas volumetric flow and superficial gas velocity in the draft tube.

3.1. Liquid circulation -velocity

The liquid velocity in riser increases with increasing gas volumetric flow. Influence of draft tube diameter on the liquid velocity in riser can be seen in Fig. 3. With decreasing diameter of draft tube the liquid velocity increases. Pretty good agreement of experiments and numerical simulations can be seen for prediction of liquid circulation velocity with increasing gas inlet volumetric flow

• type 1 CFD

2 3 4 5 6

Q -gas volumetric flow [m3/s] x10~4

Fig. 3. Liquid velocity in riser as a function of gas volumetric flow for various diameter of draft tube

Fig. 4. Liquid velocity in riser as a function of superficial gas velocity in riser for various diameter of drat tube

3.2. Gashold-up

Gas hold-up is one of important parameters of airlift reactors or bubble columns. Gas hold-up in reactor increases with increasing gas volumetric flow. Linear trend of gas hold-up in riser can be seen in Fig. 5. The difference between experimental data and CFD simulation has been caused by assumption of homogeneous distribution of bubble diameter in the reactor and neglecting of bubble coalescence, which corresponds with original setup ofbubbly flow model build in COMSOL.

Fig. 5. Gas hold-up in riser as a function of gas volumetric flow for various diameter of drat tube

Heterogeneous and Slug bubbly flow regimes are main flow regimes which occurs during experiments carried out in our apparatus. In these regimes coalescence takes place and bubble sizes are not uniform. Homogeneous and transient and heterogeneous bubble flow regimes which are typical for bubble columns [33], did not occurs at our conditions.

In small draft tube diameters we have used, the Slug bubble flow regime is dominant, it means, that resulting bubbles produced by coalescence of initial bubbles are the same size as a diameter of draft tube. Such behavior was visually observed in our experiments.

The influence of diameter of draft tube on the gas hold-up in riser can be seen in Fig. 6. With increasing of draft tube diameter gas hold-up in riser decreases for same values of superficial gas velocity in riser.

Fig. 6. Gas hold-up in riser as a function of superficial gas velocity in riser for various diameter of drat tube

Gas hold-up in downcomer increases with increasing gas volumetric flow. The influence of diameter of the draft tube to the gas hold-up in downcomer can be seen in Fig. 7, 8. Gas hold-up in downcomer increases with increasing of draft tube diameter for same values of inlet gas volumetric flow or superficial gas velocity. The differences between experiments and CFD simulation for higher values of gas volumetric flow in geometry of type 2 and 3 are caused by total circulations of bubbles. The similar problem was described by van Benthum [4].

Smaller bubbles follow liquid stream and fall down though the downcomer. When the liquid circulation velocity in downcomer is high enough (more than 0.3 m/s) the total circulation of small bubbles is presented. The used model is not capable to correctly describe such behavior. On Fig. 8 can be seen small increasing of gas hold-up in downcomer for geometry of type 1 with increasing superficial gas velocity in riser. The total circulations ofbubbles were not observed in experiments for this geometry.

Q0 - gas volumetric flow [m3/s] x | q 4

Fig. 7. Gas hold-up in downcomer as a function of gas volumetric flow for various diameter of drat tube

Fig. 8. Gas hold-up in downcomer as a function of superficial gas velocity in riser for various diameter of drat tube

4. Conclusions

Experimental investigation of geometrical parameters on to the two phase hydrodynamics effect in airlift reactor was carried out. Three different geometries of airlift with internal circulation loop were studied. Numerical simulation of experimental conditions using commercial software COMSOL was made. Experimental results and numerical simulations were compared when pretty good agreement (± 10%) was found for prediction of gas holdup and liquid circulation velocity for small inlet gas volumetric flows and medium size and large draft tube. Data for high inlet volumetric flow and smaller draft tube are loaded by significant error. This is a task for no adequate simplification of numerical model, which, on the other hand, allow us to solve such complex problem with standard personal computer almost in real

Acknowledgements

The Generous support of this study by SP2012/196 and GACR under projects 104/07/1110, 104/09/0972 is gratefully acknowledged.

References

[1] Blenke H. Loop Reactors. AdvBiochemEng 1979;13:122-213.

[2] Chisti MY. Airlift Bioreactors. London: Elsevier Applied Science; 1989.

[3] Heijnen JJ, Hols J, van der Lans RGJM, van Leeuwen HLJM. A simple hydrodynamics model for the liquid circulation velocity in full-scale two- and three-phase internal airlift reactor operating in the gas recirculation regime. Chemical Engineering Science 1997;52:2527-2540.

[4] van Benthum WAJ, van der Lans RGJM, Heijnen JJ. Bubble circulation regimes in an internal-loop airlift reactor. Chem Eng Sci 1999;54:3995-4006.

[5] Oey RS, Mudde RF, PortelaLM. Simulation of a slurryairliftusing a two-fluidmodel. Chem EngSci2001;56:673-681.

[6] van Baten JM, Ellenberer J, Krishna R. Using CFD to describe the hydrodynamics of internal air-lift reactors. The Canadian J Chem Eng 2003;81:660-668.

[7] Oey RS, Mudde RF, van der Akker HEA. Numerical simulation of an oscillating internal-loop airlift reactor. Can J Chem Eng 2003;81:684-691.

[8] Huang Q, Yang Ch, Yu G, Mao ZS. CFD simulation of hydrodynamics and mass transfer in an internal airlift loop reactor using a steady two-fluid model. Chem Eng Sci 2010;65:5527-5536.

[9] Simcik M, Mota A, Ruzicka MC, Vicente A, Teixeira J. CFD simulation and experimental measurement of gas holdup and liquid interstitial velocity in internal loop airlift reactor. Chem Eng Sci 2011;66:3268-3279.

[10] van der Lans RGJM. Hydrodynamics of a bubble column reactor. PhD. thesis. Delft; 1985.

[11] Azzopardi BJ, Mudde RF, Morvan H, Yan YY, Zhao D. Hydrodynamics of gas-liquid reactors: normal operation and upset condition. 1st ed. John Wiley & Sons Ltd; 2011.

[12] Sokolichin A, Eigenberger G. Gas-liquid flow in bubble columns and loop reactors: Part I. detailed modeling and numerical simulation. ChemEngSci 1994;49:5735-5746.

[13] Becker S, Sokolichin A, Eigenberger G. Gas-liquid flow in bubble columns and loop reactors: Part II. comparison of detailed experiments and flow simulations. Chem Eng Sci 1994;49:5747-5762.

[14] Schwarz MP, Turner WJ. Applicability of the standard k-e turbulence model to gas-stirred baths. Appl Math Model 1988;12:273-279.

[15] Cook TL, Harlow FH. Vortices in bubblytwo-phase flow. Int J Multiphase Flow 1986;12:35-61.

[16] Thomas NH, Auton TR, Sene K, Hunt JCR. Entrapment and transport of bubbles in transient large eddies in multiphase turbulent shear flow. International conference on the physical modelling ofmulti-phase flow, Coventry, UK, 1983:169-184.

[17] Delnoij E, Kuipers JAM, van Swaaij WPM. A three dimensional CFD model for gas-liquid bubble columns. Chem Eng Sc; 1999;54:2217-2226.

[18] Boisson N, Marlin MR. Numerical prediction of two-phase flow in bubble columns. Int J Num Methods Fluids 1996;23:1289-1310.

[19] Krishna R, Urseanu MI, van Baten JM, Ellenberger J. Influence of scale on the hydrodynamics of bubble columns operating in the churn-turbulent regime: experiments vs. Eulerian simulation. Chem Eng Sei 1999;54:4903-4911.

[20] Deen NG, Solberg T, Hjertager BH. Numerical simulation of the gas-liquid flow in square cross-sectioned bubble column. CHISA 14th Int. Congress ofChemieal andProeessEngineering, Prague, CZE, 2000:16.

[21] Mudde RF, van der Akker HEA., 2D and 3D simulation of an internal airlift loop reactor on the basic of a two-fluid model. ChemEngSei 2001;56:6351-6358.

[22] Joshi JB. Computational flow modelling and design of bubble column reactors. Chem Eng Sei 2001;56:5893-5933.

[23] Sokolichin A, Eigenberger G, Lapin A. Simulation of buoyancy driven bubble flow established simplification and open question. AIChEJ 2004;50:24-45.

[24] Bauer M, Eigenberger G. A concept for multi-scale modeling of bubble columns and loop reactors. Chem Eng Sei 1999;54:5109-5117.

[25] Krishna R, van Baten JM, Urseanu MI, Ellenberger J. Design and scale up of a bubble column slurry reactor for Fisher-Tropsch synthesis, Chem Eng Sei 2001;56:537-545.

[26] Johansen ST, Boysan F. Fluid dynamics in bubble stirred ladles: Part II. mathematical modelling. Metal Trans B 1988;19B:755-764.

[27] Clift R, Grace JR, Weber ME. Bubbles, Drops andPartieles. Dover Publication, Inc. New York; 1978.

[28] Schiller L, Neumann A. A drag coefficient correlation. VDI-Zeits 1933;77:318-320.

[29] Delnoij E, Kuipers JAM, van Swaaij WPM. Dynamics simulation of dispersed gas-liquid two-phase flow using a discrete bubble model. ChemEngSei 1997;52:3759-3772.

[30] Ranade VV, van der Akker. A computational snapshot of gas-liquid flow in baffled stirred reactor. Chem Eng Sei 1994;49:5175-5192.

[31] Bhanu C, Mazumdar D. Numerical prediction of melting rates in gas bubble driven systems. Trans Indian Inst Met 1997;50:249-258.

[32] Prosperetti A, Tryggvason G. Computational methods for multiphase flow. Cambridge University Press. Cambridge; 2007.

[33] Kastanek F, Sharp DH. Chemical reactorsfor gas-liquid systems. E. Horwood. Chichester; 1993.