The effect of the draft tube geometry on mixing in a reactor with an internal circulation loop – A CFD simulation

Chemical Engineering and Processing 94 (2015) 29–34

Contents lists available at ScienceDirect

Chemical Engineering and Processing: Process Intensification journal homepage: www.elsevier.com/locate/cep

The effect of the draft tube geometry on mixing in a reactor with an internal circulation loop – A CFD simulation P. Lestinsky a, * , M. Vecer b , P. Vayrynen c, K. Wichterle b a b c

Brno University of Technology, Faculty of Mechanical Engineering, Technicka 2, Brno 61669, Czech Republic VSB-Technical University of Ostrava, Faculty of Metallurgy and Material Engineering, 17. listopadu 15, Ostrava 70833, Czech Republic Aalto University, School of Science and Technology, Department of Metallurgy, Vuorimiehentie 2, Espoo, Finland

A R T I C L E I N F O

A B S T R A C T

Article history: Received 18 September 2014 Received in revised form 2 March 2015 Accepted 7 March 2015 Available online 18 March 2015

The main topic of this paper is to describe the effect of geometrical parameters on mixing time in a reactor with an internal circulation loop. The design of the draft tube inside the reactor is an important geometric parameter and has a big influence on two phase hydrodynamics as well as on mass transfer in the reactor. In the first section, the validation of the selected mathematical model is carried out. The results of experimental measurements (liquid velocity and gas hold-up) obtained on the laboratory scale reactor are compared with the CFD simulations. The CFD simulation (bubbly flow and mass transfer models) was made using COMSOL Multiphysics 3.5a. The results of the numerical simulation are in good agreement with the experimental data. In the second section, the study of mixing in the reactor is described with the new geometry of the draft tube. The standard experimental techniques for testing mixing processes are quite problematic because common tracers (soluble salts) have significant influence on two phase hydrodynamics inside the reactor. Though, an alternative nontrivial method had to be used. The split of the draft tube into two or three section has a significant effect on mixing (mass transfer) in the reactor. ã 2015 Elsevier B.V. All rights reserved.

Keywords: Internal loop reactor Mixing CFD Bubbly flow

1. Introduction Reactors with an internal circulation loop are found in an increasing number of applications in the chemical industry, biochemical fermentation and biological wastewater treatment processes. These include fermentation [1]; production of microorganism cells [2]; flue gas desulphurization [3]; water and wastewater treatment [4]; the synthesis of methanol or dimethyl ether from syngas, Fisher–Tropsch synthesis and coal liquefaction [5]; and petroleum refining [6,7]. For the design of an airlift reactor with an internal circulation loop, it is necessary to have accurate estimates of the phase holdups and velocities in the riser and downcomer. Several literature studies have focused on the estimation of these hydrodynamics parameters [8,9]. The velocities of the liquid in the downcomer and in the riser are dependent on frictional losses which are determined by the geometry of the airlift reactor and operating condition. Good mixing is essential for some possesses to enhance the chemical reaction rate and decrease possible side reactions caused by non-uniform concentration or temperature profiles.

Mixing behaviour is also crucial for scaling up and optimizating reactors [10]. Several recent publications have established the potential of CFD simulation for describing the hydrodynamics of reactors with internal circulation loops [11–14]. 2. Numerical model The simulation of an internal loop airlift reactor was carried out using COMSOL under conditions similar with previous experiments. The governing equations calculated in this case were as follows. 2.1. Hydrodynamics The momentum equation for liquid phase:

el rl

ð@ul Þ þ el rl ul rul ¼ rp þ r½el ðhl  hT Þðrul þ rul T Þ @t þ el rl g  F;

(1)

where hT is turbulent viscosity defined as: 2

k

hT ¼ rl C m ; f * Corresponding author. Tel.: +42 541144955. E-mail address: [email protected] (P. Lestinsky). http://dx.doi.org/10.1016/j.cep.2015.03.009 0255-2701/ ã 2015 Elsevier B.V. All rights reserved.

(2)

30

P. Lestinsky et al. / Chemical Engineering and Processing 94 (2015) 29–34

Nomenclature

Symbols Downcomer cross-section area (m2) Ad Ag Cross-section area to gas input (m2) Ar Riser cross-section area (m2) C Dimensionless concentration (1) c1 Homogenous concentration at the mixing time cD Drag coefficient (1) Ck Bubble induced turbulence modelling constant, Ck = 0.01 – 1 (1) Cm Turbulence modelling constant, Cm = 0.09 (1) Cf Bubble induced turbulence modelling constant, Cf = 1 – 1.9 (1) Cf1 Turbulence modelling constant, Cf1 = 1.44 (1) Cf2 Turbulence modelling constant, Cf2 = 1.92 (1) db Bubble diameter (m) din Inner draft tube diameter (m) dout Outer draft tube diameter (m) F Volume force (interaction force) (N/m3) g Gravity acceleration, g = 9.81 (m2/s) hB Distance of the draft tube from the bottom (m) hT Distance of the draft tube from the liquid level (m) HRT The overall hydraulic retention time (s) k Kinetic energy (m2 /s2) KB,T Loss coefficient at the bottom and at the top (1) L Length of draft tube (m) M Molecular weight of gas, Mair = 0.029 (kg/mol) NOC Number of liquid circulation (1) Nr Gas mass flux (kg/(m2 s)) p Hydrostatic pressure (Pa) Qg Input gas volumetric flow (m3/s) QL Liquid volumetric flow in the reactor (m3/s) r Diameter of integrated volume R Universal gas constant, R = 8.314 (J/(mol K)) Sk Additional source term (kg/(m3 s)) T Temperature of liquid, T = 288 (K) tmix Mixing time (s) ug Gas velocity (m/s) Ug Superficial gas velocity (m/s) ul Liquid velocity (m/s) Ul Superficial liquid velocity (m/s) u Slip velocity (m/s) Vreactor Volume of liquid in the reactor (m3) Dhr,d Difference of level in reverse U-tube manometer in riser and downcomer (m) Dp Pressure difference (Pa) Dzr,d Difference of height between delivery point from reverse U-tube manometer in riser and downcomer (m) eg Gas hold-up (1) egr,gd Gas hold-up in riser and downcomer (1) el Liquid hold-up (1) hl Viscosity of liquid (Pa s) hT Turbulent viscosity (Pa s) j Loss coefficient (1) rDr Density of gas-liquid dispersion in riser (kg/m3) rg Density of gas (kg/m3) rl Density of liquid (kg/m3) sk Turbulence modelling constant, s k = 1 (1) sf Turbulence modelling constant, s f = 1.3 (1) f Dissipation rate of turbulent energy (m/s3)

where k is the turbulent kinetic energy, f is the dissipation rate of turbulent energy and their values are calculated by closure equations: 
  @k h 1 rl  r h þ T rk þ rl ul rk ¼ hT ðrul þ rul T Þ2  rl w þ Sk ; 2 @t sk (3)

rl








@e h  r h þ T r’ þ rl ul r’ s’ @t 1 ’ ’2 ’ ¼ C ’1 hT ðrul þ rul T Þ2  rl C ’2 þ C ’ Sk 2

k

k

k

(4)

where Sk is the additional source term comprehensive of turbulence induced by bubble: Sk ¼ C k eg rpuslip :

(5)

The continuity equation for two phase systems is:

@ ðe r þ eg rg Þ þ rðel rl ul þ eg rg ug Þ ¼ 0: @t l l

(6)

The momentum equation for gas phase: ð@eg rg Þ þ rðeg rg ug Þ ¼ 0: @t

(7)

Invoking the equation of state for each phase:

rl ¼ const:; rg ¼

ðMpÞ ðRTÞ

(8)

and following closure relation:

el þ eg ¼ 1:

(9)

We obtained a closed system of ten differential and algebraic equations which describe the dynamics of the two-phase flow and can be solved numerically when specific conditions are reached. Another option is to provide an initial calculation with a simplified model. When a low gas fraction occurs in the system, the value of el is close to one. If we neglect the variation of el due to the presence of gas in the continuity equation for the liquid phase Eq. (6) becomes:

@el þ rðel ul Þ ¼ 0: @t

(10)

It will be transformed into (with the assumption e  1): rul ¼ 0:

(11)

The gas velocity is then calculated as a sum of liquid velocity, slip velocity and drift velocity: ug ¼ ul þ uslip þ udrift ;

(12)

where the slip velocity uslip is calculated by the Eq. (17), and drift velocity udrift is the additional condition for an employed turbulent model, calculated as: udrift ¼

hT reg rl el

(13)

COMSOL Multiphysics 3.5a uses simplification, where the Euler– Euler two phase model is simplified to the mixed model, where the continuity equation is used for calculating gas hold-up. This solution does not need computational sources and calculation is possible on a PC. 2.2. Interfacial forces, slip velocity and drag coefficient The magnitude of F describes the interaction force between the continuous and the dispersed phase. COMSOL Multiphysics 3.5a makes it possible to input the interfacial force F as a constant. The

P. Lestinsky et al. / Chemical Engineering and Processing 94 (2015) 29–34

magnitude of the interfacial force has a great influence on the resulting liquid velocity and gas hold-up F¼

ðeg ðrg  rl Þgðug  ul ÞÞ : ðug  ul Þ

(14)

The Eq. (14) was used for our simulation because it does not contain a constant as a drag coefficient and bubble diameter. These constants obtain an equation of the slip velocity and we obtained two independent variables (slip velocity and interfacial force). In simulation, we used the simplification of continuity equation; therefore, we simplified Eq. (6) to the term: F ¼ ðrg  rl Þg:

(15)

The interfacial force after substituting constants is: F ¼ 9778N=m3 :

(16)

The next parameter is slip velocity, which is the relative velocity between the gas and liquid phase. It is one of the main parameters in multiphase flows. Slip velocity, together with interfacial force, belong to variables which could affect the result of numerical simulations. The slip velocity can be calculated as follows: uslip ¼ 


0:5 ð4db DpÞ : ð3cD rl Þ

(17)

The most important parameters here are the bubble diameter – db and the drag coefficient – cD. Another parameter is the hydrostatic pressure – p, but its values cannot be affected. The bubble diameter is assumed as a constant (equivalent sphere diameter) over all the reactor. This assumption has also been used in the paper written by Oey [12]. The drag coefficient for exact bubble diameter was obtained from previous experimental work, Vecer et al. [15]. 2.3. Convection and diffusion The diffusion and convection were calculated by:

@c þ rðDrc þ cul Þ ¼ 0; @t

(18)

where ul is the liquid velocity solved in the previous calculation (in hydrodynamics). Trace was injected using Eq. (19):
  2 1 (19)  3 Þ; wherec0 ¼ 1mol=l cin ¼ c0 expð t 2

31

A time dependent segregated solver was used for calculation. Time step was 0–100 s and the simulation was in a transient state. The segregated groups are: 1 – ul,vl (liquid velocity, radial and vertical component), p (pressure), rhog eff (effective gas density); 2 – log k (logarithm of turbulent kinetic energy), log d (logarithm of turbulent dissipation rate); 3 – c (concentration). 3. Experimental The experimental apparatus (a laboratory scale reactor with an internal circulation loop) consists of a glass vessel and has a diameter of 15 cm and height of 150 cm. A gas hold-up is measured by the pressure difference in a reverse U-tube manometer. A Pitot tube was used for liquid velocity measurements. A schematic view of the experimental apparatus is shown in Fig. 1. The internal circulation loop is created by the draft tube into the reactor. Draft tubes were used with a length of 0.8 m and outer diameters of 0.06 m, 0.08 m and 0.1 m. Distances of the draft tubes from the bottom were 0.02 m, 0.03 m, 0.035 m and the distance of the draft tube from the liquid level was 0.7 m. Experiments were carried out with a two-phase air-water system. Data were measured for input gas volumetric flows from 0.11 to 0.56 dm3/s. In order to check the CFD model, data were measured (liquid velocity and gas hold-up) for three types of draft tube geometry, see Table. 1. The main measured quantity was pressure difference resulting from the gas hold-up. The gas hold-up is defined as a volume fraction of gas phase in the gas-liquid dispersion in the reactor.

egr;d ¼

Dhr;d rl ðrl  rg Þ Dzr;d

(20)

The value of liquid circulation velocity, Ulr, is one of the most important design parameters for reactors with circulation loops. Liquid circulation has an influence on the gas hold-up in the vessel, on the heat and mass transfer coefficient and on the extent of mixing in the reactors. The liquid velocity was measured by Pitot tube. The pressure difference was read from the inverse U-tube manometer: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2g rl Dhr (21) U lr ¼

rDr

2.4. Mesh grid Triangular mesh elements were used for creating a computational domain. The distance between the elements was 0.007 m. The distance between the mesh elements at problematic places (both edges of draft tube) was 0.001 m. 2.5. Numerical stabilization Isotropic diffusion or 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. The gas transport equation is stabilized by default using both terms. It is worth highlighting that a suitable scale for the effective gas density and the bubble number density must be specified. More detailed information about numerical stabilization is in the manual COMSOL.

Fig. 1. Experimental apparatus; 1-laboratory airlift reactor, 2-draft tube, 3-gas inlet and flow-meters, 4-reverse U-tube manometer to measure pressure difference (to the interpretation of gas hold-up), 5-reverse U-tube manometer to measure liquid velocity by Pitot tube, 6-input tracer, 7-gas bubbles distributor, 8-escape valve.

32

P. Lestinsky et al. / Chemical Engineering and Processing 94 (2015) 29–34

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

Type 1 Type 2 Type 3

hB (m)

hT (m)

din (m)

dout (m)

Ad/Ar (m2/m2)

L (m)

0.02 0.03 0.035

0.07 0.07 0.07

0.052 0.072 0.092

0.06 0.08 0.1

7.05 3.07 1.49

0.8 0.8 0.8

3.1. Mixing time

used for the CFD simulation for contaminated systems related to inlet gas volumetric flow. The values of liquid velocity and gas hold-up were calculated as integrated variables per volume of the riser and downcomer: R ð 2prul dU l Þ R (22) ul ¼ ð dvÞ

eg ¼

The geometry of type 3 was used for the experiment with tracer. The mixing time was measured for the input gas volumetric flow 0.056 dm3/s. Potassium permanganate was used as a tracer to measure mixing inside the reactor and timed with a stopwatch. The process of mixing was also recorded using digital camera. The 5 ml tracer solution of potassium permanganate was injected into the bottom of the reactor (0.05 m from the gas distributor). Images were recorded continually during the circulation of the colour tracer inside the reactor. The mixing time was determined as the length of time it took for the colour of the tracer to be a uniformconstant value. The images were taken by a Nikon digital camera (model D5100) which has a resolution of 4928  3264 pixels (16 megapixels) and all the images were taken under the same conditions (same exposure). The RGB pictures were converted into Grey-scale. The mixing time line of potassium permanganate in water is shown in Fig. 2. Behin [16] used a reactor with external circulation and uniform illumination. The colour intensity of the tracer was measured at several locations in the reactor. In our case, the colour intensity of the tracer cannot be measured using image analysis because the internal circulation inhibits the colour intensity of the tracer in the reactor. The experimental result confirms the result of the CFD simulation. The effect of the draft tube on the mixing was carried out only within the CFD simulation. The treatment of the experimental apparatus would be difficult and time consuming. The CFD simulations of different geometries of the draft tube significantly reduce the time and to allow testing of different geometric variants, such as internal circulation reactor with the three draft tubes. 4. Results 4.1. Checking the CFD model using two-phase hydrodynamics The model was constructed in COMSOL Multiphysics 3.5a. Bubble diameters and drag coefficients from Vecer et al. [15] were

R ð 2preg deg Þ R ð dvÞ

(23)

The liquid velocity in the riser increases with the increase of gas volumetric flow. The influence of the draft tube diameter on the liquid velocity in riser can be seen in Fig. 3. The more the diameter of the draft tube decreases, the more the liquid velocity increases. A fairly good agreement of experiments and numerical simulations can be seen for the prediction of liquid circulation velocity with increasing gas inlet volumetric flow. Gas hold-up is one of the important parameters of airlift reactors and bubble columns. Gas hold-up in reactors increases the more gas volumetric flow increases. The linear trend of gas hold-up in the riser can be seen in Fig. 3. The difference between the experimental data and the CFD simulation was caused by the assumption of a homogeneous distribution of bubble diameter in the reactor and neglecting bubble coalescence. Both assumptions are containing in the bubbly flow model built in COMSOL. A propriety check of the CFD simulation was performed and the prediction of liquid velocity agreed well with the experiment. 4.2. Checking the CFD model on the mixing The model was constructed in COMSOL Multiphysics 3.5a. Bubble diameters (db = 5.8 mm), drag coefficients for contaminated systems (cD = 2.46) and input gas volumetric flow (Qg = 0.056 dm3/ s) were used for the CFD simulation. The values of liquid velocity were calculated as integrated variables per volume of riser: R ð 2prul dU l Þ R (23) ulr ¼ dV r The value of liquid volumetric flow in the reactor was calculated using the liquid velocity in the riser and cross-section area: Q L ¼ Ulr Ar

Fig. 2. Photos of potassium permanganate in water (tmix = 70, measured using a stopwatch).

(24)

P. Lestinsky et al. / Chemical Engineering and Processing 94 (2015) 29–34

33

Fig. 3. Liquid velocity and gas hold-up in the riser as a function of gas volumetric flow for various draft tube diameters.

In the CFD simulation, the concentration values were recorded for 0.1 s in different locations of the reactor (see Fig. 4). For a given time, the average value of concentration and the standard deviation were calculated. The variation coefficient was calculated as a ratio of the standard deviation and the average value. This variation coefficient was calculated for all time steps. The homogeneous distribution of the tracer concentration in the reactor was reached when the variation coefficient was less than, or equal to, 5%. At this moment, the time was designated as the mixing time (tmix). Then, the final concentration (c1) of the tracer in the reactor was achieved at the mixing time. For better graphical rendering, a dimensionless concentration was used. The dimensionless concentration (C) was calculated as: C¼

ðci  c1 Þ c1

(25)

The dependence of the dimensionless concentrations on the mixing time can be seen in Fig. 4. The number of liquid circulation (NOC) was used as the main indicator of the influence of the draft tube geometry on the mixing: NOC ¼

tmix HRT

(26)

Where the overall hydraulic retention time (HRT) was calculated as: HRT ¼

V reactor QL

(27)

Fig. 4. The change of tracer concentration at various places in the reactor over time.

The liquid circulation was faster (NOC = 6.5) in the CFD simulation than in the experiment (NOC = 4.9) Fig. 5. The difference was due to a 25% deviation of liquid velocity between the CFD simulation and the laboratory experiment, but the mixing time was similar. The calculation of the vortex behind the bubbles and the movement of the liquid layer were not calculated in the CFD simulation. 4.3. Checking the geometry of the draft tube by mixing: The difference between the liquid velocity measured and that calculated by the CFD simulation was higher at higher gas volumetric flows Table 2. Thus, for the detecting the effect of the draft tube geometry on the mixing, a lower gas volumetric flow of Qg = 0.000167 m3/s was chosen. The draft tube was divided into the two parts. The total length of the divided draft tube was kept to a constant (0.8 m) and the distance between the parts was exactly adjusted. Results obtained from the simulation for severe cases are listed in the Table 3. It is evident, that the overall hydraulic retention time increases with increasing distance between parts of draft tube, contrary to the mixing time and number of liquid circulations which are decreases. The CFD simulation divided the draft tube into three parts. Mixing with three parts of the draft tube is better than with two parts, see Table. 4 The concentration of the tracer at point G is in the same time as at point B, but the tracer is already at point B for the second time. When the tracer is at point G for the second time, then the homogeneous concentration in the reactor has nearly been reached. This can be seen in Fig. 6.

Fig. 5. The change of the tracer concentration at various places in the reactor over time (2  DF).

34

P. Lestinsky et al. / Chemical Engineering and Processing 94 (2015) 29–34

Table 2 Mixing in the reactor; comparison of the CFD simulation and laboratory experiment.

Vreactor(m3) QG(m3/s) QL(m3/s) HRT(s) tmix(s) NOC(-)

CFD simulation

Experiment

0.0158 0.000056 0.0015 10.3 67 6.5

0.0011 14.4 70 4.9

Table 3 The effect of the draft tube geometry (distance between draft tubes) on the mixing in the reactor.

concentration was monitored at several points in the reactor. At the moment when concentration at all monitored points was constant (5% variance), the mixing time was reached. The split of the draft tube onto two or three zones has a big effect on the mixing (mass transfer) in the reactor. The mixing time decreases with the increase of draft tube segments in the reactor. On the other hand, two phase hydrodynamics was not significantly affected. A relatively simple CFD model was used for description of a relatively complex hydrodynamic situation with sufficient accuracy. There are plenty of studies dealing with CFD modelling of airlift reactor by complicated models based for instance on turbulent viscosity or other relevant parameters [17,18]. Unfortunately an application of such approaches in industry conditions is still poor.

CFD simulation 0.0158 Vreactor(m3) QG(m3/s) 0.000167 Type of geometry 1 x DF 2  DF 25 mm 2  DF 2  DF 2  DF 50 mm 75 mm 100 mm 3 QL(m /s) 0.00245 0.00216 0.00202 0.00187 0.00177 HRT(s) 6.4 7.3 7.8 8.4 8.9 54 38 28 17 16 tmix(s) NOC(-) 8.4 5.2 3.9 2 1.8

Acknowledgements The authors gratefully acknowledge the financial support of the project “Excellent young researcher at BUT” No. CZ.1.07/2.3.00/ 30.0039 of Brno University of Technology, and project ICT-National Feasibility Study (No. LO1406). References

Table. 4 The effect of the draft tube geometry (distance between draft tubes) on the mixing in the reactor. CFD simulation Vreactor(m3) QG(m3/s) Type of geometry

0.0158 0.000167 1  DF

QL(m3/s) HRT(s) tmix(s) NOC ()

0.00245 6.4 54 8.4

2  DF 50 mm 0.00202 7.8 28 3.9

3  DF 50 mm 0.0025 6.3 18 2.9

Fig. 6. The change of the tracer concentration at various places in the reactor over time (3  DF)

5. Conclusion A CFD simulation was used for the study of mixing in a reactor with a draft tube. The concept of two or three zones of the draft tube was applied to a numerical investigation of the mixing processes. This approach was simpler than the complicated experimental study. One point of interest was the mixing time. The virtual tracer was injected into the reactor and its

[1] D.J. Pollard, A.P. Ison, P.A. Shamlou, A.E. Saez, Reactor heterogeneity with Saccharopolyspora erythraea airlift fermentations, Biotechnol. Bioeng. 58 (1998) 453–463. [2] S.W. Kim, S.W. Kang, J.S. Lee, Cellulase and xylanase production by Aspergillus Niger KKS in various bioreactors, Bioresour. Technol. 96 (1997) 63–67. [3] M.R. Jafari Nasr, H.R. Bakhtiary, A. Karimi, A. Tavassoli, Single step H2S removal using chelated iron solution investigation of hydrodynamic parameters in an internal loop airlift reactor, Iranian J. Sci. Technol. 28 (2004) 643–651. [4] B. Jin, X.Q. Yan, K. Yu, H.J. van Leeuwen, Comprehensive pilot plant system for fungal biomass protein production and wastewater reclamation, Adv. Environ. Res. 6 (2002) 179–189. [5] W.G. Yang, J.F. Wang, Y. Jin, Gas–liquid mass transfer in a slurry bubble column reactor under high temperature and high pressure, Chin. J. Chem. Eng. 9 (2001) 4–9. [6] M.R. Mehrnia, B. Bonakdarpour, T. Jafar, M.M. Akbarnegad, Design and operational aspects of airlift bioreactors for petroleum biodesulfurization, Environ. Prog. 23 (2004) 206–214. [7] F.P. Shariati, B. Bonakdarpour, M.R. Mehrnia, Hydrodynamics and oxygen transfer behavior of water in diesel microemulsions in a draft tube airlift bioreactor, Chem. Eng. Proc. 46 (2007) 334–342. [8] J.J. Heijnen, J. Hols, R.G.J.M. van der Lans, H.L.J.M. van Leeuwen, 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, Chem. Eng. Sci. 52 (1997) 2527–2540. [9] R.S. Oey, R.F. Mudde, L.M. Portela, Simulation of a slurry airlift using a two-fluid model, Chem. Eng. Sci. 56 (2001) 673–681. [10] J. Behin, Modelling of modified airlift loop reactor with a concentric doubledraft tube, Chem. Eng. Res. Design. 88 (2010) 919–927. [11] J.M. van Baten, J. Ellenberer, R. Krishna, Using CFD to describe the hydrodynamics of internal air-lift reactors, Canad. J. Chem. Eng. 81 (2003) 660–668. [12] R.S. Oey, R.F. Mudde, H.E.A. van der Akker, Numerical simulation of an oscillating internal-loop airlift reactor, Canad. J. Chem. Eng. 81 (2003) 684–691. [13] Q. Huang, Ch. Yang, G. Yu, Z.S. Mao, CFD simulation of hydrodynamics and mass transfer in an internal airlift loop reactor using a steady two-fluid model, Chem Eng. Sci. 65 (2010) 5527–5536. [14] M. Simcik, A. Mota, M.C. Ruzicka, A. Vicente, J. Teixeira, CFD simulation and experimental measurement of gas holdup and liquid interstitial velocity in internal loop airlift reactor, Chem. Eng. Sci. 66 (2011) 3268–3279. [15] M. Vecer, P. Lestinsky, K. Wichterle, M.C. Ruzicka, On bubble rising in countercurrent flow, Int. J. Chem. React. Eng. 10 (2012) . [16] J. Behin, N. Farhadian, Residence time distribution measurement in a two dimensional rectangular airlift reactor by digital image processing, Exp. Therm. Fluid Sci. 51 (2013) 244–250. [17] M. Simcik, M.C. Ruzicka, A. Mota, J. Teixeira, Smart RTD for multiphase flow systems, Chem. Eng. Res. Design. 90 (2012) 1739–1749. [18] J. Chang, F. Meng, L. Wang, K. Zhang, H. Chen, Y. Yang, CFD investigation of hydrodynamics, heat transfer and cracking reaction in a heavy oil riser with bottom airlift loop mixer, Chem. Eng. Sci. 78 (2012) 128–143.