Flow field and residence time distribution simulation of a cross-flow gas–liquid wastewater treatment reactor using CFD

Chemical Engineering Science 63 (2008) 2436 – 2449 www.elsevier.com/locate/ces

Flow field and residence time distribution simulation of a cross-flow gas–liquid wastewater treatment reactor using CFD Yann Le Moullec, Olivier Potier, Caroline Gentric ∗ , Jean Pierre Leclerc Laboratoire des Sciences du Génie Chimique, UPR 6811, CNRS-ENSIC-INPL, 1, rue Grandville BP 20451, 54001 Nancy, France Received 16 January 2007; received in revised form 22 January 2008; accepted 24 January 2008 Available online 6 February 2008

Abstract A three-dimensional Eulerian–Eulerian two-phase approach has been used for the simulation of a cross-flow gas–liquid wastewater treatment reactor. Two different turbulence models have been tested: the k. and Reynolds Stress Model (RSM) models. Bubble induced turbulence source terms have been added to these models. Numerical results have been validated using Laser Doppler Velocimetry (LDV) measurements. Simulations with both turbulence models successfully predicted the hydrodynamics of the reactor. Then particle tracking with a stochastic approach has been used to calculate residence time distributions (RTD) with the flow previously simulated. It has been shown that dispersion in the reactor is primarily due to turbulence. Results have been compared with experimental RTD for various liquid and gas flowrates both on a bench scale and full scale plant. The RSM model accurately predicted the dispersion whereas the standard k. model slightly underestimated the dispersion. 䉷 2008 Elsevier Ltd. All rights reserved. Keywords: Dispersion; Multiphase flow; RTD; Simulation; Turbulence; Wastewater treatment

1. Introduction The pollution removal by microorganisms is an essential step in the biodegradable wastewater treatment plants. These biological reactions often take place in gas/liquid reactors with a very long length compared to their height and width (channel reactors). In this type of reactors, the global water flow is along the length of the reactor and gas is sparged at the bottom. In wastewater treatment reactors, apparent reaction orders are often greater than zero; therefore the efficiency of the pollution removal reaction depends on the hydrodynamics (Levin and Gealt, 1993). Therefore, in order to accurately predict conversions, the kinetics modelling must be coupled with the hydrodynamics description. But the prediction of bubbly flow hydrodynamics is quite difficult due to the complex phenomena involved. Nevertheless, with the improvement of computational power and the availability of computational fluid dynamics (CFD) codes adapted to such two-phase flows, gas–liquid ∗ Corresponding author. Tel.: +33 3 83 17 53 38; fax: +33 3 83 32 29 75.

E-mail address: [email protected] (C. Gentric). 0009-2509/$ - see front matter 䉷 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2008.01.029

reactors have been studied more extensively. In particular bubble column reactors have been extensively investigated either with the Euler–Euler approach (Jakobsen et al., 2005; Pfleger and Becker, 2001; Sokolichin et al., 2004) or with the Euler–Lagrange approach (Lain et al., 2002; Lapin and Lübbert, 1994). Two- or three-dimensional (2D and 3D) and stationary or transient simulations have been performed. Although transient 3D simulations are necessary to represent the complex flow structure of bubble columns, 2D and stationary simulations allow correct representations of the average flow field and can be sufficient in certain cases (Sokolichin et al., 2004). But simulation of bubbly flows using CFD is not straightforward. Closure terms, in particular interaction forces between phases and their impact on the simulation results have been the subject of numerous studies (Jakobsen et al., 1997). Another critical point widely discussed, is the turbulence closure terms. Mainly three modifications of single phase turbulent viscosity models have been proposed and compared in the literature: the addition of a bubble induced turbulence viscosity term (Sato and Sekoguchi, 1975), a pseudo-bubble induced turbulence tensor (Arnold, 1988) and the addition of bubble induced turbulence source terms in the k. transport

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

equations (Gu and Guo, 2005; Mudde and Simonin, 1999; Olmos et al., 2003; Pfleger and Becker, 2001). It seems that the addition of bubble induced turbulence source terms is the best compromise between precision and complexity (Sokolichin et al., 2004). Recently, efforts have been made to simulate local and global mass transfer in bubble columns, air-lift and trickle bed reactors (Ekambara et al., 2005; Gunjal et al., 2003; Vial et al., 2005). CFD simulations of liquid phase residence time distribution (RTD) in these gas–liquid reactors have also been performed with good results (Ekambara et al., 2005; Gunjal et al., 2003). With these improvements in gas–liquid simulations, some other industrial reactors can be investigated and specifically wastewater treatment reactors. The previously presented studies are limited to vertical liquid flow with counter or co-current gas flow. Wastewater treatment reactors by activated sludge are mostly cross-flow reactors in which liquid flow is horizontal and gas flow is vertical. Two types of reactors are used: closed-loop oxidation ditch reactors, also called carrousel reactors, and channel aerators. In carrousel reactors, the liquid horizontal velocity is such that gas has no real impact on the hydrodynamics and mixing is sometimes performed with impellers (Roustan and Line, 1996). CFD modelling of the velocity field and mass transfer has been carried out for these reactors (Cockx et al., 2001). In channel aerators, the liquid horizontal velocity is very low compared to carousel reactors and mixing is performed by gas injection in the reactor. Therefore an accurate simulation of the interaction between bubbles and liquid is of primary importance. The purpose of this work is to obtain an accurate and predictive hydrodynamics model of a channel aerator with respect to velocity, turbulence and RTD which is, as in most gas–liquid flows, governed by local scale phenomena (Vial et al., 2005). Validation is performed with Laser Doppler Velocimetry (LDV) measurements and RTD experiments. 2. Experimental setup and preliminary results 2.1. Bench scale reactor Experiments have been carried out in a bench scale channel reactor filled with tap water. It is built in transparent material

2437

(Plexiglas) to allow LDV measurements. Its total length is 3.6 m with a rectangular section of width and height, respectively, equal to 0.18 and 0.2 m (Fig. 1). One side of the walls of the reactor is fitted with stainless-steel tubes where 1 mm holes have been drilled every centimetre for air sparging. Gas and liquid flowrates are the variable parameters. 2.2. LDV measurements 2D LDV measurements have been performed to obtain axial (Ux ), lateral (Uy ) and vertical (Uz ) time-averaged velocity fields (see Fig. 1 for directions) as well as an estimated value of the turbulent kinetic energy (k). An acquisition time of 180 s is required in order to obtain statistically meaningful results for all measurement positions. Photographic technique has shown that the gas hold-up along the reactor length is inhomogeneous though the variations are small. In order to verify that these small variations of gas holdup are of no importance for the LDV measurements, these measurements have been carried out at two different planes with different apparent gas hold-ups (planes 1 and 2: one is located between two sparger tubes, the other is located at the middle of a sparging device, see Fig. 1). 2.3. RTD experiments Tracer experiments used for simulation validation have been performed in a previous study in the same bench scale reactor (Potier et al., 2005). Approximately 80 RTD experiments have been carried out by injecting pulses of a solution of sodium chloride (NaCl) at the inlet and monitoring its concentration with a conductimetric probe at the outlet. A plug flow with axial dispersion model (Villermaux, 1993) allows a very precise representation of the data. Mean residence times and axial dispersion coefficients have been determined by model fitting with the software DTSPRO 4.2 (PROGEPI, France). This previous study has shown the existence of a correlation between the axial dispersion coefficient on the one hand and the gas flowrate and the geometry of the reactor on the second hand:  Qg −1.99 D = (0.2032.w − 0.008569) (1) (100h)0.00476.w L

Fig. 1. Schematic representation of the bench scale reactor.

2438

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

with h, w and L, respectively the height, width and length of the reactor and Qg the gas flowrate. This correlation has been established for Qg (15.65 L min−1 ), Ql (1.5.5 L min−1 ), h (0.08.0.2 m) and w (0.05.0.2 m). Besides, this correlation has been validated on an industrial reactor successfully. In the present work, results from simulation and experiments have been compared for different gas and liquid flowrates for a given geometry: length of 3.6 m, height of 0.2 m and width of 0.18 m. 3. Hydrodynamics modelling

Table 1 Constants used in the two-phase k . model C

C1

C2

Ck

C

k



0.09

1.44

1.92

1.44

1.0

1.0

1.3

The drag coefficient is calculated from the following equation (Jamialahmadi et al., 1994): 4 (l − g ) gd g with 3 l u2∞  gd g 2l + for dg > 3 mm. u∞ = dg (l + g ) 2

Cd = Because of the high number of bubbles (approximately 7000 and more in an industrial size reactor) in the reactor and despite the low gas fraction (< 10%), an Euler–Euler approach has been used. CFD simulations have been carried out with the CFD software FLUENT.

(9)

3.2. Two-phase k. turbulence model 3.1. Euler–Euler equations The mass transfer between phases has been neglected therefore the continuity equation is expressed by j (q q ) + ∇ · (q q Uq ) = 0 jt

with q = g or l

(2)

besides: g + l = 1.

(3)

The momentum conservation equation for multiphase flows is given by

= l (t,l ∇Ul [∇Ul + (∇Ul )T ] − l l ) + k,l , (10)     t,l j ∇l (l l l ) + ∇ · (l l Ul l ) − ∇ · l l + jt 

j (q q Uq ) + ∇ · (q q Uq Uq ) jt = −q ∇p + q q g + ∇ · q (q + t,q ) + Iq with q = g or l,

(4)

where Iq is the interphase momentum exchange term and q and t,q are, respectively, the viscous stress tensor and the turbulent stress tensor, defined as follows: q = q (∇Uq + ∇Uqt ) + (q − 23 q )∇.Uq I

(5)

and t,q = t,q (∇Uq + ∇Uqt ) − 23 (kq + t,q ∇.Uq )I kq2 q

.

(11) The turbulent viscosity of the continuous phase is calculated by Eq. (7). The bubble induced turbulence source terms are calculated from (Pfleger and Becker, 2001)

,l = l (7)

and Ig = −Il .

l (C1 t,l ∇Ul [∇Ul + (∇Ul )T ] − C2 l l ) + ,l . kl

k,l = l Ck |Il | · |Ug − Ul |,

In most gas–liquid flows virtual mass force is necessary to reproduce transient flow characteristics (Leon-Becerril et al., 2002), nevertheless this force has no influence on stationary calculation results (Sokolichin et al., 2004; Zhao et al., 2006) therefore it has been neglected and only the drag force has been considered significant. Iq is thus defined as follows: 3 g l |Ug − Ul |(Ug − Ul ) Il = Cd 4 dg

= l

(6)

with t,q = q C

In order to model turbulence, a two-phase k. turbulence model can be used (Elghobashi and Abou-Arab, 1982; Mudde and Simonin, 1999). Only the continuous phase turbulence is taken into account and is considered isotropic. The dispersed phase influences the continuous phase through bubble induced turbulence terms. For the liquid phase, the transport equations for the turbulent kinetic energy k and the turbulent dissipation rate  are     t,l j ∇kl (l l kl ) + ∇ · (l l Ul kl ) − ∇ · l l + jt k

(8)

l C1 C |Il | · |Ug − Ul |. kl

(12) (13)

All the constants of the two-phase k. turbulence model are summarized in Table 1. 3.3. Two-phase RSM turbulence model The Reynolds Stress Model (RSM) (Daly and Harlow, 1970) does not assume the isotropy of turbulence and is more adapted to swirling flows than the k. turbulence model and its variants (Ranade, 2002). In this model, the required computational time is only 60% higher than with a k. model. As for the k. model, turbulence is only considered for the continuous phase.

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449 Table 2 Constants used in the two-phase RSM model

2439

Table 3 Constants used in the stochastic model

C

C 1

C2

Ck

C



CL k . model

CL RSM model

0.09

1.44

1.92

1.44

1.92

1.0

0.15

0.3

Bubble induced turbulence source terms are added to the transport equations. The transport equation for the Reynolds stress tensor component Rij is j j (l l Rij ) + (l l Uk Rij ) jt jxk     jUj jRij j jUi + l l = −l l Rik + Rj k jxk jxk jxk jxk   juj jui j    − (l l ui uj uk ) + l p + jxk jxj jxi − l l ij + R,ij

(14)

dup 3  = Cd l |up −ul |(up −ul ) dt 4 dp p with

with R,ij =

turbulent dispersion of particles is also taken into account: particles in a turbulent flow are submitted to a randomly varying flow field. This phenomenon is modelled by a stochastic approach, assuming that the particles interact with a sequence of turbulent eddies, the fluctuating velocity within each eddy being isotropic and obeying a Gaussian probability density function. The interaction time I is assumed to be sufficiently short so that the fluid velocity in a given eddy is constant during this process. The particle trajectory is therefore calculated as follows:

2 3 ij k,l

with k,l from Eq. (12)

(15)

(19)



ul =Ul +ul

and

ul

composed of

ux,l =uy,l =uz,l =

and ij = 23 ij l

with l from Eq. (17).

(16)

The third and fourth term on the right-hand side of Eq. (14) are modeled (see Fluent Inc., 2005). The turbulence dissipation rate l is calculated using the following transport equation:     t,l j j j j l  l + (l l l ) + (l l Ui l ) − l jt jxi jxj  jxj     l jui jui =−l C1 l Rik +C2 l l + ,l (17) +Rik kl jxk jxk with kl = 21 (Rii + Rjj + Rkk )

2 kl . 3 (20)

 is a normally distributed random number and Cd is calculated by the Morsi and Alexander correlation (1972). A new value of  is applied each time the interaction time is reached. This interaction time is taken as the minimum of the eddy lifetime e and the eddy crossing time ct . e = −CL

kl log( ), l

where is a uniform random number between 0 and 1.    Le ct = − ln 1 − |up − ul |

(21)

(22)

with the particle relaxation time: (18)

and ,l is given by Eq. (13) and t,l by Eq. (7). All the constants of the two-phase RSM turbulence model are summarized in Table 2. 3.4. Lagrangian particle motion model In order to perform the simulation of RTD, particle tracking has been used (Stropky et al., 2007; Thyn et al., 1998). RTD is deduced from the simulation of a pulse injection of a sufficient number of particles. The trajectory of a particle is predicted by integrating the force balance on the particle in a Lagrangian reference frame (Fluent Inc., 2005). The characteristics of these particles have been chosen such as they can be considered as perfect tracers of the fluid flow. Their diameter is 10−6 m and they have the same density as water. This ensures that they can perfectly follow all the simulated timescales. The

=

2p dg 4 3 2l Cd |up − ul |

and the eddy lifetime:  3/2 3 kl . C Le = 2 l Different values for CL are recommended depending on the turbulence model (Fluent Inc., 2005) (Table 3). 3.5. Transport equation of a tracer Another approach used to simulate RTD is the solution of a transport equation for a passive tracer with the same physical properties as the continuous phase (Talvy et al., 2007; Zhang et al., 2007). The transport equation for the concentration of a

2440

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

tracer in a turbulent flow is

   j  l Dm + t ∇Ctr (l Ctr ) + ∇.(l Ul Ctr ) = ∇. jt Sct

Therefore a refined mesh is necessary near that wall. Different grid sizes have been tested: cubic cells of 1 cm3 (130 000 cells) or 0.125 cm3 (1 000 000 cells) and cells of 1 cm3 with a refined mesh near the walls (350 000 cells) (Fig. 2). Grid independency has been verified (Fig. 3) and finally this last grid offers the best compromise between precision and computational effort. Inlets of gas and liquid have been modelled with a velocity inlet boundary condition. The turbulence boundary conditions at the inlets are given through the turbulence intensity (10%) and inlet hydraulic diameter. The outlet of liquid has also been modelled with a velocity inlet boundary condition; it allowed us to impose the outflow velocity in order to avoid backflow which causes mass balance problem when RTD simulations are performed. In order to simulate the degassing boundary condition at the top of the reactor, an outlet boundary condition has been used for the gas phase and a symmetry boundary condition has been used for the liquid phase, this degassing condition has been implemented through a user defined function. A second order discretization scheme (QUICK: Quadrative Upwind Interpolation for Convective Kinematics) has been chosen for the turbulent dissipation rate and void fraction equations in order to limit their numerical diffusion (calculation errors induced by grids and discretization schemes). A first order discretization scheme has been tested for the solution of the momentum equations but it leads to unstable results; consequently a second order discretization scheme (QUICK) has been used to solve all equations. The SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) pressure–velocity coupling scheme has been used. The under relaxation parameters used are given in Table 4.

(23)

with Sct , the turbulent Schmidt number considered constant and equal to 0.7. The simulation of the RTD is obtained following the outlet concentration of the tracer injected as a pulse at the inlet. 3.6. Grid, boundary conditions and discretization scheme It is commonly admitted in bubble induced flows that there is no need for boundary layer or advanced wall treatment (Troshko, 2006). In our case, bubbles are only located near one wall which induces a high velocity gradient at this location.

Table 4 Under relaxation parameter used for each variable

Time-averaged liquid vertical velocity Uz (m/s)

Fig. 2. Projection of the grid on a vertical plane (y, z).

Pressure

Momentum

Volume fraction

k



Rij

0.3

0.7

0.2

0.8

0.8

0.5

0.4 0.3 0.2 0.1 0 0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

-0.1 -0.2

129600 cells grid 345600 cells grid

-0.3

1036800 cells grid

-0.4 y (m)

Fig. 3. Simulated average liquid vertical velocity on a transversal line (x = 1.8 m, z = 0.1 m) for three different grid sizes.

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

2441

4. Results This CFD study will be divided into three main parts. First, a simulation of the complete turbulent hydrodynamics for different liquid and gas flowrates is carried out. The second step is the simulation of the RTD. And finally, the numerical parameters used for the bench scale reactor are validated on a plant scale reactor. A test case experiment is chosen as a reference for the validation of the simulations. It is carried out for a liquid flowrate QL of 3.6 L min−1 and a gas flowrate QG of 15 L min−1 . The bubble diameter for this air–water system has been estimated using a photographic technique: an average value of 4 mm has been used and since the gas hold-up is very low (< 10%), bubble break-up and coalescence are neglected. 4.1. Velocity profiles and turbulence characteristics Simulations have been performed with the k. and RSM turbulence models. Fig. 4 presents the projection of the liquid velocity field on a vertical plane for the experimental data and the results of simulations with the k. and RSM turbulence models. Simulated velocity fields with the k. and RSM models are superimposed, both models giving similar results. The overall agreement is good, however, some discrepancies can be observed near the bubble injection position and near the free surface. Discrepancies observed near the sparger are probably due to the simplification made for the gas inlet boundary conditions (velocity magnitude of 0.162 m s−1 and void fraction of 1) whereas those observed near the free surface are probably due to the simplified representation of the surface which is not planar in the experimental setup. Fig. 5 compares the measured and simulated profiles of the vertical liquid velocity at three different heights. These quantitative comparisons show some reasonable agreement between experiments and simulations (an average difference of 0.035 m s−1 is observed which represent 18% of the maximum liquid velocity) whereas the two simulations results are very close. The LDV measurements on the two measurement planes (defined in Section 2.1) show that there is no significant difference, with respect to the vertical liquid velocity, between these two positions despite their apparent different gas hold-ups (experimental data planes 1 and 2). As it will be discussed later, despite the discrepancies between experimental and simulated results, the liquid velocity fields calculated are satisfying for RTD calculations. LDV measurements revealed a non-uniform time-averaged axial velocity profile along the length of the reactor (Fig. 6). Simulations predict these spatial oscillations (Fig. 6) but fail to predict the period and amplitude of the phenomenon. The maximum experimental value of these oscillations (0.02 m s−1 ) is ten times the mean axial velocity (0.0017 m s−1 ). Because of these very small velocity values, measurement uncertainty is high. But, at the same time, the maximum axial velocity is still ten times smaller than the maximum vertical velocity (0.2 m s−1 ). The amplitude of these fluctuations is thus negligible compared to the highest velocities encountered but

Fig. 4. Overall representation of the experimental and simulated average velocity fields on a vertical plane (y, z) for both turbulence models.

may have an impact on dispersion; this will be investigated in Section 4.2. Table 5 presents results relative to turbulence. The experimental results for turbulent kinetic energy k are surface averages on the measurement plane (plane 1). It can be highlighted that experiments have shown an isotropic turbulence. The turbulent dissipation rate  of the liquid phase can be estimated with an energy balance (Eq. (24)), assuming that all the energy brought by the isothermal expansion of the gas is dissipated by turbulence:   Patm QG  gh l theoretical = . (24) ln 1 + l l V Patm Simulation results for the turbulent kinetic energy k are also surface averages on plane 1 and turbulent dissipation rate  is a volume average over the whole reactor. The turbulent kinetic energy k is estimated by simulation with an underestimation of around 26% by the k. model and a 33% overestimation by the RSM model. The turbulent dissipation rate  obtained with the k. model is also underestimated by 20% whereas the turbulent dissipation rate  obtained with the RSM model is more importantly underestimated by 25%. The underestimation of  could be explained by the energy dissipation at the free surface due to bubble break-up. This dissipation can reach 30% of the total energy for some configurations (Desai et al., 1995). 4.2. Residence time distribution All our experimental RTD can be described by the plug flow with axial dispersion model. Therefore the dispersion coefficient is a good estimation of the overall mixing within the reactor. The simulated dispersion coefficient can be described as the sum of four contributions, three of them are

2442

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

Fig. 5. Comparison between experimental and simulated average vertical velocities at three different heights: 5, 10 and 15 cm. 0.025 Experimental data

0.02

RSM simulation

0.015

k-eps simulation

Ux (m/s)

0.01 0.005 0 -0.005 -0.01 -0.015 -0.02 -0.025 0

5

10

15

20

25

30

35

40

45

50

x (cm)

Fig. 6. Comparison between experimental and simulated average axial velocities along a fraction of the length of the reactor (at y = 4 cm, z = 10 cm).

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

physical: molecular diffusion, dispersion due to velocity field (convection) and turbulent dispersion, and the fourth term is the numerical diffusion (Eq. (25)). Dtotal = Dmolecular + Dconvection + Dturbulent + Dnumerical . (25) There are two main methods to simulate a RTD with CFD: • Particle tracking: the motion equations for numerous particles are solved from the entrance of each particle in the reactor to its escape from the reactor (Section 3.4). A significant statistical RTD can be obtained as long as enough particles trajectories are calculated (approximately 4000 in our case). Using this method avoids solving a transport equation and, therefore, eliminates numerical diffusion because numerical diffusion is due to the discretization of a transport equation (numerical schemes and mesh size). However, it does not allow to take the molecular diffusion of a tracer into account. This has no consequence here since the molecular diffusion of our experimental tracer is negligible compared to the dispersion coefficient. This will be verified later. • Tracer transport: the transport equation for a passive tracer is solved (Section 3.4). A pulse of tracer is injected at the inlet and its concentration is followed at the water outlet. This method includes theoretically all the dispersion sources. Its implementation can be performed in two steps in order to reduce the calculation time: first the solution of the stationary flow field and second the solution of the transient transport equation for the tracer. Both methods have been tested in this work. Table 5 Comparison between experimental and simulated turbulence variables Experimental/

k .

theoretical

simulation

simulation

k (10−3 m2 s−2 )  (10−3 m2 s−3 )

2.42 (exp) 8.75 (Eq. (25))

1.78 6.84

2443

4.2.1. RTD simulation with scalar transport An integration time step of 4 s (0.2% of mean residence time) has been selected. Fig. 7 presents the RTD obtained with the RSM turbulence model and the experimental one. Simulation results obtained with the k. turbulence model are not presented because they are similar to the ones obtained from RSM simulations. The mean residence time is correctly simulated whereas the dispersion coefficient is 30% underestimated. The simulated dispersion coefficient is 4.26 × 10−4 m2 s−1 whereas the experimental one is 6.01 × 10−4 m2 s−1 . 4.2.2. Influence of the different dispersion sources on the RTD simulated via tracer transport Talvy et al. (2007) proposed a method to identify the part of the two different aspects of the dispersion: Dconvection and Dturbulent (the two other aspects are negligible in our case: Dmolecular = 1.48 × 10−10 m2 s−1 and an estimation of the numerical dispersion gives Dnumerical = 5.0 × 10−6 m2 s−1 ). Eq. (26) proposes an estimation of the turbulent dispersion (Talvy et al., 2007): Dturbulent =

t . l Sct

(26)

In our case this turbulent dispersion is equal to 3.81 × 10−4 m2 s−1 , approximately 90% of the total dispersion. The convective dispersion can be estimated from CFD simulation with Eq. (27) (Talvy et al., 2007).  ul,x Ctr y,z − ul,x y,z Ctr y,z Dconvection = . (27) jCtr y,z jx

RSM

3.31 6.10

z

This equation estimates the convective dispersion to 1.14 × 10−5 m2 s−1 , approximately 2% of the total dispersion. These results indicate that the convective dispersion is negligible compared to the turbulent dispersion. This conclusion will be discussed later. Thus the main dispersion source is turbulence.

0.0006 Experimental RTD

0.0005

Simulated RTD realised with 100 s sample count (particles tracking) Simulated RTD realised with 200 s sample count (particles tracking)

E (t)

0.0004

Simulated RTD fitted with DTS PRO (particles tracking)

0.0003

Simulated RTD with RSM model (transport equation)

0.0002

0.0001

0.0000 0

1000

2000

3000

4000

5000

6000

Time (s)

Fig. 7. Comparison between experimental and simulated RTD obtained with the RSM turbulence model and the transport equation method or the RSM turbulence model and the particle tracking method for different sample times for a liquid flowrate of 3.6 L min−1 and a gas flowrate of 15 L min−1 .

2444

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

Consequently the turbulent Schmidt number allows to adjust the dispersion coefficient. It has been verified that a turbulent Schmidt number of 0.44 would give a perfect match between experimental and simulated RTD. Table 6 Numerical experiments: Comparison of the influence of the main parameters on the simulated dispersion coefficient with the k . turbulence model Simulation method

Ref

Characteristics

D (10−4 m2 s−1 )

Particle tracking

1 2 3 4 5

Reference 2.k 2. 2.uy , 2.uz jux /jxi = 0

2.51 5.02 3.18 3.16 2.37

Transport of passive tracer

6

Dm = 10−20 m2 s−1

3.34

7

Dm = 14.8 10−10 m2 s−1

3.35

In order to verify the negligibility of the molecular diffusion, two simulations have been performed: two RTD simulations with tracer transport method, the first one with a molecular diffusion coefficient of 10−20 m2 s−1 (in order to simulate a non-diffusive tracer) and the second one with a molecular diffusion coefficient of 1.48 × 10−10 m2 s−1 (NaCl in pure water at 20 ◦ C). The results of these two experiments are presented in Table 6 (Ref 6 and 7) and they show that the molecular diffusion has no significant impact on the overall dispersion. 4.2.3. RTD simulation with particle tracking Experimental molecular diffusion is negligible compared to the turbulent and convective dispersions therefore particle tracking can also be used to estimate RTD. Because of the discrete nature of particle tracking, particles must be sampled in order to exploit data: particles escaping the reactor during the same time period are counted; RTD is obtained from the curve giving

0.0008 Experimental RTD

0.0007

Simulated RTD with the RSM model fitted with DTS PRO

0.0006

Simulated RTD with the k-eps model fitted with DTS PRO

E (t)

0.0005 0.0004 0.0003 0.0002 0.0001 0 1000

0

2000

3000 Time (s)

4000

6000

5000

Fig. 8. Comparison between experimental and simulated RTD obtained with the RSM and the k . turbulence models and the particle tracking method for a liquid flowrate of 3.6 L min−1 and a gas flowrate of 15 L min−1 .

12

45 40

D (10-4 m2/s)

35 8

30 25

6 20 4

15 Simulated dispersion coefficient Experimental dispersion coefficient Correlated dispersion coefficient (Potier 2005) Simulated residence time Experimental residence time

2

10

Residence time (min)

10

5

0

0 0

10

20

30 40 Gas flowrate (l/min)

50

60

70

Fig. 9. Experimental and simulated dispersion coefficients and mean residence times for different gas flowrates with a constant liquid flowrate of 3.6 L min−1 with the RSM model and the particle tracking method.

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

2445

50

18 16

40

14 12

30

10 8

20

6

Simulated dispersion coefficient Experimental dispersion coefficient Simulated residence time Experimental residence time

4 2 0 2

2.5

3 3.5 Liquid flowrate (l/min)

4

10

Mean residence time (min)

Dispersion coeffient (10-4 m2/s2)

20

0 4.5

Fig. 10. Experimental and simulated dispersion coefficients and mean residence times for different liquid flowrates with a constant gas flowrate of 65 L min−1 with the RSM model and the particle tracking method.

the number of particles counted during one period as a function of time. Fig. 7 shows the curves obtained with the RSM turbulence model for two different sampling times: 100 and 200 s and an experimental curve fitted with the plug flow with closed–closed axial dispersion model. As the sampling time increases, fluctuations decrease because of an integration effect. The mean residence time and dispersion coefficient D are determined using the raw data: first, the mean residence time is determined and in a second step the dispersion coefficient is optimized in order to fit a plug flow model with axial dispersion. Comparison between experimental and simulated RTD is shown in Fig. 8 for the k. and RSM turbulence models. The RSM turbulence model simulation leads to very good results and the k. turbulence model leads to an underestimated value of the dispersion coefficient of around 50%. These results are due to the stochastic model and especially the constant CL (Table 3). In fact, in the simulation of turbulent dispersion, there is a direct linear relation between CL and the dispersion coefficient D. Consequently the value of CL can be optimized to obtain a good dispersion coefficient. Empirical value of CL = 3C = 0.27, which is close to the one used for the RSM model (Table 3), can be found in literature for the k. turbulence (CD-Adapco-Group, 2003); with this value of CL both the k. and RSM turbulence models give the same results (not shown here). Additional experiments and simulations have been performed for different gas and liquid flowrates and similar good agreement between experimental and simulated RTD with the particle tracking method have been obtained. All simulations were carried out with the RSM turbulence model. Comparison between experimental and simulated dispersion coefficients and mean residence times for different gas flowrates are presented in Fig. 9. The dispersion coefficient seems to depend on Qg according to a square root law, the same as in correlation (Eq. (1)). It can be highlighted that the calculations diverge when the gas flowrate increases beyond 70 L min−1 . The same comparison, with different liquid flowrates, is presented in Fig. 10. A good agreement can also be observed for

the dispersion coefficient even if it does not vary much with the liquid flowrate in this domain. The mean residence time is also correctly simulated. 4.2.4. Influence of the different dispersion sources on the RTD simulated via particle tracking Simulated RTD with transport equation solution shows that a very high proportion of the axial dispersion is due to turbulent dispersion (90%). But the solution of the transport equation fails to accurately simulate experimental RTD. On the other hand, particle tracking allows a better RTD estimation, but the determination of the different dispersion sources is not as straightforward. With CFD it is possible to investigate the influence of one of the dispersion source at once, with the other parameters kept constant. This is impossible with physical experiments. These “numerical experiments” do not have real physical meaning but they allow a better understanding of the physical phenomena responsible for the dispersion in the reactor. These numerical experiments have been performed based on our reference case (Sections 4.1 and 4.2). Our main objective, here, is to determine which of turbulence, convection or molecular diffusion is responsible for the observed dispersion. All the following “numerical experiments” are performed with the k. turbulence model in order to avoid too long calculations. To study the impact of the velocity field, two simulations have been performed: the first one with an axial velocity set to a constant value of 1.667 mm s−1 (i.e. the average value of the axial liquid velocity) in all grid cells (jux /jxi = 0) in order to evaluate the impact of axial velocity oscillations, the second one with the velocity of the gas induced loop doubled (see Fig. 3) in order to evaluate the contribution of these loops in the dispersion (i.e. uz = 2uz,ref , uy = 2uy,ref ). Results, presented in Table 6 (Ref 5 and 4), show that the axial velocity oscillations (see Fig. 6) have negligible influence on D and gas induced loops have a quite significant impact on D (increase of approximately 25% by doubling the loop velocity).

2446

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

Fig. 11. Projection of the simulated liquid velocity field on a vertical plane for the industrial reactor (x = 51 m).

Time-averaged vertical liquid velocity Uz (m/s)

1.2 1 91936 cells grid

0.8

151776 cells grid 194208 cells grid

0.6 0.4 0.2 0 -0.2

0

1

2

3

4

5

6

7

8

9

-0.4 -0.6 -0.8 -1 y (m)

Fig. 12. Simulated average vertical liquid velocity on a transversal line (x = 51 m, z = 1.8 m) for three different grid sizes.

Finally, in order to study the impact of turbulence, a simulation with a doubled value of  and another with a doubled value of k has been performed. Results are also presented in Table 6 (ref 2 and 3). The modification of  induces a small difference on D (increase of around 25%) whereas doubled the value of k doubles the value of D (impact of 100%). D seems to be strongly linked to turbulence and especially to the value of the turbulent kinetic energy k in accordance with Eq. (26) where t is given by Eq. (7). With all these “numerical experiments”, we validate the fact that turbulence is the main cause of the overall dispersion. Turbulent kinetic energy k is the variable which mainly impacts on the overall dispersion. But the particle tracking method shows that convective dispersion is not negligible, contrary to the conclusion of the tracer method. We can conclude that the RSM turbulence model coupled with the particle tracking method gives very good RTD results. As it has been said in the RTD simulation part of Section 4.2, CL is the preponderant constant in the modelling of turbulent dispersion therefore with an empirically deduced CL such as CL = 3C , good RTD results can also be achieved with the k. turbulence model.

4.3. Scale up on industrial size reactor In order to validate the chosen hydrodynamics model, RTD simulations of a plant scale reactor have been carried out. RTD experiments (Potier et al., 2005) have been conducted on a 3300 m3 (102 m long, 9 m wide, 3.6 high) channel reactor, at the Nancy Maxéville plant (France). Lithium chloride has been used as an inert tracer and injected as a pulse near the liquid inlet. Lithium chloride concentration was measured by atomic absorption. It was verified that the tracer was not absorbed on solid matter. The liquid flowrates were 1280, 1650 and 2200 m3 h−1 and the linear gas flowrate was approximately 0.0175 m3 m−1 s−1 ; aeration takes place on both sides of the reactor which causes a double gas induced loop, as seen in the simulated flow field presented in Fig. 11. The vertical plane (x, z) in the middle of the reactor is a symmetry plane therefore a symmetry boundary condition is used in order to shorten calculation time. The grids for this kind of wastewater plant scale reactor simulations are very coarse (average cell volume of 0.016 m3 ) but give satisfactory results (Glover et al., 2006; Stropky et al., 2007). Grid independency has been verified (Fig. 12) and a 151 776 cells grid has been used with

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

2447

0.00030 experimental low flowrate simulated low flowrate experimental medium flowrate simulated medium flowrate experimental high flowrate simulated high flowrate

0.00025

E (t)

0.00020

0.00015

0.00010

0.00005

0.00000 0

4000

8000

12000

16000

20000

Time (s)

Fig. 13. Comparison between experimental and simulated RTD obtained with the RSM model and the particle tracking method on a full scale reactor for three different liquid flowrates.

boundary conditions similar to those used for the bench scale reactor. The physical properties of the mixture of wastewater and sludge are approximated to those of water, which is acceptable because of the low concentration of sludge (less than 8 g L−1 ) in the reactor and the bubble diameter is approximated to 1 cm as observed in experiments. Experimental results have been treated by curve fitting with a plug flow with axial dispersion model. Simulated results have been obtained as in Section 4.2 with a particle tracking method. This method avoids possible numerical diffusion due to the coarse grid. Fig. 13 shows very good agreement between experimental and simulated data. There is less than 5% of error on the mean residence time and the dispersion coefficient. It can be highlighted that the greater discrepancy is observed with the lowest flowrate, which is due to the difficulty to keep constant the inlet flowrate of a real wastewater treatment plant during such a long period. In fact, some small flowrate variations have been noticed. 5. Conclusion 3D Euler–Euler numerical simulations of a cross-flow gas–liquid reactor have been presented. These simulations have been carried out using the CFD code FLUENT. Grid dependence, discretization schemes and different turbulence models have been tested. It has been found that the use of the RSM turbulence model coupled with a second order discretization scheme allows to successfully predict the dispersion coefficient of the reactor. RTD simulations with particle tracking method with standard parameters give more accurate results than those with transport equation solution. But it is also possible to adjust the simulation to the experiment with the modification of key model parameters (Sct or CL ). CFD allowed to investigate the different effects responsible for the dispersion. It has been shown that the dispersion is mainly due to turbulence but the dispersion due to convection is not negligible. The poorly

estimated axial velocity does not have significant effects and the molecular diffusion is negligible. The hydrodynamics model obtained can predict accurately the dispersion for different gas and liquid flowrates in the bench scale reactor and can predict the dispersion in an industrial sized wastewater reactor with the same precision than the previously established empirical correlation (Potier et al., 2005). This correct estimation of the global hydrodynamics allows the future implementation of a kinetic model coupled with the hydrodynamics model. Notation C Cd CL C , C1 , C2 , Ck , C d D Dm g h I I k L Le p Patm Qg Ql Rij

concentration, kg m−3 drag coefficient, dimensionless constant in the Lagrangian particle motion model, dimensionless constants in the k. model or RSM model, dimensionless diameter of the dispersed phase (bubble or particle), m dispersion coefficient, m2 s−1 molecular diffusivity, m2 s−1 gravity acceleration, m s−2 reactor height, m momentum transfer term between phases, kg m−2 s−2 unit tensor, dimensionless turbulent kinetic energy, m2 s−2 reactor length, m characteristic length of turbulent eddies, m pressure, Pa atmospheric pressure, 101 300 Pa gas flowrate, m3 s−1 liquid flowrate, m3 s−1 Reynolds tensor component, m2 s−2

2448

Sct t u U u∞ V w x, y, z

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449

turbulent Schmidt number, dimensionless time, s velocity, m s−1 statistical average velocity, m s−1 terminal vertical velocity, m s−1 reactor volume, m3 reactor width, m spatial coordinates

Greek letters 

ij  l     k    ,  k ct

volume fraction, dimensionless Kronecker factor, dimensionless turbulent dissipation rate, m2 s−3 volume average of turbulent dissipation rate, m2 s−3 normally distributed random number, dimensionless dynamic viscosity, Pa s bubble induced turbulent of dissipation rate source term, kg m−1 s−4 bubble induced turbulent kinetic energy source term, kg m−1 s−3 density, kg m−3 surface tension, N m−1 constants in the k. model or RSM model, dimensionless eddy crossing time, s

e Subscripts g i, j, k l p q t tr

gas spatial coordinates liquid particle phase index: l or g turbulent tracer

Operator az ax,y

averaged value of a on a line along z averaged value of a on a (x, y) plane

eddy lifetime, s

Acknowledgements The authors would like to thank Richard Lainé for his help with the LDV experimental setup, Benoit Fiers for his help in the experimental works and Noël Midoux for his scientific contribution. References Arnold, F.C., 1988. Physical model for two-phase flow in steam injection wells. Conference Paper, American Institute of Chemical Engineers, National Meeting, New York, United States. pp. 42–76. CD-Adapco-Group, Computational Dynamics Limited, Manual – Methodology, STAR-CD version 3.2, 2003. Cockx, A., Do-Quang, Z., Audic, J.M., Liné, A., Roustan, M., 2001. Global and local mass transfer coefficients in waste water treatment process by

computational fluid dynamics. Chemical Engineering and Processing 40, 187–194. Daly, B.J., Harlow, F.H., 1970. Transport equations in turbulence. The Physics of Fluids 15 (11), 2634–2649. Desai, R.B., Kolhatkar, R.V., Joshi, J.B., Ranade, V.V., Malshelkar, R.A., 1995. Turbulence structure in bubble disengagement zone: role of polymer addition. A.I.Ch.E. Journal 41 (5), 1329–1332. Ekambara, K., Dhotre, M.T., Joshi, J.B., 2005. CFD simulations of bubble column reactors: 1D, 2D and 3D approach. Chemical Engineering Science 60, 6733–6746. Elghobashi, S.E., Abou-Arab, T.W., 1982. A two-equation turbulence model for two-phase flows. The Physics of Fluids 26 (4), 931–938. Fluent Inc., 2005. Fluent user’s guide 6.2. Discrete Phase Models 23, 1–14. Glover, G.C., Printemps, C., Essemiani, K., Meinhold, J., 2006. Modeling of wastewater treatment plants—how far shall we go with sophisticated modeling tools? Water Science and Technology 53, 79–89. Gu, H.-Y., Guo, L.-J., 2005. Modelling and simulation of the dynamic flow behaviour in a rectangular bubble column. Journal of Engineering Thermophysics 26, 72–75. Gunjal, P.R., Ranade, V.V., Chaudhari, R.V., 2003. Liquid distribution and RTD in trickle bed reactor: experiments and CFD simulations. Canadian Journal of Chemical Engineering 81, 821–830. Jakobsen, H.A., Sannas, B.H., Grevskott, S., Svendsen, H.F., 1997. Modeling of vertical bubble-driven flows. Industrial and Engineering Chemistry Research 36 (10), 4052–4074. Jakobsen, H.A., Lindborg, H., Dorao, C.A., 2005. Modeling of buble column reactors: progress and limitations. Industrial and Engineering Chemistry Research 44, 5107–5151. Jamialahmadi, M., Branch, C., Muller-Steinhagen, H., 1994. Terminal bubble rise velocity in liquids. Chemical Engineering Research and Design 72, 119–122. Lain, S., Bröder, D., Sommerfeld, M., Göz, M.F., 2002. Modeling hydrodynamics and turbulence in a bubble column using the Euler–Lagrange procedure. International Journal of Multiphase Flow 28, 1381–1407. Lapin, A., Lübbert, A., 1994. Numerical simulation of the dynamics of twophase gas–liquid flows in bubble columns. Chemical Engineering Science 49, 3661–3674. Leon-Becerril, E., Cockx, A., Liné, A., 2002. Effect of bubble deformation on stability and mixing in bubble columns. Chemical Engineering Science 57, 3283–3297. Levin, M.A., Gealt, M.A., 1993. Biotreatment of Industrial and Hazardous Waste. MCGraw-Hill, New York. 71–72. Morsi, S.A., Alexander, A.J., 1972. An investigation of particle trajectories in two-phase flows systems. Journal of Fluids Mechanics 55 (2), 193–208. Mudde, R.F., Simonin, O., 1999. Two- and three-dimensional simulations of a bubble plume using a two-fluid model. Chemical Engineering Science 54, 5061–5069. Olmos, E., Gentric, C., Midoux, N., 2003. Numerical description of flow regime transitions in bubble column reactors by multiple gas phase model. Chemical Engineering Science 58, 2113–2121. Pfleger, D., Becker, S., 2001. Modelling and simulation of the dynamic flow behaviour in a bubble column. Chemical Engineering Science 56, 1737–1747. Potier, O., Leclerc, J.P., Pons, M.N., 2005. Influence of geometrical and operational parameters on the axial dispersion in an aerated channel reactor. Water Research 39, 4454–4462. Ranade, V.V., 2002. Process system engineering. Computational Flow Modelling for Chemical Reactor Engineering, vol. 5. Academic Press, New York, pp. 57–112. Roustan, M., Line, A., 1996. Rôle du brassage dans les procédés biologiques d’épuration. Tribune de l’eau 5–6, 109–115. Sato, Y., Sekoguchi, K., 1975. Liquid velocity distribution in two-phase bubble flow. International Journal of Multiphase Flow 2 (1), 79–95. Sokolichin, A., Eigenberger, G., Lapin, A., 2004. Simulation of buoyancy driven bubbly flow: established simplifications and open questions. A.I.Ch.E. Journal 50, 24. Stropky, D., Pouqatch, K., Nowak, P., Salcudean, M., Pagorla, P., Gartshore, I., Yuan, J.W., 2007. RTD (residence time distribution) predictions in large

Y. Le Moullec et al. / Chemical Engineering Science 63 (2008) 2436 – 2449 mechanically aerated lagoons. Water Science and Technology 55 (11), 29–36. Talvy, S., Cockx, A., Line, A., 2007. Modeling hydrodynamics of gas–liquid airlift reactor. A.I.Ch.E. Journal 53 (2), 335–353. Thyn, J., Ha, J.J., Strasak, P., Zitny, R., 1998. RTD prediction, modelling and measurement of gas flow in reactor. Nukleonika 43 (1), 95–114. Troshko, A., 2006. Best practice for modelling bubble column reactor with FLUENT. www.fluentusers.com. Vial, C., Poncin, S., Wild, G., Midoux, N., 2005. Experimental and theoretical analysis of axial dispersion in the liquid phase in external-loop airlift reactors. Chemical Engineering Science 60, 5945–5954.

2449

Villermaux, J., 1993. Génie de la réaction chimique, conception et fonctionnement des réacteurs. TEC & DOC - Lavoisier, 159–193. Zhang, L., Pan, Q., Rempel, G.L., 2007. Residence time distribution in a multistage agitated contactor with Newtonian fluids: CFD prediction and experimental validation. Industrial and Engineering Chemistry Research 46 (11), 3538–3546. Zhao, Y., Chen, G., Yuan, Q., 2006. Liquid–liquid two-phase flow patterns in rectangular microchannel. A.I.Ch.E. Journal 52, 4052–4060.