ARTICLES

Dynamic three-dimensional simulation of gas-liquid flow in a cylindrical bubble column

S. Laín¹

¹ Fluid Mechanics Research Group, Energetics and Mechanics Department.
Universidad Autónoma de Occidente, Calle 25 No 115-85, Cali, Colombia
slain@uao.edu.co

Abstract- In this paper the transient threedimensional flow developing in a cylindrical laboratory bubble column is addressed from a numerical point of view. The simulation scheme combines a Large Eddy Simulation (LES) for describing the liquid phase and a Lagragian approach for the gas phase (discrete bubbles) phase. The bubble equation of motion considers all the relevant forces, i.e., buoyancy, fluid stresses, drag, added mass and transverse lift. From the calculations, the transverse lift in combination with the drag is identified as the main mechanism allowing the bubbles to spread over the column cross-section. The liquid and gas velocity profiles obtained are compared with experimental data and k - ε results. As a matter of fact, the dynamic structure of the liquid flow induced by the rising bubbles is well reproduced and also good quantitative results for all measured variables of both phases, gas and liquid, are obtained.

Keywords- Bubble Column; Numerical Simulation; Euler-Lagrange Approach; Turbulence.

I. INTRODUCTION

Bubble columns are frequently encountered in the chemical, petrochemical and biotechnological industries. For the sound design of this kind of devices in process engineering it is necessary to understand their fundamental hydrodynamic behaviour, which is determined by bubble rise, bubble-bubble and bubbleliquid interactions, bubble size and bubble size distribution, and gas hold-up. Moreover, fluctuations of velocity are induced in the liquid by the movement of the bubbles due to the shear produced in the vicinity of the bubbles, in particular due to bubble oscillations. However, while the time-averaged flow within a bubble column shows a very regular and symmetric structure, the transient flow behaviour is generally highly irregular and asymmetric. As the bubbles react to local and instantaneous flow patterns, the dynamic interactions among the bubbles and between bubbles and liquid affect the performance of the column.

From the engineering point of view two approaches are mainly used to simulate flow in bubble columns and in multiphase flow in general (for a review the reader is referred for instance to Jakobsen et al., 1997). Euler-Euler approaches or two-fluid models treat the liquid and gas as two interpenetrating media, where conservation equations are required for each phase together with interphase exchange terms; the main advantage of this method is the relatively low computational demands since two sets of conservation equations with similar structure are solved. However, sophisticated closures are required to describe bubblefluid interaction as well as interactions between bubbles. For polydisperse bubbles (which is usually the case) moreover several sets of conservation equations have to be solved, which will remarkably increase computational demands.

The Euler-Lagrange procedure, on the other hand, treats the continuous phase in an Eulerian frame of reference and the dispersed phase (i.e. particles, droplets and bubbles) is simulated by tracking a large number of particles through the flow field by accounting for all relevant forces acting on the particle. The influence of the dispersed phase on the continuous phase is accounted for by source terms in the continuous phase conservation equations describing the momentum transfer. These source terms and the particle phase properties are obtained by ensamble averaging in the Lagrangian tracking procedure. For obtaining a converged soluting of the coupled system Eulerian and Lagrangian computations have to be performed successively. Therefore, this method is computationally more expensive than Euler-Euler approach but it has the advantage that bubble-bubble interactions and polydispersity are taken into account in a natural manner. In both strategies, the liquid phase is treated as a continuum so the bubble induced velocity fluctuations in the liquid have to be modeled somehow. Generally the computations have been performed in a time-dependent way in order to resolve the inherent unsteady flow structure in bubble columns. These approaches in connection with an appropriate turbulence model are called URANS (unsteady Reynolds-averaged Navier-Stokes). During the last years the k - ε model has been applied by various authors (e.g., Sanyal et al., 1999; Laín et al., 2002; Bourtloutski and Sommerfeld, 2002) to describe the fluctuating structure of the liquid velocity field getting promisingly good qualitative and quantitative agreement with experimental data obtained by PDA (Phase-Doppler Anemometry), PIV-PTV (Particle Image Velocimetry -Particle Tracking Velocimetry) or CARPT (Computer Assisted Radioactive Particle Tracking) techniques. Nevertheless, because of the low Reynolds number of the liquid flow, a Large Eddy Simulation (LES) description of the underlying liquid motion is conceptually more attractive than RANS turbulence models, originally developed for turbulent, high Reynolds number flows.

In fact, the use of LES for vertical bubble driven flows has been suggested since several years ago (e.g., Jakobsen et al., 1997). However, only recently bubble column simulations employing LES formulations have appeared in the literature under the Euler-Euler framework (Deen et al., 2001; Milelli et al., 2001) and using the Euler-Lagrange approach (van den Hengel et al., 2005; Hu, 2005). The main conclusions of these studies were: first, relatively coarse grids could be used without losing any fundamental characteristics of the flow; second, the effect of bubbles on the liquid effective viscosity, the so-called bubble induced turbulence, could be neglected and, finally, the transversal dispersion of the bubbles was mainly driven by the transverse lift force.

This paper focuses on transient simulations of the three-dimensional flow developing in a cylindrical laboratory bubble column. The liquid flow is described by means of Large Eddy Simulation whilst for the bubble phase a Lagrangian tracking is employed, i.e., an Euler-Lagrange scheme is adopted. The main difference to previously performed LES is that the bubbles are present in all the liquid domain instead of being a bubble plume. The results of the simulations are compared with PIV experimental measurements (Bröder and Sommerfeld, 2002) and k - ε calculations (Laín et al., 2001).

II. MODELLING AND NUMERICAL APPROACH

The three-dimensional dynamic simulations of the flow evolving in the bubble column have been performed using the Euler/Lagrange approach. The liquid flow is calculated by using the time-dependent filtered Navier-Stokes equations, i.e., Large Eddy Simulation (LES). Therefore, the three-dimensional equations of continuity and momentum extended by accounting for the effects of the dispersed (gas) phase are considered. They can be written in tensorial notation as (where the comma followed by a subscript means partial derivative and summation is performed over repeated indexes):

	(1)
	(2)

where the hat ( ˆ ) indicates filtered variables. Here, is the liquid density, u is the liquid velocity which is decomposed in a resolved part ú_j and the sub-grid SGS scale (SGS) part , p is the liquid pressure and gj is the respective component of the gravity acceleration. In Equation (2), µ_{ef f} is the effective viscosity for the liquid which is composed of two contributions: the molecular viscosity µ_l and the so-called sub-grid scale turbulent viscosity µ_t:

µ_{ef f} = µ_l + µ_t (3)

The turbulent viscosity accounts for the contribution of the subgrid scales and in this work it is described by the Smagorinsky model (Smagorinsky, 1963), which is essentially an eddy viscosity model. Therefore, this sub-grid scale viscosity is expressed as:

(4)

where is the symmetric deviatoric part of strain tensor of the resolved scales. It is written as:

(5)

In Eq. (4) C_s is the Smagorinsky constant whose value ranges in the literature from 0.065 up to 0.2, and Δ is the filter width. The last term in Eq. (2) FiB, the momentum interaction term, represents the action of the bubble phase on the liquid phase through the interface forces which will be discussed below. Since at present only low void fractions are considered, lower than 2% in all the considered cases, the liquid density is assumed to be unaffected by the presence of the bubbles. Therefore, hereafter we write .

Time integration of the Eulerian equations is performed by a second-order full implicit method in order to decrease the truncation error in time, whilst a central differencing scheme is used for spatial discretisation of the liquid phase equations.

The simulation of the bubble phase by the Lagrangian method requires the solution of the equations of motion for each computational bubble (representing a parcel of a number of real bubbles with identical properties). Hence, the bubble motion after injection is calculated by solving the following set of ordinary differential equations:

	(6)
	(7)

with : inertia force acting on the bubble due to its acceleration; : force due to drag; : force due to gravity and buoyancy; : force due to fluid stresses resulting from the pressure term (Crowe et al., 1998); : force due to virtual mass effect; : force due to transverse lift.

Other forces such as the Basset history term are assumed to be negligible. Here, x_Bi are the coordinates of the bubble position, u_Bi are the velocity components, D_B is the bubble diameter, and is the gas density which is assumed to be constant at present. The symbol means the derivative following the fluid element. are the components of the Levi-Civita pseudo-tensor, which are used to express the curl or the cross product of vectors.

The drag coefficient C_D is calculated using the empirical correlations for a fluid bubble (Lain et al., 2002),

(8)

where is the bubble Reynolds number.

The transverse lift coefficient C_TL appearing in the bubble equation of motion (7) is based on the correlation proposed by Tomiyama (1998):

(9)

where Eo is the Eötvös number defined as:

(10)

σ being the air/water surface tension coefficient and

(11)

The above equation (9) gives 0 < C_TL < 0.288 for small bubbles that migrate towards the wall and negative values for large distorted bubbles. Following Tomiyama, for an air-water system under atmospheric pressure and room temperature, C_TL changes its sign for a bubble size of D_B = 5.6 mm, which is well enough above the maximum bubble diameter considered in this work, i.e. D_B = 2.6 mm. It is necessary to point out that the values of C_TL are below the factor of 0.5 used in some other works (e.g., Deen et al., 2001; van den Hengel, et al., 2005).

Also, Legendre and Magnaudet (1997), carrying out a theoretical analysis of forces acting on a spherical bubble in a low Reynolds number shear flow, found values of C_TL above those provided by Eq. (9) which tends to 0.5 for bubble Reynolds numbers beyond 20.

In the considered air-water system Eo < 0.14 for bubbles with diameter below 1 mm, which means that we are in the surface tension dominant regime and, therefore, the shape of the bubbles is approximately spherical. In consequence, for the smaller bubbles (i.e., below 1 mm) the transverse lift coefficient is taken as 0.5, in accordance with the previous cited works, but for bubbles of larger diameter (presumably nonspherical) the Tomiyama correlation (9) will be used.

Moreover, C_VM is the virtual mass coefficient, which is taken equal to 0.5 through this paper.

It is necessary to point out that although the simulation strategy considers point-like bubbles, the effect of shape and size of the bubbles is indireC_TLy taken into account in the drag and transversal lift coefficients. Moreover, in the experiments, special care was taken to prevent bubble coalescence (Bröder and Sommerfeld, 2002) which implies that bubbles hardly modify their size inside the (small) bubble column. On the other hand, as it has been said, the low values of Eo (below 0.14) imply that bubble shape is not very far of the spherical one (i.e., spheroids) except, perhaps, the higher diameter class bubbles (2.6 mm) considered in this study, whose number density is very low.

The bubble equation of motion (7) is analytically integrated by assuming that the forces such as gravity and buoyancy, virtual mass, transverse lift and fluid stresses term are constant during the time step. The numerical solution requires that the time step of integration (i.e. the Lagrangian time step Δt_L) is sufficiently smaller than all relevant time scales for the bubble motion which are the integral time scale of turbulence, the bubble response time and the time required for a bubble to cross a grid cell. In the considered vertical bubble-driven flow, the limiting time scale is mostly the bubble response time scale:

(12)

In order to avoid numerical instabilities, the time step was limited to be 25 % of the minimum of the above time scale (Laín and Göz, 2001). For improving numerical efficiency the Lagrangian time step was not fixed, but allowed to vary along the bubble trajectory.

III. EFFECT OF BUBBLES ON LIQUID FLOW

Since the liquid flow in the bubble column is driven by the bubble rise, the source terms due to the bubble phase are essential. As both phases are computed time-dependent, the evaluation of the source terms and the coupling between the phases requires some special treatment in order to yield reasonable averages of the source terms for each control volume, in which bubbles are present. This implies also that all bubbles beeing present in the flow field are tracked simultaneously.

Since the time scale for resolving the fluid flow is generally much larger than that required for bubble tracking a quasi unsteady approach is adopted. This implies a sequential solution of the Eulerian and Lagrangian part with different time steps. Hence, during the Eulerian time step the bubbles are viewing a frozen flow field. The selected Eulerian time step (Δt_E) determines the temporal resolution of the flow fluctuations and its upper limit is given by the constraint that the maximum Courant-Friedrichs-Levy (CFL) number must be less than one. The Lagrangian time step (Δt_L) for calculating the bubble trajectories is generally much smaller than Δt_E and is automatically adjusted along the trajectory according to the criteria mentioned above. Typically about 50 tracking time steps are located within one Eulerian time step.

The calculation of the interaction terms is realised by means of the Particle-source-in-cell (PSI-cell) approximation of Crowe et al. (1977). In the LES scheme adopted here, this model considers the dispersed phase as a local source of momentum. In this context, the expression for the momentum equation source term due to the bubbles is obtained by timeand ensemble-averaging in the following form (Gouesbet and Berlemont, 1999):

(13)

where the sum over n indicates averaging of the instantaneous momentum contributions along the bubble trajectory (i.e. time averaging) and the sum over k is related to the number of computational bubbles passing through the considered cell of size V_cv. The mass of an individual bubble is given by m_k, while N_k is the number of real bubbles contained in one computational bubble. In (13) only the interfacial forces have to be taken into account, so the external forces need to be substracted. Moreover, it should be mentioned that in the two-way coupling procedure no under-relaxation should be applied.

Following the PSI-Cell strategy, the coupling terms are introduced only within the cell where the centre of gravity of the bubbles is located. Let us remark that in Eq. (13) the temporal change of the instantaneous bubble velocity is taken instead of the forces acting on such bubble, because from a Lagrangian perspective this is easier and automatically all forces are accounted for.

IV. NUMERICAL SIMULATIONS

A. Background

At this point some results and conclusions obtained by other authors in similar configurations of bubbly flows employing the LES strategy for the liquid phase are summarised.

Regarding the grid size, Deen et al. (2001) and van den Hengel et al. (2005) conclude that, for vertical bubble-driving flows, relative coarse grids can be used without losing any fundamental characteristics of the flow. This is true in both, Euler-Euler and Euler-Lagrange approaches. For instance, Deen et al. (2001) report calculations in a square cross-sectioned bubble column of 0.15 × 0.15 × 1 m³ musing meshes of 15 × 15 × 100 and 32 × 32 × 100 cells without any significative differences between them. The studies of Milelli et al. (2001) (Euler-Euler), on a bubble plume developing in a cylindrical bubble column, and Hu (2005) (Euler-Lagrange), on the locally aerated bubble column of Becker et al. (1994), present similar conclusions regarding the grid dependence. In addition Milelli et al. (2001) show that for the case of a shear layer laden with bubbles it was also possible to provide an optimum filter width 1.2 < Δ⁄D_B < 1.6. However, this result was not fully supported for the bubble plume, where a coarser grid improved the results in the vertical direction. The constraint imposed on the ratio Δ⁄D_B implied that the interaction of bubbles with the smallest resolved scales is captured without additional approximation. In practice, this condition conflicts with the fact that the computational grid should be large enough when compared to the bubble diameter in order to provide significant statistical samples (Tran, 1997).

The approach adopted here is to establish a grid size larger than bubble diameter but such that the ratio Δ⁄D_B is not very far away from 1.6. Consequently, two grids have been considered in this work: a coarse one with 30 × 30 × 50 cells and a refined one with 45 × 45 × 150 cells which gives for the maximum bubble diameter of 2.6 mm a ratio Δ⁄D_B 1.8 (bubble column diameter 140 mm and height 650 mm).

Regarding the bubble tracking, in the general case, the instantaneous fluid velocity at the bubble location occurring in Eq. (7) is determined from the local resolved fluid velocity linearly interpolated from the neighbouring grid points and a sub-grid root mean square velocity contribution . Van den Hengel et al. (2005) found that for the case of a monodispersed bubble plume developing in the previously mentioned square cross-section bubble column the effect of considering in the bubble equation of motion was negligible in the liquid variables. The same approach is adopted in Hu (2005). This fact was due to the low effective Reynolds number of the liquid flow developing in this kind of bubble-driven flow systems. Therefore, nearly all the liquid energy is contained in the resolved scales which are provided by LES. Unfortunately, no information is given in those works about the influence of the liquid subgrid rms velocity on the bubble variables. Similar behaviour was encountered by Deen et al. (2001) and Milelli et al. (2001) using an Euler-Euler or two-fluid approach.

Therefore, in order to get some insight of the effect of on the bubble variables, in this paper the following simple model for the consideration of the subgrid scale fluctuating liquid velocity effects has been tested:

(14)

where ξ_j is an independent random number sampled from a Gaussian distribution with zero mean and unit variance. On the other hand, the turbulent kinetic subgrid scale k^SGS is taken following the estimation of Lilly (1967):

(15)

A further conclusion of the previously mentioned research is that the values adopted for C_TL influence significatively the performance of the LES simulation. If the transverse lift was neglected, the bubbles did not experience any transversal spreading. This fact could be expected because in absence of transport by fluid fluctuations (or when this effect is not significant), the only way of achieving bubble dispersion is due to the transverse lift force. Nevertherless, the values for C_TL are different: 0.5 in Deen et al. (2001) and van den Hengel et al. (2005), and 0.25 in Milelli et al. (2001). Unfortunately, the work of Hu (2005) does not mention explicitly the values employed for the coefficient C_TL. In the section of results the dependence of the bubble dispersion on the value of C_TL will be illustrated.

B. Flow configuration

The flow configuration consists of a cylindrical laboratory bubble column with a diameter of 140 mm and a height of 650 mm (i.e., water level in the bubble column) corresponding to the experimental set up described in Bröder and Sommerfeld (2002). In such experiment aeration was performed by means of a porous membrane with a diameter of 100 mm and pore sizes of 0.7 µm pretending to establish a homogeneous aeration over the cross-section of the aerator. The gas flow rate could be varied through the supply pressure. For analysing the bubble swarm behaviour and simultaneously evaluating the flow structure and bubbleinduced turbulence in a bubble column Particle Image Velocimetry (PIV) was applied. Full details about the experimental facility and measurement technique can be found in Bröder and Sommerfeld (2002).

C. Simulation set-up

As it has been previously stated, the cylindrical bubble column of diameter 140 mm and height 650 mm was discretised by employing two grids: a coarse one with 30 × 30 × 50 cells and a refined one with 45 × 45 × 150 cells, the latter to check the influence of grid size on the results. The boundary conditions employed for the liquid phase computations were:

• no-slip boundary conditions at the walls and the bottom of the column

• the free surface of the bubble column was also specified as a wall boundary condition with zero stress, which implies a free slip condition.

The bubbles were injected just above the bottom of the bubble column (3 mm above the aerator) over a cross-section with a diameter of 100 mm according to the experiments. The gas phase mass flux was constant across the aerator. The size of the bubbles was sampled stochastically from the measured size distribution. The initial bubble velocity was sampled from a Gaussian distribution with a mean-and rms-value corresponding to the measurements. At the free surface the bubbles are leaving the computational domain. Bubble-wall interactions are modeled as perfectly C_TL elastic collisions.

The calculation procedure is briefly summarised in the following. Firstly, the bubbles are randomly injected at the inlet area and tracked in quiescent liquid for a duration corresponding to the Eulerian time step in order to evaluate the source terms described above and the bubble phase properties. During this first step of the start-up period of course the bubbles do not yet reach the surface of the column. The source terms are used to calculate the fluid flow field until a converged solution is achieved. Then the Lagrangian tracking with injecting new bubbles in each Eulerian time step is continued starting with the bubble locations from the previous Lagrangian computation. With the new source terms the flow field at the next time level is recalculated and so forth. It should be mentioned that with this unsteady procedure no under-relaxation of the source terms should be used.

The experimental case considered is the case 1 presented in Laín et al. (2001). The superficial gas velocity was 0.272 cm/s, the volume flow rate 151 l/h, the average air volume fraction in the column 1.46 % and the range of measured bubble diameters was [0.2, 2.6] mm with a mean bubble diameter of 0.92 mm.

This case has been chosen because the extension of the k - ε model presented in Laín et al. (2002) was not able to provide simultaneously good agreement for the liquid mean vertical velocity and turbulent kinetic energy (Laín et al., 2001) despite the fact that the model performed good enough for operating conditions with smaller bubble sizes. The main idea of using LES for this kind of bubbly flows results from the expectation that the large scale motions (which carry most of the energy) would be primarily responsible for the macroscopic influence of the liquid velocity fluctuations on the bubble motion, including dispersion, whereas small-scale fluctuations would be less important, affecting more the localized bubble oscillations. Therefore, there is a hope that the statistics of velocity fluctuations induced in the liquid by the bubbles motion could be reasonably reproduced.

The choice of the Eulerian time step, Δt_E, is determined by the criterium that the maximum CourantFriedrichs-Levy (CFL) number must be less that one. Because of the low velocities induced in the liquid by the rising bubbles it has been checked that a value of Δt_E = 50 ms is sufficiently small to satisfy the requirement CFL < 1. However, during this investigation it was observed that the magnitude of the liquid fluctuating velocity depends on the value of Δt_E . Asa consequence, reducing the Eulerian time step results in an increase of the fluctuating liquid velocity. Fortunately, the changes of this variable when Δt_E was reduced below 5 ms were very small, which was explicitely checked using Δt_E = 2.5 ms.

The evolution of the quasi-steady flow in the bubble column begins when the first bubbles leave the column. Normally this situation is reached after about 10 s. After that, the data are time averaged up to the end of the simulation. At this stage typically 60,000 to 80,000 computational bubbles are included in the entire liquid domain. The flow is simulated usually for 250 s, but some simulations have been run up to 500 s to be sure that the statistical steady state had been reached. In this case, around one month of CPU time was necessary in a PC Pentium 4 at 1.8 GHz with one processor. Time average (denoted by an overbar) of the liquid velocity field components is performed in each time step as:

(16)

where the averaging is started at time step n0, corresponding in this work to t = 10 s. After that, the velocities at each vertical cross-section are averaged in the angular coordinate, providing a velocity profile versus the radial coordinate.

After performing the sensibility analysis, it was determined that when an Eulerian time step of Δt_E = 5 ms is considered the results for the vertical fluid mean velocity are independent of the grid size. On the other hand, the changes in the velocity profile when the Eulerian time step is halved (i.e., Δt_E = 2.5 ms), using the coarse grid, are quite small. In all the cases the requirement of CFL < 1 is satisfied. Therefore, the following simulations have been performed on the coarse grid with an Eulerian time step of Δt_E = 5 ms. In all these cases the value for the Smagorinsky constant C_s = 0.1 has been adopted.

V. RESULTS

Figure 1 illustrates the temporal evolution of the resolved fluid vertical mean velocity at a monitoring point located at r = 0.01 mm and z = 250 mm above the aerator. It can be seen that the resolved instantaneous velocity (continuous line) evolves in a sort of non-deterministic fashion; however, its cumulative average (dashed line) attains a quite stable value after 100 s of averaging. The typical temporal evolution of the flow structure in the bubble column is shown in Fig. 2 by plotting contours of the (resolved) instantaneous vertical velocity of the fluid at three different equidistant times. The red color in the slices indicates upwards velocities and blue color downwards velocities. In general, the over-all flow pattern consists of a large recirculating loop where the flow is directed upwards in the core of the column and downwards near of the walls. It can be seen in some areas, which change with time, that the liquid may also rise near the wall, especially in the region immediately above the aerator. Eventually, the average over a long enough time period gives the usual recirculation pattern where liquid rises in the center and goes down near the walls (Fig. 2, right).

Figure 1: Instantaneous resolved fluid velocity (continuous line) and cumulative averaged fluid velocity (dashed line) at monitoring point (z = 250 mm above the aerator, r = 0.01 mm).

Figure 2: Instantaneous liquid vertical velocity con-tours for three different times: 125 s, 150 s and 175 s (from left to right) and averaged vertical fluid velocity (right).

However, the averaged fluid flow pattern in the bubble column is not always completely axisymmetric, especially in the zone just above the aerator. This fact has been already experimentally noticed by Sanyal et al. (1999) using the CARPT (Computer Automated Radioactive Particle Tracking) technique. The asymmetry of averaged fluid velocities found in the experiments just above the aerator is evident in the left part of Fig. 3, where the flow crosses the centre line. The average fluid velocity field is, however, axisymmetric C_TL. Mean vertical fluid velocities (top) and vertical bubble velocities (bottom) (z = 495 mm). Calculation with the experimental bubble size distribution.

beyond the cross-section of 100 mm. Interestingly, the numerical computations are able to capture such behaviour.

Figure 3: Experimental (left) and numerical (right) vector field of mean fluid velocity.

In the following, the time and azimuthally averaged velocity profiles versus radius, for both phases at the z = 495 mm section above the aerator, are considered. Figure 4 shows the results of LES in combination with several transverse lift coefficients C_TL in order to illustrate the role of the transverse lift force in the bubble dispersion through the column cross section. Three cases have been considered: C_TL = 0, the Tomiyama expression (9) and finally the values obtained multiplying (9) by a factor of five in order to get a large value for C_TL. The mean vertical liquid velocities are in the upper graph and the bubble velocities are shown in the lower plot. The bubble mean velocities show very clearly the effect of the transverse lift coefficient value. If no transverse lift force is considered, the bubbles rise up to the top of the bubble column without significant spreading, something noticed previously by Deen et al. (2001), for instance. Therefore, the mean velocity keeps a constant value up to a radius of 0.055 m, which is slightly larger than the aerator size, and then decreases fast up to zero. The liquid velocities induced by the bubble rise in this case presents a nearly flat profile in the core region and then the velocity decreases gradually to negative values in the zone of downflow and reaches eventually the zero value on the wall. The fluctuations of the liquid velocity observed for some of the profiles might be caused by small scale oscillations of the flow.

Figure 4: Influence of the transverse lift coefficient C_TL. Mean vertical fluid velocities (top) and vertical bubble velocities (bottom) (z = 495 mm). Calculation with the experimental bubble size distribution.

When a value C_TL 0 is used, the bubbles experience a transversal spreading towards the wall, as the bottom part of Fig. 4 shows. Obviously, increasing the value of C_TL implies stronger lateral dispersion of the bubbles. It should be noted that the lift force in the transverse direction depends on the value of the slip velocity in the vertical component, which is usually much larger than the horizontal components, and that, on the other hand, the Tomiyama correlation Eq.(9) provides larger values of C_TL for the bigger bubbles. These two effects combine to produce a more efficient spreading of the larger diameter bubbles than of the smaller diameter ones (i.e., the bubble mean diameter increases towards the wall), resulting in the occurrence of a bubble velocity peak near the wall. Computations with a constant lift coefficient for all bubble sizes, which are not shown here, revealed that the peak near the wall still exists. Such a peak is due to the fact that rising velocity of the bigger bubbles is larger than that of the smaller ones. As C_TL in-creases, the peak displaces towards the wall and the center values of mean and rms bubble velocities de-crease. Moreover, this fact implies a smaller momentum transfer from the bubbles to the liquid and also the mean values of liquid velocities are decreased (Fig. 4 top).

This behaviour reveals that the lateral spreading of the bubbles is governed exclusively by the transverse lift force, when the magnitude of subgrid rms velocities is negligible compared to the resolved velocities. If, for instance, a k - ε turbulence model is employed, the dispersion process determined by the magnitude of the turbulent kinetic energy provides a mechanism for spreading the bubbles which turns to be more efficient than the transverse lift force.

However, the experiments show a monotonically decreasing of the bubble mean velocity profile from the symmetry axis towards the wall. Therefore, the bubble velocity peak in the vicinity of the wall is not realistic and should be solved somehow. Here, two strategies have been implemented in order to improve the spreading of the bubbles across the bubble column cross-section: the first one is to include the effect of the sub-grid root mean square velocity in the bubble equation of motion through the simplified expression (14). The second strategy is to implement the proposed correlation for the transverse lift coefficient C_TL (assuming a value of C_TL = 0.5 for spherical bubbles with diameter below 1 mm) and the Tomiyama correlation (9) for larger bubbles.

The effect of considering and the C_TL formulation on bubble dispersion is demonstrated in figure 5. The upper graph shows the mean bubble diameter distribution (azimuthally and time averaged) across the section z = 495 mm above the aerator, whereas the lower plot shows the average gas void fraction. It is readily seen that the effect of on the mean bubble diameter and gas void fraction profiles is marginal (compare continuous and dotted curves). The bubble diameter near the wall is remarkably higher than in the core region and the decreasing of gas volume fraction towards the wall is small. On the other hand, the performance of the proposed C_TL formulation is much more promising (compare dashed and dotted curves). Due to the enhancement of the transverse lift force for the spherical smaller bubbles, the bubble diameter profile is nearly flat and the slope of the gas void fraction profile is increased implying a higher gas mass flux in the core of the column when compared with that provided by the original Tomiyama correlation (9). Consequently, the proposed C_TL correlation will be used hereafter in the simulations.

Figure 5: Influence of the transverse lift coefficient formulation and effect of the subgrid root mean squarevelocity uon the distribution of bubble diameter (up) and gas void fraction (down) over the crosssection z = 495 mm above the aerator, calculations with the experimetnal size distribution.

In the following, the time and azimuthally averaged velocity profiles versus radius, for both phases at the cross-section z = 495 mm above the aerator, are compared with the experimental data and k - ε results.Figure 6 shows the mean velocity profiles for liquid (top) and bubbles (bottom). The liquid velocities induced by the bubble rising in this case presents a nearly flat profile near of the centre-line and then velocities decrease gradually to negative values in the zone of downflow and reach eventually the zero value on the wall. As it can be readily seen in that figure, the agreement of the liquid velocity with the measurements is reasonably good at z = 495 mm. The discrepancies in the zone of downflow are due to the fact that the fluid velocity profile does not satisfy overall mass conservation because the measurements were only possible for a single vertical plane through the center of the column. Therefore, an experimental averaging of the liquid velocity profile in the circumferential direction could not be realised.

Figure 6: Comparison of results for the k - ε model and LES against experimental data at section z = 495 mm above the aerator. Fluid variables (top) and bubble velocities (bottom). Calculations with bubble size distribution, lift coefficient according to proposed correlation and fluid bubble drag.

On the other hand, the bubbles mean velocity agrees well with the measurements for the section at z = 495 mm above the aerator. It is important to emphasise that no bubble velocity peak is observed near the wall, mainly due to the fact that the mean bubble diameter is roughly constant over the entire cross-section.

Moreover, Fig. 6 shows the liquid fluctuating kinetic energy (top, red curves) and the bubble fluctuating vertical velocity (bottom). In the scheme followed in this work, the liquid velocity fluctuations are calculated during the simulation as:

(17)

Considering the upper plot of that figure, it can be seen that the simulated mean liquid fluctuating kinetic energy, , captures fairly well the tendence of the experimental PIV measurements. Also, the bubble vertical fluctuating velocity compares reasonably well with the experimental data, showing a flat profile as the experiments.

Finally, also in Figure 6 can be seen the quantitative comparison between the LES, experimental data and the k - ε two-dimensional simulations reported in Laín et al. (2001) at section z = 495 mm above the aerator. The upper part of Fig. 6 shows that the mean vertical liquid velocities (left vertical axis) provided by the 3D LES calculations are in good agreement with both, experimental data and results of the k - ε turbulence model. Also, the fluctuating kinetic energy of the liquid is shown in the same plot (right vertical axis). Here, it is observed that the LES simulations provide a very good estimation of that variable, improving drastically the result of the 2D simulations performed with the k - ε model. On the other hand, the bubble vertical mean velocity and its fluctuating component are quite similar to the values provided by the k - ε model and the experiments, although LES gives a better prediction of the mean velocity.

VI. CONCLUSIONS

In this work results on the simulation of the threedimensional transient flow developing in a cylindrical laboratory bubble column have been presented. The calculation scheme is based on the combination of Large Eddy Simulation for the liquid phase and a Lagrangian approach for the dispersed gas phase including coupling through the momentum source terms. Previous conclusions available in the literature have been also reproduced here. Interestingly, the consideration of a simple model to account for the so-called bubble induced turbulence, linked to the subgrid liquid fluctuating velocity, did not modify significantly the performance of the simulation, regarding both fluid and gas variables, confirming that this effect does not influence bubble dispersion in this kind of flow configuration.

As in previous works, a strong dependency of the bubble dispersion in the column on the value of transverse lift force coefficient, C_TL has been found. It is concluded that the transverse lift, which depends on the bubble-liquid relative velocity, is the main mechanism responsible for the spreading of the bubbles across the column cross-section.

As a result, the three-dimensional LES reproduces fairly well the dynamic structure of the liquid flow pattern. Moreover, the long term steady state recirculating flow configuration provided also good quantitative comparison with the experiments.

The main drawback of the considered approach is the long CPU times needed. However, this fact can be partially paliated by parallelizing the bubble tracking module, a task that is currently being carried out.

REFERENCES

1. Becker, S., A. Sokolichin and G. Eigenberger, "Gas-liquid flow in bubble columns and loop reactors: Part II. Comparison of detailed experiments and flow simulations," Chemical Engineering Science, 49, 5747-5762 (1994).         [ Links ]

2. Bourloutski, E. and M. Sommerfeld, "Parameter studies on the three-dimensional calculation of bubble colums." Joint US ASME/European Fluids Engineering Summer Conference, Montreal, Paper No. FEDSM 2002-31218 (2002).        [ Links ]

3. Bröder, D. and M. Sommerfeld, "An advanced LIFPLV system for analysing the hydrodynamics in a laboratory bubble column at higher void fraction," Experiments in Fluids, 33, 826-837 (2002).         [ Links ]

4. Crowe, C.T., M.P. Sharma and D.E. Stock, "The Particle-Source-in-Cell (PSI-Cell) method for gasdroplet flows," J. Fluids Eng., 99, 325-332 (1977).        [ Links ]

5. Crowe, C.T., M. Sommerfeld and Y. Tsuji, Multiphase flows with droplets and particles, CRC Press (1998).        [ Links ]

6. Deen, N.G., T. Solberg and B.H. Hjertager, "Large eddy simulation of the gas-liquid flow in a square cross-sectioned bubble column," Chem. Eng. Sci., 56, 6341-6349 (2001).        [ Links ]

7. Gouesbet, G. and A. Berlemont, "Eulerian and Lagrangian approaches for predicting the behaviour of discrete particles in turbulent flows," Prog. Energy Combust. Sci., 25, 133-159 (1999).        [ Links ]

8. Hu, G., Towards Large Eddy Simulation of dispersed gas-liquid two-phase turbulent flows, Ph.D. Thesis, West Virginia University, Morgantown WV (2005).         [ Links ]

9. Jakobsen, H.A., B.H. Sannaes, S. Grevskott and H.F. Svendsen, "Modeling of vertical bubble-driven flows," Industrial and Eng. Chemistry Res., 36, 4052-4074 (1997).         [ Links ]

10. Laín, S., D. Bröder and M. Sommerfeld, "Numerical simulations of the hydrodynamics in a bubble column: Quantitative Comparisons with experiments," Proc. of the 4th Int. Conf. on Multiphase Flow, May 27 -June 1, New Orleans, USA (2001).         [ Links ]

11. Laín, S. and M.F. Göz, "Numerical instabilities in bubble tracking in two-phase flow simulations," Int. J. Bifurcation and Chaos, 11, 1169-1181 (2001).        [ Links ]

12. Laín, S., D. Bröder, M. Sommerfeld and M.F. Göz, "Modelling hydrodynamics and turbulence in a bubble column using the Euler-Lagrange procedure," Int. J. Multiphase Flow, 28, 1381-1407 (2002).        [ Links ]

13. Lilly, D.K., "The representation of small scale turbulence in numerical simulation experiments," Proc. IBM Scientific Computing Symp. on Environmental Sciences, 320 195-210 (1967).         [ Links ]

14. Legendre, D. and J. Magnaudet, "A note on the lift force on a spherical bubble or drop in a low Reynolds number shear flow," Phys. Fluids, 9, 3572-3574 (1997).        [ Links ]

15. Milelli, M., B.L. Smith and D. Lakehal, "Some new approaches to bubble plume modeling using CFD," Proc. Int. Mechanical Eng. Congress and Exposition, New York (USA), November 2001.        [ Links ]

16. Sanyal, J., S. Vásquez, S. Roy and M.P. Dudukovic, "Numerical simulations of gas-liquid dynamics in cylindrical bubble column reactors," Chem. Eng. Sci., 54, 5071-5083 (1999).        [ Links ]

17. Smagorinsky, J., "General circulation experiments with the primitive equations," Monthly Weather Review, 91, 99-165 (1963).        [ Links ]

18. Tomiyama, A., "Struggle with computational bubble dynamics," Proc. of the 3th Int. Conf. on Multiphase Flow, June 8-12, Lyon, France (1998).        [ Links ]

19. Tran, M.L., Modélisation instationnaire de la distribution spatiale des phases dans les écoulements diphasiques en régimes à bulles, (in French) Ph.D. Thesis, University Claude Bernard, Lyon (1997).        [ Links ]

20. van den Hengel, E.I.V., N.G. Deen and J.A.M. Kuipers, "Application of coalescence and breakup models in a discrete bubble model for bubble columns," Ind. Eng. Chem. Res., 44, 5233-5245 (2005).         [ Links ]

Received: September 25, 2008
Accepted: October 31, 2008
Recommended by Subject Editor: Eduardo Dvorkin