Euler–Lagrange CFD Simulation of a Gas–Liquid Fluidized Bed

ARTICLE pubs.acs.org/IECR

Euler
Lagrange CFD Simulation of a Gas
Liquid Fluidized Bed Reactor for the Mineralization of High-Strength Phenolic Wastewaters Rodrigo J.G. Lopes,* M.L.N. Perdigoto, and Rosa M. Quinta-Ferreira Centro de Investigac-~ao em Engenharia dos Processos Químicos e Produtos da Floresta (CIEPQPF) GERSE, Group on Environmental, Reaction and Separation Engineering Department of Chemical Engineering, University of Coimbra Rua Sílvio Lima, Polo II, Pinhal de Marrocos, 3030-790 Coimbra, Portugal ABSTRACT: A state-of-the-art Euler
Lagrange model was developed to simulate the ozonation of phenol-like pollutants in a bubble column reactor. First, several numerical simulations were performed to evaluate on how the bubble velocity and oxidant concentration can improve the detoxification of liquid effluents by noncatalytic ozonation. Second, the effect of inlet ozone velocity as well as the influence of inlet ozone concentration has been investigated comparatively under different process conditions. We found that the Eulerian
Lagrangian computations have correctly handled the experimental observations in the quasi-homogeneous flow regime both in terms of the gas
liquid velocity distributions and normalized pollutant concentration. The numerical confidence exhibited by the CFD simulations underlined the ozonation-based technology as one promising application to improve the environmental performance of bubble columns, which is typically operated under the low-interaction regimes, especially when the mass transfer of ozone is rate controlling and affects the mineralization rate. The interstitial flow maps have been successfully correlated with total organic carbon concentration profiles as function of inlet ozone velocities and concentrations. Moreover, the multiphase CFD framework gathered positively the mixing degree induced by different inlet bubble velocities as demonstrated by the total organic carbon concentration mappings and experimental conversion data.

1. INTRODUCTION Among the extensive range of multiphase reactors, bubble columns are the most pervasive gas
liquid reactors in industry mainly because of their distinctive heat and mass transfer characteristics, plain construction, and low operating costs. Typical applications differ from petrochemical, pharmaceutical, and biochemical industries to water and gas pollution abatement technologies.1 Notwithstanding the voluminous literature devoted to their investigation and the complex bubble-induced hydrodynamics, the design and scale-up of bubble columns depends overly on semiempirical correlations. During the past decade, an integrated approach to investigate bubble column hydrodynamics was accomplished both experimentally and computationally for scale-up and design considerations. As the gas holdup, bubble size, and bubble rise velocity dictates predominantly the overall performance of bubble column reactors, there are also important factors that play a major role on the reaction yield, such as the bubble
bubble interactions and mixing rate.2,3 Here, advanced experimental characterization methodologies have been proposed based on magnetic resonance imaging and electrical capacitance/impedance tomography techniques, but since they are relatively expensive and cumbersome to operate, validated computational fluid dynamics (CFD) multiphase models appeared as a more cost-effective tool to gain further insights into the interstitial flow field and chemical reaction occurring in bubble column reactors. In the realm of multiphase reaction engineering, state-of-theart CFD codes have been used to perform numerical simulations of bubble columns either employing Eulerian
Eulerian models,4
7 Eulerian
Lagrangian models,8,9 or volume of fluid (VOF) models.10 Whereas the Eulerian
Eulerian model treats r 2011 American Chemical Society

gas bubbles and liquid phases as interpenetrating continua, it also describes the motion for gas and liquid phases in an Eulerian frame of reference. Using this mathematical concept for the dispersed and continuous phases, two-dimensional simulations of gas
liquid flows for a rectangular bubble column were investigated by Sokolichin and Eigenberger,7 Pan and Dudukovic,5 and Monahan et al.4 The Eulerian
Eulerian approach has been also applied to predict time-averaged axial liquid velocity revealing that the liquid primarily traveled up the column center and reduced gas holdup with liquid downflow at the walls. Similarly, the Eulerian
Lagrangian framework describes the liquid phase according to an Eulerian representation, while the gas is treated as discrete bubbles and each bubble is tracked by solving the Newtonian equations of motion for individual bubbles. This mathematical formulation is able to simulate detailed bubble
bubble interactions along with interfacial forces in the low interaction bubbly flow regime. The interstitial multiphase flow can be qualitatively identified in terms of bubble plume size and shape, time-dependent phase velocities and volume fractions, which allow ascertaining about the multiple-staggered vortex mode of liquid-phase circulations within the bubble column.11 The constitutive flow equations can be further solved by means of high spatial and temporal resolution methods, and Special Issue: CAMURE 8 and ISMR 7 Received: September 17, 2011 Accepted: November 21, 2011 Revised: November 17, 2011 Published: November 21, 2011 8891

dx.doi.org/10.1021/ie2021366 | Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research here, we have the multiphase VOF model to obtain the gas and liquid flow fields with better accuracy. The interstitial flow patterns are computed to provide time-dependent behavior of a dispersed bubbling flow and to account for the coupling effects of the pressure field and the liquid velocity on the bubble motion. A volume-tracking scheme is commonly used to follow the development of the gas
liquid interface.10 As long as state-of-the-art CFD codes make conceivable the interstitial flow simulation of multiphase physical phenomena, here we investigate the hydrodynamics of an ozonation bubble column to evaluate the reactor performance by means of an Euler
Lagrange framework that accounts for individual phase interactions. The theoretical framework accounts for uncoupled contributions arising from gravity, pressure, virtual mass, drag, lift forces, and chemical kinetics of ozonation is incorporated into the multiphase CFD model and the mineralization efficiencies are investigated comparatively under different operating conditions.

2. PREVIOUS WORK Svendsen et al.12 have reviewed the modeling of vertical bubble-driven flows. They focused on case-studies, where the bubble movement itself was the main source of momentum to the flow field and were often characterized by low superficial liquid velocities, relatively high superficial gas velocities. Both classical and current modeling approaches were presented in terms of gravity, buoyancy, centrifugal forces, conventional Magnus and Saffman forces, form and friction drag, and added mass, as well as turbulent migration and other instability mechanisms. The application of CFD to model gas-driven bubbly flows was evaluated for dynamic and steady-state descriptions and Euler
Lagrange versus Euler
Euler formulations. They have come to the conclusion that all the two-dimensional models, whether steady state or dynamic, lack the ability to predict the relatively slow rotational movement of bubble clusters observed experimentally in bubble columns. Three-dimensional models solved numerically on a very fine grid were proposed, which apart from adding strongly to the computer-resource demand, it also allows evaluating the cumbersome interstitial physical phenomena. More recently, the VOF approach has been used to investigate the passage of a gas bubble at a plan liquid
liquid interface.13 The numerical tracking methodology on the transient nature of the bubble characteristics gave interesting insight into the dynamic behavior of the interface. The authors have compared the VOF simulations with the experimental results showing the effect of the bubble size as well as the rheological properties of the heavy phase on the bubble’s retention time at the liquid
liquid interface. They found a qualitative agreement between the experimental data and the preliminary numerical results obtained by the VOF model. Aiming to derive new insights on flow maps for gas
liquid and liquid
solid risers with Lagrangian frameworks, Roy and Dudukovic14 have presented a fluid dynamics study on cold-flow circulating fluidized-bed riser using noninvasive methods. The time-averaged solids holdup distribution was determined by γray computed tomography, whereas the instantaneous solids and ensemble-averaged velocity patterns, as well as the solids residence time distribution in the riser, was obtained by computerautomated particle tracking. Accordingly, these experimental data were used for validation of a two-fluid Euler
Lagrange model, which is coupled with appropriate closures including the kinetic theory of granular solids. The authors have demonstrated

ARTICLE

the aptitude of the multiphase CFD model on the prediction of liquid and solids residence time distributions in the riser, as well as the solids velocity and holdup pattern. Moreover, the quantitative predictions from the two-dimensional axisymmetric simulation for the solids volume fraction and the time-averaged solids axial velocity radial distributions agreed well with the experimental data. Later, Diakov and Varma15 have reported an optimization analysis of packed-bed membrane reactor by modifying the oxygen feed distribution for the catalytic methanol oxidative dehydrogenation to formaldehyde. An Eulerian
Lagrangian was found to replicate identical distribution patterns and demonstrate the superiority of variable versus uniform permeability membranes. Higher productivities and increased reaction rates can be attained by means of the multiphase flow model that also demonstrated the most significant factor leading to the nonuniformity in optimal feed distribution of the reactants. An Eulerian
Lagrangian method has been also used to characterize the behavior of a heterogeneous cell population in a stirred-tank bioreactor with nonideal mixing.16 The authors have found that the modeling approach and the numerical method combined with a fractional-step method allowed for a stable, accurate, and numerically efficient solution of the underlying equations. The heterogeneity present in real reactors in both the fluid and cellular phases was accounted accurately and the theoretical approach was suggested to similar problems requiring a segregated and high-dimensional description of the internal structure of the dispersed phase. Following the flow regime diagnosis in bubble columns via pressure fluctuations and computer-assisted radioactive particle tracking measurements, Fraguío et al.17 presented a diagnostic analysis by noninvasive techniques to assist the validation of Eulerian
Lagrangian models. They proposed that experimental methodology to map the hydrodynamics in multiphase systems in a Lagrangian sense as long as the tracer particle trajectory time series was used to establish the optimal set of parameters. The time series led to successful flow regime identification by means of pressure fluctuation signals detected at various locations, which validated the method to detect and compare flow regimes to the ones computed by multiphase CFD codes. On the environmental and energy applications, Wang et al.18 have reported a numerical analysis of solid
liquid two-phase flow by means of an Eulerian
Lagrangian framework. The interfacial momentum transfer, buoyant force, and the collision force were considered, and the CFD model was successfully employed to determine the particle trajectories, the flow patterns, and the removal efficiency of the particles. A good agreement with the experimental results was claimed on the system productivity giving rise to the sensitivity of particle behavior to the collision force. Finally, and experimental and computational investigation of the hydrodynamics induced by a bubble plume in two-phase gas
liquid flows have been carried out by Ali and Pushpavanam.19 As the gas
liquid flow in such systems is inherently unsteady, particle image velocimetry (PIV) was used to experimentally determine the transient velocity fields in the multiphase system. Both the fluctuating and mean liquid velocities were compared effectively with a two fluid Eulerian
Lagrangian model. They reported that the numerical predictions from the 3D simulations captured the oscillating behavior found using 3D PIV, whereas the 2D simulations predicted a significantly higher value of turbulent viscosity, which emphasized the requirement of three-dimensional simulations to capture the oscillating behavior. 8892

dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research

ARTICLE

From the above literature survey, none of the CFD models published so far were used to simulate the ozonation of highstrength wastewaters. This work is devoted to the CFD simulations and experimental validation to verify the theoretical predictions on the mineralization of organic matter at different operating conditions. Here, we will evaluate several experiments at bubbly flow conditions querying the influence of superficial gas velocities and inlet ozone concentrations. Hence, an Euler
Lagrange CFD model coupling the transport processes such as mass transfer and ozonation kinetics is first developed accounting for the dependency of flow regimes and extent of interaction among gas and liquid phases. The aim of this work is to evaluate the ozonation performance of bubble column reactor without lost the discrete character of the underlying process as in the Eulerian
Eulerian approach. The remainder of this paper is outlined as follows. The multiphase CFD framework is presented with the constitutive equations and the boundary conditions. The experimental procedure is described after the simulation setup. The succeeding section shows the optimization of the numerical solution parameters and the results are organized by qualitative and quantitative comparisons of computed and experimental gas and liquid superficial velocities. The morphological features of interstitial flow maps for hydrodynamic and reaction variables are correlated under distinct process conditions.

FG ¼ Fb Vb g FP ¼
Vb ∇P 1 FD ¼
CD FL πRb 2 jub
uL jðub
uL Þ 2 FL ¼
CL FL Vb ðub
uL Þ  ∇  uL
FVM ¼
CVM FL Vb

Fb

dVb ¼ ðm_ L, b
m_ b, L Þ dt

dub ¼ Fb Vb dt

∑


 dVb F
Fb ub dt

ð1Þ

∂ ðαL FL Þ þ ∇ 3 αL FL u ¼ m_ b, L
m_ L, b ∂t

∑ F ¼ FG þ FP þ FD þ FL þ FVM

ð4Þ

∂ ðαL FL uÞ þ ∇ 3 αL FL uu ¼
αL ∇P
∇ 3 αL τL ∂t þ αL FL g þ Sm

ð5Þ

where the stress tensor τL is derived assuming the Newtonian dynamics for the liquid phase flow accounting for the effective viscosity μeff,L composed by the molecular viscosity and the turbulent viscosity.
 2 T ð6Þ τL ¼
μeff , L ∇u þ ∇u
I∇ 3 u 3 μeff , L ¼ μlam, L þ μturb, L In our case-study, the multiphase k
ε model has been used to compute turbulence properties in the range of the ozonation operating conditions T = 20 °C, P = 1 atm, ub e 0.4 cm/s. 3.2. Species Continuity Equations. The species continuity equation is expressed in eq 7 where the concentration of chemical species i in the liquid phase is accounted by CL,i and Si is the source term describing the production or consumption of species i because of homogeneous chemical reaction.

 μeff , L ∂ ðαL FL CL, i Þ þ ∇ αL FL uL CL, i
∇CL, i ∂t Sci ¼ m_ b, L Cb, i
m_ L, b CL, i þ αL Si

ð2Þ

and the specific forces for an individual bubble are computed by separate contributions of gravity, pressure, drag, lift, and virtual mass terms as shown in eq 3.



State of art closure models have been used to account for the net force acting on a bubble without conferring additional accuracy on the ozonation of organic compounds20 so we have calculated the closure terms as CD = 2/3(E€o)1/2 and CL = CVM = 0.5. The multiphase formulation of the volume-averaged Navier
Stokes equation was developed for the evaluation of liquid phase hydrodynamics. The bubbly phase is expressed by the liquid phase volume fraction αL. The source terms account for the interphase mass transfer m_ and the total interphase momentum transfer due to forces, Sm. The continuity and momentum conservation laws are represented by eqs 4 and 5, respectively.

3. COMPUTATIONAL FLOW MODEL 3.1. Governing Flow Equations. In the Euler
Lagrange model, we treat the mixture of the liquid and bubbles as a continuous phase whose motion is simulated using the Eulerian description. This multiphase framework is applicable to quasihomogeneous flow regimes characterized by moderately small spherical bubbles that do not coalesce nor break-up with low dispersed-phase velocities. The gas
liquid mass transfer is computed for each bubble following the surface renewal theory which comprises both physical and chemically enhanced mass transfer. The multidimensional distributions of chemical species in the interstitial space are derived from the coupled species conservation equations implemented in the Eulerian methodology. Subsequently, the numerical details of individual bubble
bubble interaction and the two-way coupling between phases will be presented along with ozonation chemical kinetics. The following governing equations are used for the numerical simulation where the motion of for each individual bubble is computed from the bubble mass and momentum equations including the interphase mass transfer rate m_ and the net force ΣF, respectively. The liquid contributions for an incompressible bubble are then expressed by eqs 1 and 2.

Dub DuL
Dt Dt

ð7Þ

The species mass concentration is computed from the overall species balance using eq 8 as we need to solve ns
1 transport equations for a mixture involving ns chemical species. ns

∑ Ci ¼ 1 i¼1

ð3Þ 8893

ð8Þ dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research

ARTICLE

Additionally, the liquid density and viscosity have been calculated as the average of properties of each species as shown in eq 9. 1 ¼ FL

ns

ns

∑ L, i ; μlam, L ¼ i∑¼ 1 CL, i μL, i i ¼ 1 FL, i C

ð9Þ

3.3. Interphase Mass Transfer. The mass concentration of chemical species i in the liquid and bubble phases are expressed by CL,i and Cb,i, respectively. At the bubble-liquid interface, both quantities are given by C/L,i and C/b,i. Taking into account a bubble with radius Rb, the interphase mass transfer is described as shown in eq 10.

m_ b, i ¼ EkL, i Ab FL ðCL, i
CL, i Þ

kL dB ¼ 2 þ 0:642½ReSc0:5 D O3

ð11Þ

The enhancement factor was calculated using the relation given by Westerterp 8 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 > Ha4 E∞ Ha2 <
Ha þ þ 1 E∞ > 1 2 þ 2ðE∞
1Þ E∞
1 E¼ 4ðE∞
1Þ > :1 E∞ e 1

ð12Þ

E∞ ¼

CTOC DTOC 1 þ 2DO3 HCO3

!rffiffiffiffiffiffiffiffiffiffiffi DO 3 DTOC

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kDO3 CTOC Ha ¼ kL

ð13Þ

ð14Þ

The Henry constant is then applied to calculate the mass concentration on the liquid interface side as long as the mass transfer resistance lies in the liquid phase. The total mass transfer is computed as the sum of the mass transfer rates for each single specie as expressed in eq 15. CL, i ¼ Hi Cb, i

Fb ; m_ b ¼ FL

ns

∑ m_ b, i i¼1

ð15Þ

3.4. Ozonation Kinetics. The ozonation kinetics of phenollike compounds have been derived considering a three-step mechanism with the compounds lumped into three groups: easier degraded pollutants (A), intermediates with difficult degradation (B), and the desired products (C).22 The reaction rates are given by eqs 13 and 14.


rTOCA ¼
dCTOCA =dt ¼ ðk01 þ k02 ÞCTOCA

ð16Þ


rTOCB ¼
dCTOCB =dt ¼ k03 CTOCB
k02 CTOCA

ð17Þ

CgO3/gO3/Nm3

k01 (min
1)

k02 (min
1)

k03 (min
1)

20

0.0048

0.011

1.47  10
10

40 88

0.0080 0.0143

0.022 0.180

0.0038 0.0053

After numerical integration, the normalized total organic carbon concentration is expressed as shown in eq 18. k01
k03 CTOC k02 0 0
k03 t ¼ 0 e þ e
ðk1 þ k2 Þt 0 0 0 0 0 CTOC0 k1 þ k2
k3 k1 þ k2
k3

ð10Þ

where E is the enhancement factor due to chemical reactions, Ab is the surface area of the bubble, and kL,i is the mass transfer coefficient for species i. The Sherwood correlation has been used to compute the latter parameter considering eq 11.21 Sh ¼

Table 1. Ozonation Kinetic Parameters for Different Oxidant Concentrations

ð18Þ This model encompasses two different routes for reactivity of oxidizable compounds: a direct conversion to end-products (k10 ) and a final oxidation preceded by formation of intermediates (k02, k03), and has been calibrated with the experimental conversion data obtained in the range of inlet ozone gas concentration used in our case-study. Table 1 summarizes the kinetic parameters of the apparent reaction rate dependent on ozone [k0 = f(CO3)] for each one of the three reaction steps involved in the use of different ozone inlet gas concentrations.

4. NUMERICAL SIMULATION AND BOUNDARY CONDITIONS 4.1. Multiphase Flow Fields and Interphase Coupling. Having described in detail how the Euler
Lagrange model has been developed to describe the ozonation bubbly flow, here the details of the numerical implementation and boundary conditions are presented accordingly. On the characteristic time scales, different time stepping methodologies were investigated to deal with the time-dependent motion of the bubbles and the liquid phase as well as the phenol-like compounds. The macroscopic liquid flow field is obtained through the numerical solution of the multiphase Navier
Stokes equations with longer time steps while the computation of interphase mass and momentum transfer is based on the bubble time step. In addition, shorter time steps are required for bubble
bubble collisions and for the sake of ozonation kinetics which encompasses a relatively fast chemical reaction. The binary bubble interactions are tracked within the computational domain after the establishment of the volume fraction of a single bubble included in a given node. The mass and momentum balance equations for each bubble are solved in a segregated fashion where the mass transfer rate is calculated explicitly following the equations already presented in section 3.3. This method evaluates the bubble volume rate of change given the numerical integration of the volume fraction of the bubbly phase. The bubble volume rate of change is iteratively computed to produce a new bubble size and consequently the bubble momentum equation is solved after the explicit estimation of the interphase momentum and mass transfer terms. Equation 2 is used to obtain the acceleration of each individual bubble and the current bubble velocity is achieved through the integration of Newton’s law. In a hard-sphere system the trajectories of the bubbles are determined by momentumconserving binary collisions. The bubble interactions are considered to be pairwise additive and instantaneous. In our numerical Euler
Lagrange simulations, the collisions were processed one 8894

dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research by one according to the order in which the events occur. Bearing in mind that for not too dense systems, the hard-sphere models are considerably faster than the soft-sphere models, the hard sphere model has been used following the work of Hoomans et al. where a constant time step is applied to account for the net force in order to deal with the binary collisions between bubbles.23 The velocity of bubbles has been assumed to change only due to binary collisions between individual bubbles within this time scale. The collision time is derived using the so-called bubble neighbor region that defines only the encapsulated bubbles for the collisions event. The liquid flow field is obtained from the discretization of the governing equations. The volume-averaged Navier
Stokes equations have been solved with a semi-implicit method for pressure linked equations (SIMPLE). The temporal derivatives have been integrated by first- or second-order explicit schemes and the convective terms were discretized either by first-, second-, or third-order schemes (MUSCL) as described by Lopes and Quinta-Ferreira24 for trickling flow conditions. The pressure gradient was computed implicitly and the interphase mass and momentum transfer was calculated explicitly. Apart from using a fully implicit algorithm for the discretization of the species continuity equation, the interphase mass transfer and ozonation reaction rate source terms were resolved explicitly. Whereas the mass and momentum balance equations for each bubble were solved sequentially, the resulting discretized species continuity equation was computed concurrently and subjected to the overall species balance expressed in eq 8. For a general time integrable variable j the formula can be written as jt+1 = jt + (dj/dt)tδtb. Using the numerical scheme described in the previous equation, the bubble volume rate of change is integrated to obtain the new bubble size. The acceleration of each individual bubble is then obtained in a straightforward manner using eq 2 and the bubble velocity is calculated by integrating the acceleration using this equation. The interphase coupling between the bubbles and the liquid phase is accomplished through the volume fraction αL, the interphase momentum transfer rate Sm, as well as the mass transfer rate from and to the bubbles. A mapping technique has been used to encompass the Lagrangian and Eulerian reference frames for the bubbles and liquid phase, respectively. Conversely to the literature works that presented translation techniques for exclusive bubble sizes smaller or larger than the computational grid aperture, we have implemented a discrete bubble model addressing deviations in bubble size because of coalescence and mass transfer processes. The mixed formulation allows mapping the Lagrangian variables to the Eulerian mesh and the Eulerian parameters to the Lagrangian framework following the work performed by Tomiyama et al.25 4.2. Boundary Conditions. A bubble column with similar characteristic dimensions to the geometry investigated by Pfleger et al. has been set up in our laboratory.26 Ozone with different superficial gas velocities in the range 0.1
0.4 cm/s was fed into the reactor by means of a perforate plate. The sparger configuration has 60 holes with a diameter of 1 mm being located in the center of the plate at square pitch of 5 mm and it is located at the bottom of the column. The nozzles are arranged in six equal groups of 10 needles. Various mesh densities were generated and different time steps were used in preliminary optimization studies to produce a parameter-independent solution. In addition, the Eulerian
Lagrangian model replicates the dimensions of the physical domain as follows 0.1 m  0.1 m  1 m.

ARTICLE

Appropriate boundary conditions were imposed according to the corresponding cell types. For interior domain cells and corner cells, no boundary conditions were specified and at impermeable walls no slip boundary-type condition was used. At the top a free slip boundary has been stipulated as well as for prescribed pressure cells. This methodology is required to compensate for the change of liquid volume due to bubbles entering and leaving the column so they were considered as inlet and outlet channels. Additionally, a Neumann-type boundary condition was used in all of the boundary cells for the species transport equation. The sparger holes have been simulated in such a way that ozone bubbles with specific size enter the column with a fixed velocity. 3-mm bubbles were considered for the ozonation of phenolic wastewaters and the typical distance between contiguous injections was specified to 3-fold of the bubble radius in order to avoid redundant binary collisions.

5. EXPERIMENTAL SECTION The synthetic phenolic wastewater was prepared using 100 ppm of each of the six phenolic acids (obtained from SigmaAldrich) corresponding to TOC0 = 370 mg of C/L, TPh0 = 350 mg of gallic acid/L, and COD0 = 970 mg of O2/L. The reaction experiments have been performed in the quasi-homogeneous flow regime and the gas oxidant was continuously fed to the reaction medium. Ozone was generated in situ by means of pure oxygen (99.999%, Praxair, Porto, Portugal) in an ozone generator (802N, BMT, Berlin, Germany) that operates based on the corona effect. The ozone gas concentration was monitored by means of an ozone gas analyzer (BMT 963 vent, BMT, Berlin, Germany). The ozone flow rate was typically 500 cm3/min with an inlet ozone concentration of 40 g/Nm3. The injection of ozone bubbles underneath the bottom plate was accomplished in such way that none of the bubbles enters the column at the same time. This method has been used to avoid pressure fluctuations at the top of the bubble column reactor as well as the artificial pulsing flow generated by incoming bubbles, conversely to the simultaneous introduction of ozone bubbles in all sparger holes. Liquid samples were withdrawn during the reaction and further analyzed in term of total organic carbon concentration. Total organic carbon (TOC) was measured with a Shimadzu 5000 Analyzer based on the combustion/nondispersive infrared gas analysis method. Each sample was run in triplicate in order to minimize the experimental errors. The deviation between the same sample runs was always lower than 2% for TOC, respectively. High-performance liquid chromatography was used to measure the concentrations of the individual compounds of the model effluent, as well as some of the reaction intermediates. The samples were injected via autosampler (Knauer Smartline Autosampler 3800), the mobile phase (20% of methanol in water slightly acidified) was pumped by a Knauer WellChrom K-1001 pump at a flow rate of 1 mL/min through an Eurokat H column at 85 °C, and detection was performed at 210 nm. 6. RESULTS AND DISCUSSION 6.1. Optimization of Mesh Density and Time Step. As long as the Eulerian-Lagrangian framework is able to predict each individual bubble trajectory, allowing consideration of additional effects related to bubble
bubble and bubble-liquid interactions, the geometric and numerical parameters have to be previously optimized to give independent results. We begin to address the 8895

dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research

ARTICLE

Figure 1. Grid sensitivity optimization of Euler
Lagrange predictions on the time-averaged bubble velocity profile for different mesh resolutions (uG = 0.3 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm).

Figure 2. Grid sensitivity optimization of Euler
Lagrange predictions on the time-averaged liquid velocity profile for different mesh resolutions (uG = 0.3 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm).

computational grid aperture and the time stepping strategy to advance the numerical computations on mass transfer with and without chemical reaction, bubble coalescence, and redispersion. The multiphase flow transport equations for each cell of the system have been solved by numerical schemes that enable the estimation of the scalar values at the faces of these cells from their centers values. First-order upwind schemes were initially applied to the upstream cell central values regarding the corresponding face expressed in accordance with the velocity gradients. This numerical scheme ensured fast convergence, but simulation results were not accurate and reliable. Subsequently, the secondorder upwind scheme has been used to estimate the value at a given face by combining the value of the upstream cell with the scalar gradient which was calculated by means of the GreenGauss cell based method of that same upstream cell to achieve second-order accuracy. The computed results agreed reasonably well with experimental data but so far the better accuracy was obtained with the third-order Monotone Upstream-centered Schemes for Conservation Laws scheme that is a total variation diminishing scheme that defines the scalar value at a given face by combining second-order upwind method with central differencing scheme. Consequently, the remainder of the CFD simulations was carried out with third-order accuracy which is wellknown recommended for the study of multiphase flow regimes.24 Henceforward, the Eulerian-Lagrangian approach has been used to evaluate the dynamics of gas
liquid flows with different grid refinements. The computational mesh was refined in eight levels: four levels for the coarser meshes with 3  105, 3.4  105, 5.6  105, 7.8  105 cells and four levels for the finer meshes with 2.1  106, 2.3  106, 2.7  106 and 3  106 cells. Figure 1 displays the effect of the number of computational cells on the time-averaged bubble velocity profile when the ozonation of phenolic wastewaters was simulated at uG = 0.3 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm. As can be seen from Figure 1, the increase of grid cells from 3  105, 3.4  105, 5.6  105, 7.8  105 to 2.1  106, 2.3  106, 2.7  106, 3  106 of computational nodes led to a better agreement between the computed results and experimental bubble velocities. The bubble velocity was roughly overpredicted in the center region of the bubble column reactor. Additionally, when mass transfer and ozonation kinetics

was simultaneously taken into account for both CFD calculations and experimental data, the averaged bubble velocity roughly decreased if compared to nonreacting flow conditions. A further comparison between computed results and experimental bubble velocities revealed that the time-averaged gas velocity profile predicted by the Euler
Lagrange model was somewhat higher in the center region of the bubble column reactor. The effect of grid refinement on time-averaged liquid velocity profile is shown in Figure 2 at uG = 0.3 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm. The time-averaged liquid velocity profile exhibited similar trends as the time-averaged bubble velocity counterpart. Indeed, the maximum upflow liquid velocity was found in the center of the ozonation column while downflow velocities are observed close to the wall region. Moreover, the ozonation reaction of phenolic wastewaters has the outcome to reduce the time-averaged liquid velocities and the assessment between the computed liquid and ozone velocity profiles unveiled a slip velocity almost crosswise invariant around the mean value. After reaching such an asymptotic solution for grid-size independence, the discrete bubble
bubble collision model which accounts for momentum exchange due to bubble collisions through four-way coupling can be used to evaluate the effect of lift force on the instantaneous and time-averaged flow properties. In this ambit, the enclosure of lift force with a lift coefficient of 0.5 led to higher lateral dispersion levels of bubbles in the upper part of the column. This fact translated in the breakdown of the recirculatory flow in the uppermost part of the bubble column reactor due to the lateral dispersion. Having mentioned the imperative validation for the dynamic characteristics under reacting flow conditions, we should accordingly perform unsteady simulations to identify how the timeaveraged properties such as liquid and bubble velocity agree with the experimental data. The noncatalytic ozonation of phenolic compounds were simulated at unsteady conditions to confer a realistic way of gas
liquid flows in environmental applications. The progress development of the Euler
Lagrange model was subsequently extended to evaluate the time-averaged gas and liquid velocities and phase holdup together with dynamic reaction conditions that have been compared with the experimental data. The influence of different time steps on time-averaged bubble velocity profile is shown in Figure 3 at uG = 0.3 cm/s, 8896

dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research

Figure 3. Parametric evaluation of time stepping methodology on the time-averaged bubble velocity profile (uG = 0.3 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm).

Figure 4. Parametric evaluation of time stepping methodology on the time-averaged liquid velocity profile (uG = 0.3 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm).

[O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm. When ozone is introduced into the bubble column filled with phenolic wastewater, the ozone bubbles formed at the sparger holes rise upward through the interstitial liquid exhibiting distinct length and time scales. The bubbly flow is intrinsically unsteady and the individual bubbles depending on the spatial coordinate make gross recirculatory convection streams. As can be seen form Figure 3, the predictions made by the multiphase CFD framework have shown that the meandering motion of bubbles was satisfactorily captured by the Euler
Lagrange model. Likewise, Figure 4 shows the effect of the time stepping strategy on the time-averaged liquid velocity profile when the ozonation of phenolic wastewaters was simulated at uG = 0.3 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm. Bearing in mind the time step range 10
3, 5  10
4, 10
4 and 10
5s, the optimized value between computational expenditure and numerical accuracy was accomplished with the shortest time step. 6.2. Gas
Liquid Flow Patterns. Aiming to evaluate the effect of the superficial gas velocity on the detoxification of organic

ARTICLE

Figure 5. Euler
Lagrange predictions on the temporal evolution of normalized total organic carbon converstion at the head of bubble column under different inlet bubble velocities ([O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm).

matter, three-dimensional computations have been carried out in terms of the total organic carbon conversion as a function of reaction time at different inlet ozone velocities. Figure 5 displays the temporal evolution of normalized total organic carbon concentration at the head of bubble column under different inlet bubble velocities when the ozonation of phenolic wastewaters was simulated at [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm. As can be seen, the total organic carbon conversion was 18.7, 22.6, 25.1, and 27.4% at 0.15, 0.25, 0.35, and 0.45 cm/s after 30 min and became significantly higher after 1 h of reaction time to reach 22.4, 26.3, 28.7, and 31.2% for 0.15, 0.25, 0.35, and 0.45 cm/s, respectively. Several gas
liquid two-phase simulations were performed to investigate the ability of the Euler
Lagrange model to capture the low frequency oscillations and its sensitivity to superficial gas velocity. We found that an average bubble size between 3 and 5 mm was observed for superficial gas velocities less than 0.05 cm/s. A distribution of bubble size varying from 1 mm to 20 mm was obtained if one further increases the ozone throughput which allowed inferring on how the multiphase CFD framework dealt with bubble size distribution. Both bubble coalescence and break-up processes contribute to the overall mineralization efficiencies so that methodology allowed us to understand the role of convection on the evolution of bubble size distribution within the bubble column reactor under reacting flow conditions. Figure 6 portrays typical snapshots mapping the hydrodynamic flow patterns resulting from CFD simulations for the ozonation of phenolic wastewater systematized spatially by the bubble size distribution, gas velocity, liquid velocity and ozone concentration field at uG = 0.15 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm. Here, a noteworthy ozone concentration field was typically found in the vicinity of the bubbles. As long as the ozonation reaction products were mainly transported by the convective streams originated from the multiphase flow field, the intermediate species and end products of ozonation are inherently conveyed with a similar liquid velocity field, thereby promoting the mineralization of phenolic compounds as shown by the temporal evolution of the normalized total organic carbon concentration, see Figure 5. According to velocity 8897

dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research

Figure 6. CFD snapshots mapping the hydrodynamic flow patterns resulting from Euler
Lagrange simulations organized spatially by the bubble size distribution, gas velocity, liquid velocity, and ozone concentration field (uG = 0.15 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm).

Figure 7. CFD snapshots mapping the hydrodynamic flow patterns resulting from Euler
Lagrange simulations organized spatially by the bubble size distribution, gas velocity, liquid velocity, and ozone concentration field (uG = 0.35 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm).

field maps, an improved uniformity of gas holdup distribution with slight local flow recirculations was found at the lowest superficial gas velocity (uG = 0.15 cm/s), which is in accordance with the experimental evidence. Conversely, if the superficial

ARTICLE

Figure 8. CFD snapshots mapping the hydrodynamic flow patterns resulting from Euler
Lagrange simulations organized spatially by the bubble size distribution, gas velocity, liquid velocity and ozone concentration field (uG = 0.45 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm).

ozone velocity increased up to uG = 0.35 cm/s at [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm as shown in Figure 7, the local recirculatory liquid flow and meandering bubble plume-like structures became more intense. At this stage, the Eulerian
Lagrangian model has been quantitatively reproduced the gas
liquid dynamics as revealed by the subsequent and good agreement between the computed results and experimental total organic carbon conversions. In fact, the three-dimensional transient simulations have been extended to evaluate higher superficial ozone velocities. The instantaneous snapshots of bubble size distribution, gas velocity, liquid velocity, and ozone concentration field taken at uG = 0.45 cm/s, [O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm were illustrated in Figure 8. As can be seen, the mixing degree between the gas and liquid phases was reinforced by stronger phase interactions so there was a more pronounced effect of lower gas volume fractions near the reactor wall and liquid ascending in the center as well as higher gas volume fractions in the center and resilient recirculating flow downward alongside the wall. In addition to the numerical predictions of time-averaged axial liquid velocity at different superficial gas velocities, one should also evaluate the gas holdup profiles that were obtained by time averaging over approximately 10, 250, 500, 1000, and 2000 s of real time. The predicted time averaged gas hold-up profiles at different liquid heights and superficial gas velocity are shown in Figure 9 and 10. Preliminary simulations were carried out to study the effect of grid resolution on predicted characteristic time scale of low frequency oscillations and time-averaged flow variables. The simulations were carried out using coarse and fine density grids. In order to obtain time-averaged profiles, the averaging process was carried out for different time intervals of real time, thereby illustrating on how computationally attractive can be accomplished for the high-density grid simulations. Here, we found that 8898

dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research

Figure 9. Euler
Lagrange predictions on the time-averaged gas holdup at the head of bubble column under different inlet bubble velocities ([O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm).

Figure 10. Euler
Lagrange predictions on the time-averaged gas holdup at the half-height of bubble column under different inlet bubble velocities ([O3] = 40 g/Nm3, T = 20 °C, and P = 1 atm).

the time-averaged gas holdup profiles obtained using the optimized grids resemble to the gas holdup measurements presented in the literature. Mainly due to increased computational requirements to carry out several cycles of averaging, Eulerian-Lagrangian computations with further grid refinement were quantitatively addressed for both hydrodynamic variables. The behavior of the gas holdup profiles is in accordance with moderate density dispersed gas
liquid multiphase systems as it simulated accurately the key dynamic and time-averaged flow properties satisfactorily over the range of system parameters studied for the ozonation of phenolic compounds. The radial profiles of time-averaged liquid and gas velocity and gas holdup were then used to evaluate the drag coefficient to provide a rigorous validation of the underlying model. As portrayed by Figures 9 and 10, the radial distribution of gas holdup exhibited the characteristic time scales of low frequency oscillations and the time-averaged flow properties have become insensitive to column height/diameter ratio, thereby reinforcing the experimental findings for bubbly flows. 6.3. Hydrodynamic and Mineralization Parameters. Apart from the hydrodynamic integral variables, the mineralization

ARTICLE

Figure 11. Effect of inlet ozone concentration on the mineralization efficiency predicted by the Euler
Lagrange model at the head of bubble column under different inlet bubble velocities (T = 20 °C and P = 1 atm).

parameters play a major role on the feasibility of the detoxification process. In this regard, the effect of inlet ozone concentration on the mineralization efficiency was investigated in the range of 20
95 g/Nm3. The transient profiles of the total organic carbon conversion are shown in Figure 11 at the head of bubble column under different inlet ozone concentrations when the ozonation of phenolic wastewaters was simulated at T = 20 °C and P = 1 atm. As can be seen, the degradation of organic matter followed a similar profile as the one we obtained for different inlet gas velocities. In fact, the total organic carbon conversion was 20.1, 25.2, 26.9, and 27.4% at [O3] = 20, 40, 88, and 95 g/Nm3, after 30 min and increased after 1 h of reaction time to 25.3, 28.7, 30.3, and 31.2% for [O3] = 20, 40, 88, and 95 g/Nm3, respectively. The unsteady development of ozone concentration has been spatially averaged over a grid scale which is roughly coarse than the individual bubble size. On the basis of the chemical oxidation mechanism, phenol-like pollutants were found to be quickly oxidized by ozone as this is extremely reactive with compounds comprising high electronic density sites. The reaction started rapidly and developed a steep total organic carbon concentration profile in the first fifteen minutes as depicted by Figure 11. Nevertheless, the organic intermediates were only partially oxidized without significant mineralization mainly after 30 min. As the reaction time evolved, the refractory compounds were slowly oxidized revealing lower detoxification rates. Indeed, the increase of ozone from [O3] = 20 to 40 g/Nm3 has not the counterpart effect for the first fifteen minutes of reaction comparing to the [O3] = 88 and 95 g/Nm3. The bubble column reactor had considerable different oxidation levels so the higher the ozone concentration fed into the reactor, the lower normalized total organic carbon concentrations at the end of the experimental run. The reaction system revealed to be considerably dependent on the ozone load, so the overall performance of the bubble column reactor can be improved for such a dynamic process behavior. 6.4. Interstitial Flow Mappings. To gain further insights on how the detoxification reaction is correlated with hydrodynamic variables, the Eulerian
Lagrangian model has been used to evaluate the interstitial reactive flow patterns. In fact, the ozonation of liquid pollutants is directly related with the strength of 8899

dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research

ARTICLE

Figure 12. Effect of inlet ozone concentration on the mineralization efficiency predicted by the Euler
Lagrange model at the half-height of bubble column (T = 20 °C and P = 1 atm).

Figure 14. CFD snapshots mapping the normalized total organic carbon concentration resulting from Euler
Lagrange simulations at different inlet bubble velocities uG = 0.1, 0.3, and 0.4 cm/s ([O3] = 40 g/Nm3, T = 20 °C, P = 1 atm).

Figure 13. Effect of inlet bubble velocities on the mineralization efficiency predicted by the Euler
Lagrange model at the half-height of bubble column (T = 20 °C, and P = 1 atm).

ozone mass transfer so we aim to analyze different computational mappings of total organic carbon attained at different heights of the bubble column reactor. The temporal evolution of total organic carbon conversion at different inlet ozone concentrations is shown in Figure 12 at the half-height of the column when the ozonation of phenolic wastewaters was simulated at T = 20 °C and P = 1 atm. After 30 min of reaction time, the computed results revealed the following conversions 20.1, 25.2, 26.9, and 27.4% for [O3] = 20, 40, 88, and 95 g/Nm3, whereas at the head of the bubble column we obtained 20.1, 25.7, 27.2, and 28.1% at [O3] = 20, 40, 88, and 95 g/Nm3, respectively. The EulerianLagrangian predictions are in reasonable agreement with experimental results so they are after 1 h being possible to improve the mineralization efficiencies up to 22.4, 26.3, 28.7, and 31.2% for [O3] = 20, 40, 88, and 95 g/Nm3, as can be seen in Figure 11. Figure 13 shows the temporal evolution of the total organic carbon conversion at half-height of the bubble column under different inlet bubble velocities when the ozonation of phenolic

wastewaters was simulated at T = 20 °C and P = 1 atm. As can be seen, the total organic carbon conversion was 17.7, 19.1, 20.4, and 22.7% for uG = 0.15, 0.25, 0.35, and 0.45 cm/s, respectively, after 1 h of reaction time. Here, we found comparable decontamination efficiencies, so the multiphase reaction system strongly depends on the homogeneous regime that has been observed with lower frequency for ozonation of low- to moderate-polluted wastewaters. The representative flow mappings of the normalized total organic carbon concentration are depicted in Figure 14 at different inlet bubble velocities uG = 0.15, 0.35, and 0.45 cm/s, [O3] = 40 g/Nm3, T = 20 °C, P = 1 atm. These computed results demonstrated diverse organic matter depuration levels along the bubble column reactor. Moreover, having analyzed the ozone concentration fields shown in Figures 6
8, the Euler
Lagrange model indicated a rapid intensification in ozone mass transfer at the bottom of the bubble column reactor mainly due to the higher degree of turbulence and higher shear rates. In fact, the characteristic dimensions of ozone bubbles were promoting the mass transport mechanism which translated to enhanced interfacial area available for the ozone mass transfer from the gas into the liquid phase. Consequently, both gas and liquid velocity plots were intimately related with the mineralization efficiencies that enable to advance the decontamination of organic liquid pollutants as demonstrated in Figure 14 by the interstitial flow patterns of the normalized total organic carbon concentration.

7. CONCLUSIONS The contemporary status of computational flow modeling is encouraging the application of modern CFD codes to the design and investigation of multiphase reactors. Aiming to accomplish 8900

dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research novel and stringent environmental regulations on the decontamination of high-strength wastewaters, we presented a comprehensive discrete bubble model for the ozonation of liquid pollutants. As long as a bubble column reactor can operate on a multitude of hydrodynamic regimes, we developed an Euler
Lagrange CFD model coupling the transport processes such as momentum and mass transfer, and ozonation kinetics. The model has been reasonably accounted for the dependency of flow regimes and extent of interaction among gas and liquid phases. The quantitative assessment between computed and experimental gas and liquid superficial velocities indicated that ozone velocity and concentration play a prominent role on the mineralization rate. We found that the process efficiency has been increased until it reveals a quasi-steady behavior which indicated that the reaction system happens to be noticeably reliant on oxidant concentrations under bubbly flow conditions. Finally, interstitial flow maps were effectively corroborated with total organic carbon concentration profiles as function of inlet bubble velocities and ozone concentrations.

’ AUTHOR INFORMATION Corresponding Author

*Phone: +351-239798723. Fax: +351-239798703. E-mail: rodrigo@ eq.uc.pt.

’ ACKNOWLEDGMENT The authors gratefully acknowledged the financial support of Fundac-~ao para a Ci^encia e Tecnologia, Portugal. ’ NOMENCLATURE A = interfacial area, m2 C = specie mass concentration, ppm d = diameter, m D = diffusivity, m2s
1 E = enhancement factor, dimensionless E€o = E€otv€os number, E€o=(FL-Fb)gdb2/σ, dimensionless F = force vector, N g = gravity acceleration, ms
2 H = Henry constant, dimensionless Ha = Hatta number, dimensionless k0 = reaction rate constant, min
1 kL,i = mass transfer coefficient, m s
1 I = unit tensor, dimensionless m = mass, kg m_ = mass transfer, kg s
1 P = pressure, Nm2
R = radius, m r = reaction rate Re = Reynolds number, Re=FL|uB-uL|db/μL, dimensionless Si = source term in the species balance equation Sc = Schmidt number, Sc= μL/(FLD), dimensionless Sm = volume averaged momentum transfer due to interphase forces Sh = Sherwood number, Sh = kLdb/D, dimensionless t = time, s uL = liquid velocity vector, ms
1 ub = bubble velocity vector, ms
1 V = volume, m3 Greek Letters

α = volume fraction, dimensionless μ = viscosity, kgm
1s
1

ARTICLE

F = density, kgm
3 σ = interfacial tension, Nm1
τ = stress tensor, Nm2
Subscripts

aq = aqueous b = bubble D = drag eff = effective G = gravity i = ith species L = liquid, lift P = pressure s = superficial T = turbulent VM = virtual mass * = interfacial equilibrium value

’ REFERENCES (1) Sanyal, J.; Vasquez, S.; Roy, S.; Dudukovic, M. P. Numerical simulation of gas-liquid dynamics in cylindrical bubble column reactors. Chem. Eng. Sci. 1999, 54, 5071–5083. (2) Joshi, J. B. Computational flow modeling and design of bubble column reactors. Chem. Eng. Sci. 2001, 56, 5893–5933. (3) Joshi, J. B.; Vitankar, V. S.; Kulkarni, A. A.; Dhotre, M. T.; Ekambara, K. Coherent flow structures in bubble column reactors. Chem. Eng. Sci. 2002, 57, 3157–3183. (4) Monahan, S. M.; Vitankar, V. S.; Fox, R. O. CFD predictions for flow-regime transitions in bubble columns. AIChE J. 2005, 51, 1897– 1923. (5) Pan, Y.; Dudukovic, M. P. Numerical investigation of gas-driven flow in 2-D bubble columns. AIChE J. 2000, 46, 434–449. (6) Rampure, R. M.; Buwa, V. V.; Ranade, V. V. Modeling of gas
liquid/gas
liquid
solid flows in bubble columns: Experiments and CFD simulations. Can. J. Chem. Eng. 2003, 81, 692–706. (7) Sokolichin, A.; Eigenberger, G. Gas
liquid flow in bubble columns and loop reactors: I. Detailed modeling and numerical simulation. Chem. Eng. Sci. 1994, 49, 5735–5746. (8) Delnoij, E.; Lammers, F. A.; Kuipers, J. A. M.; van Swaaij, W. P. M. Dynamic simulation of dispersed gas
liquid two-phase flow using a discrete bubble model. Chem. Eng. Sci. 1997, 52, 1429–1458. (9) Delnoij, E.; Kuipers, J. A. M.; van Swaaij, W. P. M. Dynamic simulation of gas
liquid two-phase flow: effect of column aspect ratio on the flow structure. Chem. Eng. Sci. 1997, 52, 3759–3772. (10) Lin, T.-J.; Reese, J.; Hong, T.; Fan, L.-S. Quantitative analysis and computation of two-dimensional bubble columns. AIChE J. 1996, 42, 301–318. (11) Delnoij, E.; Kuipers, J. A. M.; van Swaaij, W. P. M. Computational fluid dynamics applied to gas
liquid contactors. Chem. Eng. Sci. 1997, 52, 3623–3638. (12) Jakobsen, H. A.; Sannæs, B. H.; Grevskott, S.; Svendsen, H. F. Modeling of vertical bubble-driven flows. Ind. Eng. Chem. Res. 1997, 36, 4052–4074. (13) Kemiha, M.; Olmos, E.; Fei, W.; Poncin, S.; Li, H. Z. Passage of a gas bubble through a liquid
liquid interface. Ind. Eng. Chem. Res. 2007, 46, 6099–6104. (14) Roy, S.; Dudukovic, M. P. Flow mapping and modeling of liquid
solid risers. Ind. Eng. Chem. Res. 2001, 40, 5440–5454. (15) Diakov, V.; Varma, A. Optimal feed distribution in a packedbed membrane reactor: The case of methanol oxidative dehydrogenation. Ind. Eng. Chem. Res. 2004, 43, 309–314. (16) Lapin, A.; M€uller, D.; Reuss, M. Dynamic behavior of microbial populations in stirred bioreactors simulated with Euler
Lagrange methods: Traveling along the lifelines of single cells. Ind. Eng. Chem. Res. 2004, 43, 4647–4656. 8901

dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902

Industrial & Engineering Chemistry Research

ARTICLE

(17) Fraguío, M. S.; Cassanello, M. C.; Degaleesan, S.; Dudukovic, M. Flow regime diagnosis in bubble columns via pressure fluctuations and computer-assisted radioactive particle tracking measurements. Ind. Eng. Chem. Res. 2009, 48, 1072–1080. (18) Wang, X. L.; Li, T.; Lang, J.; Zhou, S. S.; Zhang, L. L.; Chen, M. X. Numerical analysis of solid
liquid two-phase flow on sandstone wastewater of hydropower stations in a rectangular sedimentation tank. Ind. Eng. Chem. Res. 2010, 49, 11714–11723. (19) Ali, B. A.; Pushpavanam, S. Experimental and computational investigation of two phase gas
liquid flows: Point source injection at the center. Ind. Eng. Chem. Res. 2011, DOI: 10.1021/ie2002268. (20) Chiang, P.-C.; Ko, Y.-W.; Liang, C.-H.; Chang, E.-E. Modeling an ozone bubble column for predicting its disinfection efficiency and control of DBP formation. Chemosphere 1999, 39, 55–70. (21) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena, 2nd ed.; Wiley: New York, 2002. (22) Martins, R. C.; Lopes, R. J. G.; Quinta-Ferreira, R. M. Lumped kinetic models for single ozonation of phenolic effluents. Chem. Eng. J. 2010, 165, 678–685. (23) Hoomans, B. P. B.; Kuipers, J. A. M.; Briels, W. J.; Swaaij, W. P. M. V. Discrete particle simulation of bubble and slug formation in a twodimensional gas-fluidised bed: A hard-sphere approach. Chem. Eng. Sci. 1996, 51, 99–118. (24) Lopes, R. J. G.; Quinta-Ferreira, R. M. Turbulence modelling of multiphase flow in high-pressure trickle-bed reactors. Chem. Eng. Sci. 2009, 64, 1806–1819. (25) Tomiyama, A.; Zun, I.; Higaki, H.; Makino, Y.; Sakaguchi, T. A three-dimensional particle tracking method for bubble flow simulation. Nucl. Eng. Des. 1997, 175, 77–86. (26) Pfleger, D.; Gomes, S.; Gilbert, N.; Wagner, H.-G. Hydrodynamic simulations of laboratory scale bubble columns fundamental studies of the Eulerian
Eulerian modeling approach. Chem. Eng. Sci. 1999, 54, 5091–5099.

8902

dx.doi.org/10.1021/ie2021366 |Ind. Eng. Chem. Res. 2012, 51, 8891–8902