PIV measurements and Eulerian–Lagrangian simulations of the unsteady gas–liquid flow in a needle sparger rectangular bubble column

Chemical Engineering Science 126 (2015) 560–572

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

PIV measurements and Eulerian–Lagrangian simulations of the unsteady gas–liquid flow in a needle sparger rectangular bubble column S. Besbes a,n, M. El Hajem b, H. Ben Aissia a, J.Y. Champagne b, J. Jay c a

Metrology Research Unit and Energy Systems (URMSE), National School of Engineers of Monastir, Road Ouerdanine, 5000 Monastir, Tunisia Laboratory of Fluid Mechanics and Acoustics (LMFA), Lyon, France c Thermal Center of Lyon (CETHIL), INSA Lyon, 20 Av. A. Einstein, 69621 Villeurbanne Cedex, France b

H I G H L I G H T S







Experiments with particle image velocimetry and Euler–Lagrange simulations in a needle sparger rectangular bubble column. Effects of gas flow rates on the dynamic and time-averaged flow properties were studied. The simulated vertical liquid velocity was found to be satisfactory to the PIV measurements for low flow rates. The stronger bubble plume velocity oscillations were located near the entrance zone.

art ic l e i nf o

a b s t r a c t

Article history: Received 16 June 2014 Received in revised form 24 October 2014 Accepted 17 December 2014 Available online 24 December 2014

In this work, we conducted a detailed experimental as well as a computational study for the hydrodynamic characterization of the flow in a needle sparger rectangular bubble column. Particle image velocimetry technique (PIV) was used to simultaneously capture the images of both bubbles and seeding polyamide tracers of 50 μm in diameter, the tracer particles were selected to provide sufficient contrast and were neutrally buoyant. A three-dimensional computational model based on an Euler– Lagrange (E–L) approach was performed to simulate the dynamic characteristics of the oscillating bubble plume. The continuous phase velocity field was obtained by solving the unsteady Reynolds averaged Navier–Stokes equations along with the k–ε turbulence model. Bubble tracking was achieved by solving the Newton equations of motion taking into account drag force, pressure, buoyancy and gravity. Twoway coupling between the liquid and gas phases are accounted for in the continuous phase momentum equation. The effect of gas flow rate on the dynamic and time-averaged flow properties were studied, the flow was found to have an unsteady structure with three vortices moving in the upward direction and the bubble plume between them oscillated. The simulated results were found to be in accordance with the PIV measurements for a gas flow rate up to 0.1 l/min. The simulated turbulent kinetic energy indicates that the stronger bubble plume velocity oscillations are located near the entrance zone and are caused by the addition of a shear induced turbulence produced by the presence of an oscillating bubble plume. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Bubble column Hydrodynamics Particle image velocimetry PIV CFD

1. Introduction Reactors, where the movement of the liquid phase is driven by a gas flow, called bubble columns are more and more used in the modern bioprocess research and development, over a broad spectrum of processes ranging from the production of very expensive n Correspondence to: Metrology Research Unit and Energy Systems, Department Energetic of National School of Engineers of Monastir, Road Ouerdanine, 5000 Monastir, Tunisia. Tel.: þ 216 73500826/97664041. E-mail address: [email protected] (S. Besbes).

http://dx.doi.org/10.1016/j.ces.2014.12.046 0009-2509/& 2015 Elsevier Ltd. All rights reserved.

biochemical’s to waste water treatment. The performance of bubble column reactors depends on the hydrodynamic behavior where local flow, turbulence and gas holdup distribution are interrelated in a complex way. It is established that two regimes, homogeneous and heterogeneous, occur in bubble columns. The homogeneous regime occurs at relatively low gas flow rates and is characterized by a uniform size of nearly spherical bubbles. The heterogeneous regime occurs at relatively high gas flow rates and is characterized by a wider bubble size distribution and the presence of large-scale vortical liquid circulation structures forcing an unsteady flow (Chen et al., 1994). However, Mudde et al. (1997a) and Mudde and Simonin (1999)

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

Nomenclature CD dB FD FG FL FP FVM g H k Ml mB P Q Re TL t uB ul u0 , v0 , w0 Ul VB Vcell W X, Y, Z

drag coefficient bubble diameter (m) drag force (N/m3) gravity–buoyancy force (N/m3) lift force (N/m3) pressure force (N/m3) virtual mass force (N/m3) acceleration due to gravity (m/s2) liquid height (m) turbulent kinetic energy of liquid phase (m2/s2) inter-phase momentum transfer (N/m3) bubble mass (kg) pressure (Pa) gas flow rate (l/min) bubble Reynolds number integral time scale of turbulence (s) time (s) bubble velocity (m/s) instantaneous liquid velocity (m/s) fluctuating liquid velocities (m/s) mean liquid velocity (m/s) bubble volume (m3) cell volume (m3) column width (m) Cartesian coordinates (width, height and depth, respectively) (m)

showed that it is possible to find vortical structures even at low superficial gas flow rates. The size of the bubble column, the length scale of hydrodynamic phenomena and the operating mode wherein the bubble column is used determine the most appropriate model to study the problem. Dispersed gas–liquid flow in bubble columns is unsteady and can be simulated using either Euler–Euler or Euler– Lagrange approaches. In the Euler–Euler approach the phases are treated as interpenetrating continua and interact through the source terms in the momentum equation which is calculated as the summation of the various inter-phase forces (e.g. drag, lift and virtual mass forces) (Sokolichin and Eigenberger, 1994; Torvik and Svendsen, 1990; Krishna et al., 2000). In the Euler–Lagrange approach the continuous phase is treated in a usual Eulerian description, whereas the dispersed phase is described in a Lagrangian way, this means that the dispersed phase is tracked through the flow domain from inlet to outlet. Sokolichin et al. (1997) have clearly shown that, provided that appropriate numerical discretization schemes are used, the Euler– Lagrange and Euler–Euler model leads to essentially the same results for the hydrodynamic behavior. The added value for the Euler– Lagrange model is the fact that the positions of all the bubbles are known, this has the great advantage of assessing the interactions between bubbles with a deterministic manner. It is clear that due to its discrete character the Euler–Lagrange model is better suited to model bubble size distributions. However, according to its relative high computational load this model has been limited to the situation where the gas velocity and gas holdup are relatively low. Lapin and Lübbert (1994), Delnoij et al. (1997a,b), Lain et al. (1999, 2002) and Buwa et al. (2006) used the Euler–Lagrange model to simulate the gas liquid flow in a rectangular bubble column. It should be noted that in the majority of these simulations, the liquid phase flow was assumed to be laminar except for Lain et al. (1999, 2002) who carried out two-dimensional turbulent simulations of gas–liquid flow in cylindrical bubble columns. Buwa et al. (2006) used the standard k–ε model to simulate the

Yn

561

dimensionless vertical coordinate

Greek letters

α ΔtE ΔtL Δtn ε λ μ ξ ρ τ τB τe

volume fraction Eulerian time step (s) Lagrangian time step (s) characteristic time (s) dissipation rate (m2/s3) step length factor viscosity (Pa s) normally distributed random number density (kg/m3) viscous stress tensor (Pa) bubble relaxation time (s) eddy life time (s)

Subscripts B eff g l L T

bubble effective gas phase liquid phase laminar turbulent

effective viscosity of liquid phase. Several researchers (such as Pfleger et al., 1999) have reported that two-dimensional simulations over-predict the turbulent viscosity and therefore damp the low frequency oscillations. Buwa et al. (2006) used the Lagrangian approach to simulate the dynamic characteristics of the oscillating bubble plume in a rectangular bubble column. The effect of the superficial gas velocity and aerated liquid height-to-column width (H/W) ratio on the dynamic and time-averaged flow properties was studied. Their results showed that the predicted plume oscillation period was not affected by the lift force for an H/W ratio of 2.25 and a superficial gas velocity up to 0.59 cm/s. For higher H/W ratio, the numerical diffusion and incorrect magnitude of the lift force lead to gas accumulation near the column walls. Since liquid velocity measurements are not available for the highest values of the superficial gas velocity, at which the lift force is considered to have the greatest effect, all quantitative validations with experimental measurements were represented in terms of plume oscillation period and time-averaged gas hold-up. Darmana et al. (2009) used the Euler–Lagrange code (developed by Darmana et al., 2006) to simulate the experimental cases performed by Harteveld et al. (2004) and Harteveld (2005) in a multiple orifice needle sparger rectangular bubble column. Their results showed that large-scale vortical structure was formed in the bottom region of the column when the non-aerated region was located near the side walls, whereas when the non-aerated region was in the center region of the inlet area no large-scale vortical structure was formed. Julia et al. (2007) also investigated the effects of different flow conditions and aeration pattern in a rectangular bubble column with an inlet system similar to those in Harteveld et al. (2003), Harteveld (2005) and Darmana et al. (2009). As the non-aerated regions introduced near the side walls expanded vortical structure was formed near the inlet region. If the non-aerated region was large the vortices grew in the upward direction and the bubble plume between them oscillated.

562

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

However, if the non-aerated region was not large the vortical structure was almost static. While there is an extensive literature based on the Euler–Lagrange simulation models for bubbly flows, it should be noted that most of the simulation studies were realized using two-dimensional geometry, the liquid phase was assumed to be laminar. Studies of a single needle sparger rectangular bubble column at low gas flow rates which are important for validating numerical models are not as common. On the other hand, experimental studies of multiple orifice needle sparger (such as Harteveld et al., 2004; Harteveld, 2005; Mayur et al., 2013) or of porous sparger (such as Ashraf et al., 2011) using the PIV technique were performed at a restricted region of field of view where gas bubbles rise at their terminal velocities covering half of the liquid height. To our knowledge, PIV measurement technique has not been applied to determine the liquid phase velocity field but was rather restricted to find the bubble velocity map or to give a qualitative impression of the large flow structure even at low flow rates (Delnoij et al., 2000; Harteveld et al., 2004; Harteveld, 2005). In the present work, we studied two major purposes, the first one consists of presenting a detailed experimental hydrodynamic characterization of a needle sparger rectangular bubble column operating in homogeneous regime at low gas flow rates. In this way, the influence of the superficial gas velocity on the column hydrodynamics is studied by PIV. In order to provide the liquid phase velocity map over the entire liquid height (0 oyo 500 mm) two images was recorded in two steps each step imaging an area of 200  260 mm2 with an overlapping region of 10 mm. The second purpose of the work focus on the numerical simulations using the Euler–Lagrange approach along with the standard k–ε turbulence model to examine whether or not it can reproduce the characteristics of the liquid phase velocity field induced by the oscillating bubble plume at low gas flow rates.

2. The experimental set-up The experimental study of gas plume hydrodynamics in a 3D bubble column are further complicated by the swirling motion of the bubble plume. A pseudo two-dimensional rectangular bubble column with a large width to depth ratio restricts the degrees of motion of the plume, preventing it from swirling. However, this ratio makes that the wall effect becomes important. In case the lateral dimensions of the rectangular bubble columns in the horizontal directions exceed 3–4 times the initial bubble diameter, the terminal rise velocities were not

affected by the size of the computational domain (Van Sint Annaland et al., 2005). The experimental assembly developed in the LMFA (Fig. 1) consists of a transparent rectangular column of 267 mm(width)  600 mm (height)  15 mm (depth). The column is made of Plexiglas to enable an optical access to the measurement zone. Air and distilled water are used as gas and liquid phases, respectively. The column was initially filled with water up to a height of 500 mm. The gas sparger is a medical needle of 0.4 mm diameter located at the bottom center of the column. A regulator placed between the circuit of compressed air and a rotameter adjusted the air flow rate before it enters the column, three values of gas flow- rates are considered as shown in Table 1. PIV technique was applied to study the instantaneous and mean velocity field of liquid flow structures induced by the rising bubbles. This technique of visualization allows illuminating a plane of the flow by a double pulsed Nd:YAG (532 nm, pulse duration 10 ns, 30 mJ energy per pulse) laser. The liquid is seeded with Polyamide Particles (mean diameter 50 μm and density 1.03 g/cm3). This particle size was found to scatter enough light to be imaged by the CCD camera and ensures that images contain a sufficient number of seeding particles to accurately use cross correlations within the interrogation areas of 32  32 pixel with 25% overlap. The CCD camera is an 8 bit Hi-Sense with a resolution of 1280  1024 pixel placed perpendicular to the laser light sheet it registers double frame images in a sequential way. Since bubbles are much bigger and brighter than tracer particles the intensity of light reflected from the bubble’s surface could saturate the CCD camera. In order to avoid damaging the CCD and to capture simultaneously the images of both bubbles and tracer particles we acts on the camera diaphragm aperture to reduce the laser light intensity scattered from bubbles while ensuring that the tracer particles scatter enough light to be imaged by the CCD camera. Ideally, the use of fluorescent tracer particles and long wave pass filter in front of camera should completely block the laser light reflected by bubbles, however, it was very difficult to detect bubbles since their intensity was very close to background intensity fluctuation (Mayur et al., 2013). A synchronization device was used for the measurements. Time interval between the laser pulses can be adjusted according to the velocity range being measured and is chosen such that the particle displacements are typically less than 1/4 of interrogation region. During this experiment, the pulse separation was 15 ms, 14 ms and 13 ms corresponding respectively to three flow rates (0.05 l/min; 0.1 l/ min; 0.2 l/min). The displacement between laser pulses of particles in

Contol valve Flow meter

Compressed Air

z y

x

Laser Sheet Generation optics 532 nm Pulsed Yag Laser

Water column (267x600x15 mm3) With a needle for gaz injection

8 bits Hi-Sense camera 1280 x 1024

synchronisation and image acquisition unit

PC with FlowManager software for images storage and analysis

Fig. 1. Top view of the experimental setup and PIV system.

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

the flow is analyzed to determine the velocity field over the entire liquid height (0oyo500 mm). Images were recorded in two steps, each step imaging an area of 200  260 mm2 (with 10 mm overlapping region) providing a resolution of 0.2 mm/pixel. The recorded images are analyzed using commercial software “Flow Map” provided by Dantec. After the instantaneous liquid velocity fields are calculated they are validated to identify and remove “false” vectors from the statics calculation.

3. Euler–Lagrange model The homogenous or dispersed bubble regime characterized by relatively low gas velocities and small spherical bubbles that do not coalesce is described by a three-dimensional discrete bubble model. In the Euler–Lagrange model the liquid phase is treated as a continuum by solving the volume averaged Navier–Stokes equations, while the dispersed phase is solved by tracking a large number of bubbles through the calculated flow field using the Newtonian equation of motion. Two-way coupling between the liquid phase and the bubbles are accounted for in the continuous phase momentum equation. 3.1. Liquid-phase hydrodynamics The liquid phase hydrodynamics is represented by volume averaged, incompressible, transient Navier–Stokes equations, which consist of the continuity and momentum equations: ∂ ðα ρ Þ þ ∇  ðαl ρl U l Þ ¼ 0 ∂t l l

ð1Þ

∂ ðα ρ U Þ þ∇  ðαl ρl U l U l Þ ¼  αl ∇P  ∇  ðαl τl Þ þ αl ρl g þ M l ∂t l l l

ð2Þ

where αl is the liquid volume fraction, ρl is the liquid density, Ul is the mean liquid velocity, P is the pressure shared by both gas and liquid phases, τl is the viscous tensor, g is the gravity constant and Ml accounts for the inter-phase momentum transfer from all the gas bubbles to the liquid phase per unit volume and per unit time. The liquid phase viscous stress tensor τl is assumed to obey the general form for a Newtonian fluid: 2 3

τl ¼ αl μef f ;l ð∇U l þ ∇U l T Þ  αl μef f ;l δij ð∇:U l Þ

ð3Þ

where meff,l is the effective viscosity due to both, the molecular viscosity and the turbulent viscosity:

μef f ;l ¼ μL;l þ μT;l

ð4Þ

Table 1 Hydrodynamic parameters and bubble sizes in the experiments. Case number

Flow rate Bubble diameter (l/min) (mm)

Superficial gas velocity (cm/ s)

C1 C2 C3

0.05 0.1 0.2

0.02 0.04 0.08

1.5 2 2.5

563

3.2. Bubble dynamics The bubbles velocity can be computed by solving the Newton second law for each individual bubble: mB

duB ¼ F G þ F P þ F D þ F L þ F VM dt

ð5Þ

where mB is the bubble mass, uB is the bubble velocity and FG, FP, FD, FL and FVM, are the forces due to gravity and buoyancy, pressure gradient, drag, lift and virtual mass force, respectively. The bubble trajectories can be computed from the bubble velocities as dxB ¼ uB dt

ð6Þ

In many applications, for instance in the case of dilute flows with spherical particles, it is possible to simplify the general particle momentum equation and consider only the gravity–buoyancy, pressure and drag forces. In comparison with the drag force, the magnitude of the other two inter-phase forces, lift and virtual mass are small. Rampure et al. (2003) and Diaz et al. (2008) indicate that inclusion of virtual mass force does not result in any significant difference in dynamic and time-averaged flow properties. Delnoij et al. (1997a) found that the inclusion of lift force with a lift coefficient of 0.53 led to a higher dispersion of bubbles in the upper region of the column and a breakdown of the re-circulatory flow. Bhole et al. (2008) and Silva et al. (2012) found that at low and high superficial gas velocities the inclusion of lift effects in the model resulted in worse predictions for the gas holdup. On the other hand, Zhang et al. (2013) showed that consideration of virtual mass force and lift force improved prediction results of gas hold-up in a region which was much closer to the sparger, but there was not improvement for the velocity prediction. Tabib et al. (2008) was found that using a lift model at low superficial gas velocity has not significant effect on the flow pattern in comparison with high superficial velocities of different bubble sizes. However, using the proper lift model corresponding to the bubble size significantly improves the numerical result. Recently, Pourtousi et al. (2014) found that the lift model improves the accuracy in prediction of the axial liquid velocity, turbulent kinetic energy and radial gas hold-up. The sum of forces due to gravity–buoyancy and pressure gradient in the liquid phase is calculated as
 ρ mB F G þ F P ¼ mB g 1  l þ ∇P ð7Þ

ρB

ρB

The drag force exerted by the liquid on the bubble is given by: FD ¼

τB ¼

mB

τB

ðul  uB Þ

ð8Þ

4ρB dB 1 ρ dB jul uB j where Re ¼ l 3μl ReC D μl 2

ð9Þ

here, ul is the instantaneous liquid velocity, ρB is the bubble density, dB the bubble diameter, ml the liquid viscosity and g the constant of gravity. τB represents the relaxation time of the bubble, Re the Reynolds number based on the relative velocity and CD the drag coefficient which can be taken from Morsi and Alexander (1972):

Table 2 Boundary conditions for the continuous phase model (CPM) and the discrete phase model (DPM). Boundary conditions Needle orifice Top surface Wall

Velocity inlet Surface injection Wall Wall

CPM

DPM

Velocity, turbulent intensity and hydraulic diameter No setting Stationary wall, no slip Stationary wall, no slip

Escape Bubble diameter, velocity magnitude and flow rate (Table 1) Escape Reflect

564

C D ¼ a1 þ

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

a2 a3 þ Re Re2

ð10Þ

where a1, a2 and a3 are constants that apply to smooth spherical particles over several ranges of Reynolds number given by Morsi and Alexander (1972): For 00 o Reo 1000;

a1 ¼ 0:3664;

a2 ¼ 98:33

and

a3 ¼  2778:

The instantaneous liquid velocity ul at the bubble location accuring in Eq. (8) is estimated as the sum of the mean liquid velocity Ul and a fluctuating velocity component u0 as ul ¼ U l þ u0 ðtÞ

ð11Þ

For the estimation of the fluctuating liquid velocity we used the discrete random walk (DRW) model, which is also called the “eddy lifetime model”, the fluctuating velocity is added to the mean liquid velocity for a constant time interval given by the characteristic lifetime of liquid eddies. The fluctuating velocity components u0 , v0 and w0 are assumed to have a Gaussian distribution during the lifetime of an eddy and are estimated as pffiffiffiffiffiffi u0 ¼ ξ u02

ð12Þ

pffiffiffiffiffiffi where ξ is a normally distributed random number and u02 is the root mean square (RMS) of local velocity fluctuations. Since the local turbulent kinetic energy k is known, the RMS fluctuating

Fig. 2. PIV images of oscillating bubble plume, (a) Q¼ 0.05 l/min, (b) Q¼ 0.1 l/min, (c and d) Q ¼0.2 l/min.

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

3.4. Numerical simulations

components are estimated from: pffiffiffiffiffiffi pffiffiffiffiffiffi pffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffi u02 ¼ v02 ¼ w02 ¼ 2k=3

ð13Þ

The local turbulent kinetic energy k and its dissipation rate ε are obtained using the standard k–ε turbulence model. The characteristic eddy life time is defined as a random variation about TL:

τe ¼  T L log ðrÞ

ð14Þ

where r is a random number from 0 to 1, TL is the integral time scale of turbulence and is estimated from the turbulent kinetic energy k and its dissipation rate ε: T L ¼ 0:15

k

ð15Þ

ε

565

Over this eddy lifetime τe, the bubbles are assumed to interact with the liquid phase. When the eddy life time is reached, a new value of the instantaneous velocity was obtained by applying a new value of ξ in Eq. (12).

3.3. Boundary and initial conditions The simulations of the rectangular bubble column were carried out in a transient manner to catch the highly dynamic flow structures in the apparatus. The dimensions of the simulated column are the same as the experimental one (see Fig. 1) except that the simulated height is the clear liquid height, i.e. 500 mm. The effect of mesh size on the simulation results was investigated and a uniform grid of 31(width)  100(height)  5 (depth) cells ensures that all simulations were grid independent. Starting at t¼0.0 s gas was fed to the column through a needle of 0.4 mm diameter, numerically the needle orifice was modeled as a velocity inlet. It has been assumed that all the bubbles are spherical, their average diameter and rise velocity at the inlet were deduced from the experimental conditions (see Table 1). No slip boundary conditions were used at all the impermeable walls for the continuous phase, where bubbles are allowed to reflect off these impermeable walls. The top surface of the column was specified as a stationary wall boundary condition, which implies a no-slip condition for the continuous phase and escape boundary condition for the discrete phase, with the assumption that gas bubbles leave the computational domain at the top surface. The settings used for the boundary conditions are summarized in Table 2, for the initial conditions the gas phase velocity at the needle orifice was set from Table 1.

The resulting set of Equations for the liquid phase was solved using the commercial solver Fluent based on a finite volume method. The pressure implicit with splitting of operator (PISO) algorithm was used for the pressure–velocity coupling in the momentum equation. The time derivatives were descretized using a first order implicit method, while the diffusive and convective terms for velocities, kinetic energy (k) and its dissipation rate (ε ) were descretized using the higher-order Quick discretization scheme. It should be noted that the transient simulations need to be carried out for several plume oscillation cycles in order to obtain a reasonably well-developed time-averaged flow after discarding initial transients. An Eulerian time step of 0.01 s was used in all simulations. The equation of bubble motion was analytically integrated by assuming that the forces such as gravity–buoyancy, pressure and drag are constant during the time step. The resulting equation was descretized and solved using the implicit and trapezoidal numerical schemes, in combination with automated tracking scheme selection. The numerical solution requires that the time step of integration (Lagrangian time step ΔtL) has to be smaller than all relevant time scales for the bubble motion, namely the time required for a bubble to traverse the continuous phase control volume (Δtn), the bubble response time scale (τB) and the integral time scale of turbulence (TL). Therefore, the Lagrangian time step ΔtL was computed based on the characteristic time Δtn as

Δt L ¼

Δt n λ

ð16Þ

where λ is the step length factor, which is roughly equal to the number of time steps required to traverse the continuous phase control volume. A step length factor of 100 was used for all simulations. The coupling between the liquid phase and the bubbles (two-way coupling) was accomplished through both the liquid volume fraction αl and the inter-phase momentum transfer Ml from bubbles to liquid calculated in accordance with the number of bubbles present in a computational cell. The local gas-phase volume fraction in a computational cell was derived from the volume occupied by the bubbles in that cell and the total volume of that cell. P VB αg ¼ n ð17Þ V cell where n is the number of bubbles. The momentum transferred from the bubbles to the liquid was estimated for each bubble as the sum of instantaneous momentum contribution along the bubble trajectory in

Fig. 3. Close-up PIV images of oscillating bubble plume, (a) Q¼ 0.05 l/min, (b) Q¼ 0.1 l/min and (c) Q ¼ 0.2 l/min.

566

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

Fig. 4. Time-averaged liquid velocity fields from PIV measurements, (a) Q¼ 0.05 l/min, (b) Q¼ 0.1 l/min, (c) Q ¼ 0.2 l/min.

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

each control volume under the Lagrangian time step Ml ¼ 

XX 1 ðF D ÞΔt L V cell Δt E n k

ΔtL as ð18Þ

where Vcell is the cell volume, ΔtE and ΔtL are Eulerian and Lagrangian time steps, respectively, FD is the momentum transferred through drag

567

force. The summation over k represents the sum of the instantaneous momentum contributions along the bubble trajectory in the control volume and the summation over n accounts for the number of bubbles passing through the control volume under consideration. Once, the liquid-phase velocity field is calculated, the bubble velocities and positions are computed; thereafter the forces acting on the bubbles are evaluated and transferred into the momentum

Fig. 5. Snapshots of PIV and E–L simulation of oscillating bubble plume and liquid velocity fields for Q ¼0.2 l/min.

568

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

equation for the liquid phase. The new liquid velocity field is then determined. Typically an average about 100 to 150 bubbles are present in the domain for the range of superficial gas velocities considered. This amount of bubbles is resulting from the balance between bubbles coming into the column through the needle and the bubbles leaving the column from the top surface. During the coupled calculation the bubble positions are updated every each iteration of the continuous phase. The determination of the time averaged liquid and bubble phase properties begins when the flow pattern is established and a single recirculation loop exists in the bubble column. The averaging time for the properties of both phases is about 210 s.

gradient at the edges of the bubble and bubbles diameter varies according to the flow rate considered (see Table 1). For the measurement of the liquid velocity with PIV it is assumed that the tracer particles follow accurately the liquid phase. The time averaged liquid velocity field obtained from PIV measurements along the entire liquid height (0 oYn o1), as represented in Fig. 4, shows that the liquid is entrained upward in the wake of the bubbles and downwards near the column walls. This liquid circulation mode commonly referred as “Gulf-stream” or “cooling tower” flow regime (Chen et al., 1989) differs considerably from the instantaneous flow regime experimentally observed. It is therefore essential to simulate the flow in bubble columns in a transient manner to capture the oscillations of the

4. Results and discussion Experimental results were first examined for the three cases considered (see Table 1), at low gas flow rates the column operates in the dispersed (homogenous) bubbly regime characterized by the absence of bubble coalescence or breakup. Instantaneous PIV images of the oscillating bubble plume for the lower area (0oYn o0.52) of the column are shown in Fig. 2, where Yn is a dimensionless axial coordinate with respect to the liquid height (Yn ¼Y/H). When the flow rate is less than 0.2 l/min (Fig. 2(a and b)) a little quantity of spherical gas bubbles is rising from the needle and the bubble plume oscillations are low, for a flow rate of 0.2 l/min (Fig. 2(c and d)) the number of bubbles rising from the needle increases significantly, the direction of the lower part of the bubble plume is stable however the upper part changes its appearance and location corresponding to the transient liquid circulation flows. Close-up PIV images of oscillating bubble plume are shown in Fig. 3. As can be seen from this figure the bubble size is relatively small and its distribution is narrow, the bubble contour is not well captured due to an increased brightness level

Fig. 7. E–L simulation of time-averaged vertical liquid velocity profiles for Q¼ 0.2 l/ min at two different times.

Fig. 6. E–L simulation for Q¼ 0.2 l/min. (a) mean liquid velocity field, (b) corresponding bubble plume.

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

569

bubble plume and to accurately describe the hydrodynamics as well as mixing and transfer process in bubble columns. Typical snapshots of oscillating bubble plume and liquid velocity field obtained from Eulerian–Lagrangian simulations at gas flow rate of 0.2 l/min are shown in Fig. 5(c–f), respectively. As can be seen from these figures the oscillations of the bubble plume are captured in satisfactory qualitative agreement with the experiment (Fig. 5(a and b)), the form of undulation of the bubble plume is found to differ slightly from that observed in experiment however, the stable lower part of the bubble plume is correctly reproduced in the simulation. In Fig. 5(e and f) the flow is found to have an unsteady structure with staggered rows of vortices moving downwards in a periodic way. As found by Borchers et al. (1999) and Diaz et al. (2006) the number of vortices depends on the aspect ratio (H/W) and increases from two at H/W¼1.5 to three at H/W¼ 2 and four at H/W¼3. In our case the aspect ratio (H/W) is 1.87 and three vortices are observed however, the third vortex which is located in the top left (Fig. 5(e)) or in the top right (Fig. 5(f)) region of the column is not well developed. The time averaged liquid velocity field and the bubble plume are represented for a flow rate of 0.2 l/min in Fig. 6(a and b), respectively. This time averaged flow pattern resembles the classical “cooling tower” flow pattern with liquid up flow in the center of the column and liquid down flow near the walls of the column as shown in Fig. 6(a). The evolution of the axial mean liquid velocity is analyzed for a flow rate of 0.2 l/min at the mid-depth plane of the column at two positions Yn ¼ 0.25 and 0.75 above the needle and at two different times. As shown in Fig. 7 the variations of the mean liquid velocity after 210 s are very small, the flow never reaches steady-state but time

averaging produces reproducible stationary patterns obtained at 210 s and 606 s. Quantitative comparison between the measured and the simulated vertical and horizontal time-averaged liquid velocity profiles at the mid-depth plane of the column are represented respectively in Figs. 8(a and b) and 9(a and b) at two positions Yn ¼0.25 and 0.75 above the needle. The results from Fig. 8(a and b) indicate that the numerical and experimental vertical velocity profiles exhibit a reasonable agreement when the flow rate was increased up to 0.1 l/min illustrating as well the strong up flow at the column center and the down flow near the column walls. Furthermore, when the flow rate was increased the structure of the liquid velocity profiles was preserved. However, a flow rate of 0.2 l/min causes an enlargement of the rising liquid area due to a greater quantity of bubbles which was not well reproduced by the simulation and the discrepancy is clearly observed in Fig. 8(b) at Yn ¼0.75 where the simulated velocity gradient shows a steep increase to the velocity maximum not observed in PIV measurements. It should be noted that the final averaged velocity fields obtained from PIV measurements were performed from 2000 images for a flow rate of 0.05 l/min; and from 1000 images for a flow rate of 0.1 and 0.2 l/ min, respectively. Furthermore, the slight asymmetry in the velocity profiles for the flow rates of 0.1 and 0.2 l/min as seen in Fig. 8 is a result of an insufficient number of acquisition images. As shown in Fig. 9(a) the simulated horizontal liquid velocity profiles exhibit a reasonable agreement with PIV measurements for the lower region of the column at Yn ¼0.25. However, at Yn ¼0.75 (Fig. 9(b)), the discrepancy between simulated and experimental

Fig. 8. PIV and E–L simulation of time-averaged vertical liquid velocity profiles at two positions above the needle (a) Y* ¼ 0.25 and (b) Y* ¼ 0.75.

Fig. 9. PIV and E–L simulation of time-averaged horizontal liquid velocity profiles at two positions above the needle (a) Y* ¼0.25 and (b) Y* ¼ 0.75.

570

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

Fig. 10. E–L simulation of time-averaged turbulent kinetic energy profiles according to the axial direction at the mid-depth plane of the column.

horizontal liquid velocity is clearly observed for the flow rate of 0.2 l/min. To study the location and strength of the bubble plume oscillations we are represented the time averaged turbulent kinetic energy for various column heights and air flow rates at the mid-depth plane of the column. As observed in Fig. 10 the profiles of the turbulent kinetic energy show a peak located at a closer and closer distance of the entrance zone of the column as the flow rate increases. The higher turbulent kinetic energy value about 0.07 m2/s2 was obtained for Q¼0.2 l/min at Yn ¼0.20. This higher value indicates that the stronger bubble plume velocity oscillations are located near the entrance zone and are caused by the addition of a shear induced turbulence produced by the presence of an oscillating bubble plume. Simulated time series of the vertical liquid velocity over a measurement time of 100 s are shown in Fig. 11 at different gas flow rates at point (X¼0, Yn ¼0.5, Z¼0). High-frequency turbulent fluctuations around the mean velocity value can be observed for the three flow rates considered, in addition low-frequency oscillations of the liquid velocity corresponding to the periodic movement of the bubble

Fig. 11. E–L simulation of time series of vertical liquid velocities and corresponding spectral density at different gas flow rates, at point (X ¼0, Yn ¼0.5, Z¼ 0).

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

plume can be seen only at Q¼ 0.2 l/min. To calculate the number of low-frequency oscillations the data set of the vertical liquid velocity is transformed from the time domain to the frequency domain by means of the calculation of the fast Fourier transform. The bubble plume oscillation frequency for different gas flow rates in the rectangular needle sparger column are shown in Fig. 11. When the flow rate is up to 0.1 l/min the oscillations of the bubble plume are weak however, for Q¼0.2 l/min the frequency of oscillations is about 0.125 Hz. This frequency value is equal to that obtained by Julia et al. (2007) in a multiple needle sparger rectangular bubble column for the F16 aeration pattern which consists of 25 needles (inner diameters of 0.8 mm located at the bottom center of the column) and corresponds to the free bubble plume picture where the oscillations are not affected by the column walls.

5. Conclusion The hydrodynamic characterization of the gas–liquid flow in a needle sparger rectangular bubble column operated in homogenous regime at low flow rates was investigated successfully using PIV measurements and E–L simulations. PIV technique was applied to visualize simultaneously bubbles and seeding polyamide tracers of 50 μm in diameter, this particle size was found to scatter enough light to be imaged by the CCD camera. A fully three-dimensional simulation using the standard k–ε turbulence model is necessary to capture the transient flow structure in a sufficient fine grid. The flow was found to have an unsteady structure with an oscillatory movement of bubble plume that causes the formation of liquid vortices of changing size and position and three vortices moving in the upward direction were obtained for Q¼0.2 l/min. The predicted time-averaged vertical liquid velocity profiles with E–L simulations were found to be satisfactory to the PIV measurements for a flow rate up to 0.1 l/min, however for the flow rate of 0.2 l/min the concordance with experimental results can be improved by adding to the inter-phase momentum transfer the lift force which acts on the dispersion of bubbles over the cross section mainly in the upper region of the column. The evolution of the turbulent kinetic energy according to the axial direction at the middepth plane indicates that the stronger bubble plume velocity oscillations are located near the entrance zone and are caused by the addition of a shear induced turbulence produced by the presence of an oscillating bubble plume. In the absence of reaction, mass transfer is a relatively important factor in applications such as water treatment in aquaculture and plays an important role in the performance of the airlift. One of the authors has already investigated experimentally the dissolution of CO2 when it is this gas that is injected (Valiorgue et al., 2013). Exploring the evolution of the concentration of CO2 near the surface of the bubble is described and mass transfer coefficient is determined (Souzy, 2014). The model established in the present study should be completed with a resolution of the transport equations for predicting these transfer coefficients and test the accuracy of their estimation compared to experimental ones.

Acknowledgment If this work was done, it is thanks to international mobility program CMIRA supported by the French Rhône-Alpes region for which the authors would address all their thanks. References Ashraf, A.B., Pushpavanam, S., 2011. Analysis of unsteady gas–liquid flows in a rectangular tank: Comparison of Euler– Eulerian and Euler– Lagrangian simulations. Original Research Article, Int. J. Multiphase Flow 37, 268–277.

571

Bhole, M., Joshi, J., Ramkrishna, D., 2008. CFD simulation of bubble columns incorporating population balance modeling. Chem. Eng. Sci. 63, 2267–2282. Borchers, O., Busch, C., Sokolichin, A., Eigenberger, G., 1999. Applicability of the standard k–ε turbulence model to the dynamic simulation of bubble columns. Part II: Comparison of detailed experiments and flow simulation. Chem. Eng. Sci. 54, 5927–5935. Buwa, V.V., Deo, D.S., Ranade, V.V., 2006. Eulerian–Lagrangian simulations of unsteady gas–liquid flows in bubble columns. Int. J. Multiphase Flow 32, 864–885. Chen, J.J.J., Jamialahmadi, M., Li, S.M., 1989. Effect of liquid depth on circulation in bubble columns: a visual study. Chem. Eng. Res. Des. 67, 203–207. Chen, R.C., Reese, J., Fan, L.S., 1994. Flow structure in a three-dimensional bubble column and three phase fluidized bed. AIChE J. 40, 1093–1104. Darmana, D., Deen, N.G., Kuipers, J.A.M., 2006. Parallelization of an Euler–Lagrange model using mixed domain decomposition and mirror domain technique: application to dispersed gas–liquid two-phase flow. J. Comput. Phys. 220, 216–248. Darmana, D., Deen, N.G., Kuipers, J.A.M., Harteveld, W.K., Mudde, R.F., 2009. Numerical study of homogeneous bubbly flow: influence of the inlet conditions to the hydrodynamic behavior. Int. J. Multiphase Flow 35, 1077–1099. Delnoij, E., Lammers, F.A., Kuipers, J.A.M., van Swaaij, W.P.M., 1997a. Dynamic simulation of dispersed gas–liquid two-phase flow using a discrete bubble model. Chem. Eng. Sci. 52, 1429–1458. Delnoij, E., Lammers, F.A., Kuipers, J.A.M., van Swaaij, W.P.M., 1997b. Dynamic simulation of gas–liquid two-phase flow: effect of column aspect ratio on the flow structure. Chem. Eng. Sci. 52, 3759–3772. Delnoij, E., Kuipers, J.A.M., van Swaaij, W.P.M., 2000. Measurement of gas–liquid two-phase flow in bubble columns using ensemble correlation PIV. Chem. Eng. Sci. 55 (17), 3385–3395. Diaz, M.E., Montes, F.J., Galan, M.A., 2006. Influence of aspect ratio and superficial gas velocity on the evolution of unsteady flow structures and flow transitions in a rectangular two-dimensional bubble column. Ind. Eng. Chem. Res. 45, 7301–7312. Diaz, M.E., Montes, F.J., Galan, M.A., 2008. Numerical simulation of the gas–liquid flow in a laboratory scale bubble column: influence of bubble size distribution and non-drag forces. Chem. Eng. J. 139, 363–379. Harteveld, W.K., Mudde, R.F., vandenAkker, H.E.A., 2003. Dynamics of a bubble column: influence of gas distribution on coherent structures. Can. J. Chem. Eng. 81, 389–394. Harteveld, W.K., Julia, J.E., Mudde, R.F., van den Akker, H.E.A., 2004. Large scale vortical structures in bubble columns for gas fractions in the range of 5–25%. In: Proceedings of CHISA 16th International Congress of Chemical and Process Engineering, Prague, Czech Republic. Harteveld, W.K., 2005. Bubble Columns-Structure or Stability? (Ph.D. Thesis). Delft University of Technology. Julia, J.E., Hernandez, L., Chiva, S., Vela, A., 2007. Hydrodynamic characterization of a needle sparger rectangular bubble column: homogeneous flow, static bubble plume and oscillating bubble. Chem. Eng. Sci. 62, 6361–6377. Krishna, R., van Baten, J.M., Urseanu, M.I., 2000. Three-phase Eulerian simulations of bubble column reactors operating in the churn-turbulent regime: a scale up strategy. Chem. Eng. Sci. 55, 3275–3286. Lain, S., Broder, D., Sommerfeld, M., 1999. Experimental and numerical studies of the hydrodynamics in a bubble column. Chem. Eng. Sci. 54, 4913–4920. Lain, S., Broder, D., Sommerfeld, M., Goz, M.F., 2002. Modelling hydrodynamics and turbulence in a bubble column using the Euler–Lagrange procedure. Int. J. Multiphase Flow 28, 1381–1407. Lapin, A., Lübbert, A., 1994. Numerical simulations of the dynamics of two-phase gas–liquid flows in bubble columns. Chem. Eng. Sci. 49, 3661–3674. Mayur, S., Jyeshtharaj, J., Geoffrey, E., 2013. Characterization of turbulence in rectangular bubble column. Chem. Eng. Sci. 100, 52–68. Morsi, S.A., Alexander, A.J., 1972. An investigation of particle trajectories in twophase systems. J. Fluid Mech. 55, 193–208. Mudde, R.F., Groen, J.S., van den Akker, H.E.A., 1997a. Liquid velocity field in a bubble column: LDA experiments. In: Chem. Eng. Sci. Proceedings of the Third International Conference on Gas–Liquid–Solid Reaction Engineering, vol. 52, pp. 4217-4224. Mudde, R.F., Simonin, O., 1999. Two and three dimensional simulation of a bubble plume using a two fluid model. Chem. Eng. Sci. 54, 5061. Pfleger, D., Gomes, S., Gilbert, N., Wagner, H.G., 1999. Hydrodynamic simulation of laboratory scale bubble columns: fundamental studies of Eulerian–Eulerian modeling approach. Chem. Eng. Sci. 54, 5091–5099. Pourtousi, M., Sahub, J.N., Ganesan, P., 2014. Effect of interfacial forces and turbulence models on predicting flow pattern inside the bubble column. Chem. Eng. Process.: Process Intensif. 75, 38–47. Rampure, M.R., Buwa, V.V., Ranade, V.V., 2003. Modelling of gas–liquid/gas–liquid– solid flows in bubble columns: experiments and CFD simulations. Can. J. Chem. Eng. 81, 692–706. Silva, M.D., D’Ávila, M., Mori, M., 2012. Study of the interfacial forces and turbulence models in a bubble column. Comput. Chem. Eng. 44, 34–44. Sokolichin, A., Eigenberger, G., 1994. Gas–liquid flow in bubble columns and loop reactors: Part I. Detailed modelling and numerical simulation. Chem. Eng. Sci. 49, 5735–5746. Sokolichin, A., Eigenberger, G., Lapin, A., Lubbert, A., 1997. Dynamic numerical simulation of gas–liquid two-phase flows Euler/Euler versus Euler/Lagrange. Chem. Eng. Sci. 52 (4), 61–626.

572

S. Besbes et al. / Chemical Engineering Science 126 (2015) 560–572

Tabib, M.V., Roy, S.A., Joshi, J.B., 2008. CFD simulation of bubble column an analysis of interphase forces and turbulence models. Chem. Eng. J. 139, 589–614. Souzy, N., 2014. Thesis 14 October 2014, INSA. Torvik, R., Svendsen, H.F., 1990. Modeling of slurry reactors: a fundamental approach. Chem. Eng. Sci. 45, 2325. Valiorgue, P., Souzy, N., El Hajem, M., Ben Hadid, H., Simoëns, S., 2013. Concentration measurement in the wake of a free rising bubble using planar laser-

induced fluorescence (PLIF) with a calibration taking into account fluorescence extinction variations. Exp. Fluids 54, 4. Van Sint Annaland, M., Deen, N.G., Kuipers, J.A.M., 2005. Numerical simulation of gas bubbles behavior using a three-dimensional volume of fluid method. Chem. Eng. Sci. 60, 2999–3011. Zhang, Y., Bai, Y., Wang, H., 2013. CFD analysis of inter-phase forces in a bubble stirred vessel. Chem. Eng. Res. Des. 91 (1), 29–35.