Chemical Engineering Journal 139 (2008) 589–614 CFD simulation of bubble column—An analysis of interphase forces and turbulence models Mandar V. Tabib, Swarnendu A. Roy, Jyeshtharaj B. Joshi ∗ Institute of Chemical Technology, University of Mumbai, Matunga, Mumbai 400019, India Received 25 June 2007; received in revised form 8 September 2007; accepted 10 September 2007 Abstract 3D transient CFD simulations of bubble column have been performed for a wide range of superficial gas velocity on an industrially relevant cylindrical column and the CFD predictions have been compared with the experiments of Menzel et al. [T. Menzel, T. Weide, O. Staudacher, U. Onken, Reynolds stress model for bubble column reactor, Ind. Eng. Chem. Res. 29 (1990) 988–994]. Simulations have also been performed to understand the sensitivity of different interphase forces (drag, lift, turbulent dispersion and added mass). This work highlights the importance of choosing the CL value and the drag law in accordance with the bubble size. Further, a laboratory scale bubble column with three different spargers (perforated plate, sintered plate and single hole) has been simulated using three different turbulence closure (k–ε, RSM and LES) models, with the purpose of critically comparing their predictions with experimental data [M.R. Bhole, S. Roy, J.B. Joshi, Laser doppler anemometer measurements in bubble column: effect of sparger, Ind. Eng. Chem. Res. 45 (26) (2006) 9201–9207; A.A. Kulkarni, K. Ekambara, J.B. Joshi, On the development of flow pattern in a bubble column reactor: experiments and CFD, Chem. Eng. Sci. 62 (2007) 1049–1061]. It has been found that the RSM shows better agreement than the k–ε model in predicting the turbulent kinetic energy profiles. Comparatively, the LES has been successful in capturing the averaged behavior of the flow, while at some locations; it slightly over predicts the kinetic energy profiles. Further, it has been able to simulate the instantaneous vortical-spiral flow regime in case of sieve plate column, as well as, the bubble plume dynamics in case of single hole sparger. Thus, LES can be effectively used for study of flow structures and instantaneous flow profiles. © 2007 Elsevier B.V. All rights reserved. Keywords: Bubble column; CFD; Sparger; Interphase force; Turbulence models; k–ε; LES 1. Introduction The processes involving reactions between liquid and gas phases are ubiquitous feature of the chemical process industry. Many types of reactors are used for carrying out such reactions, and amongst them, the bubble column reactors find wide spread applicability, a summary of which can be seen in Shah et al. [4], Doraiswamy and Sharma [5], and Deckwer [6]. In a bubble column reactor, the density difference between the gas and liquid induces the required stirring action, thus offering an attractive way to carry out gas–liquid and gas–liquid–solid reactions. Moreover, the use of bubble column reactor is advantageous because of its simple construction, absence of any moving parts and small floor space requirements. However, because of the simple construction, bubble column reactors also have an ∗ Corresponding author. Tel.: +91 22 2414 5616; fax: +91 22 2414 5614. E-mail address: jbj@udct.org (J.B. Joshi). 1385-8947/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cej.2007.09.015 inherent limitation of having fewer degrees of freedom available to tailor the performance characteristics. In this type of reactor, local flow, turbulence and gas hold-up distribution are interrelated in a complex way with the operating and design variables; hence the extensive knowledge of prevailing hydrodynamics is crucial. Development of detailed fluid dynamic model is therefore essential to understand these complex interactions, which is beneficial for the reliable and efficient design of these reactors. The complex hydrodynamics of bubble columns has inhibited the development of design procedures from first principle. Hence, during the past 45 years, vigorous attempts have been made to understand the flow pattern using various techniques of flow visualization and mathematical modeling in bubble columns. Joshi [7] and Sokolichin and Eigenberger [8] have critically analyzed the investigations and have made a coherent presentation of the current status. Joshi [7] has classified the entire modeling effort into three phases. Phase I and II models were made simple by making many assumptions which ultimately incorporate fair degree of empiricism. In the Phase III 590 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Nomenclature C0 , C1 CB CD CL Cs CTD CVM Cε1 C␧2 Cµ Cµ,BI dB do D Eo EI g G H HD k MI MD ML MTD MVM P P′ r R ReB S t u u′ uinst uL uSGS UV vL V VS VT W z drift-flux constants (dimensionless) interface energy transfer factor (dimensionless) drag force coefficient (dimensionless) lift force coefficient (dimensionless) Smagorinsky constant (dimensionless) turbulent dispersion coefficient (dimensionless) added mass force coefficient (dimensionless) model parameter in turbulent dissipation energy equation (dimensionless) model parameter in turbulent dissipation energy equation (dimensionless) constant in k–ε model (dimensionless) constant in bubble induced turbulence model (dimensionless) bubble diameter (m) sparger hole diameter (m) diameter of the  column (m)  Eotvos number = g(ρL −ρG )dB2 σ (dimensionless) rate of energy supply from the gas phase to the liquid phase (kg m2 /s3 ) gravitational constant (m/s2 ) generation term (kg/m s2 ) height (m) height of dispersion (m) turbulent kinetic energy per unit mass (m2 /s2 ) total interfacial force acting between two phases (N/m3 ) drag force (N/m3 ) lift force (N/m3 ) turbulent dispersion force (N/m3 ) added mass force acting (N/m3 ) pressure (N/m2 ) exact production term (kg/ms3 ) radial distance (m) column radius (m) Reynolds number (= dB VS /ν) (dimensionless) strain rate (1/s) time (s) velocity vector (m/s) fluctuating velocity (m/s) instantaneous velocity (m/s) average axial liquid velocity (m/s) sub-grid scale velocity (m/s) axial-radial Reynolds stress centerline axial liquid velocity (m/s) superficial velocity (m/s) axial slip velocity between gas and liquid (m2 /s) terminal rise velocity (m2 /s) width of the rectangular column (m) axial distance along the column (m) Greek symbols ∆ Grid size (m) t x y z r θ z ε φj µ µBI µeff µT θ ρ ρu′i u′j σ σε σk σθ2 τk ǫ ǭ time step used in simulation (s) grid spacing in x direction (m) grid spacing in y direction (m) grid spacing in z direction (m) grid spacing in radial direction (m) grid spacing in tangential direction (m) grid spacing in axial direction (m) turbulent energy dissipation rate per unit mass (m2 /s3 ) pressure-strain correlation (kg/m s3 ) molecular viscosity (Pa s) bubble induced viscosity (Pa s) effective viscosity (Pa s) turbulent viscosity (Pa s) tangential co-ordinate (m) density (kg/m3 ) Reynolds stresses (N/m2 ) surface tension (N/m) Prandtl number for turbulent energy dissipation rate (dimensionless) Prandtl number for turbulent kinetic energy (dimensionless) dimensionless variance, defined in Eq. (14) shear stress of phase k (Pa) fractional phase hold-up (dimensionless) average fractional phase hold-up (dimensionless) Subscripts k phase (k = G: gas phase, k = L: liquid phase) models, the foremost emphasis has been given on the reduction of empiricism by ensuring completeness of the formulation of continuity and momentum equations. Moreover, the comprehensive developments in the computational fluid dynamics (CFD) in the past few decades have strengthened the Phase III models, which has generated hope of complete understanding of fluid mechanics with a caution of the arduous task ahead thus giving a better knowledge of existing hydrodynamics. It is well known that a major problem in realizing a CFD code is to capture the physics behind the phenomena occurring in bubble columns. The interaction between the dispersed gas phase and the continuous liquid phase affects the interphase forces (e.g. drag force, lift force and added mass force) and turbulence in the column. Therefore, the correct modeling of interphase forces and turbulence is of prime importance for capturing the physics correctly. Several models for interphase forces have been reported in the literature, a detailed account of which is provided by Joshi [7] and recently by Rafique et al. [9]. Another critical issue in the CFD modeling of bubble column is the proper description of turbulence in the continuous phase. Researchers have tested different closures [standard k–ε model, Reynolds Stress Model (RSM), Large Eddy Simulation (LES)] for turbulence, and amongst these, the standard k–ε model is the most adopted one owing to its simplicity and lesser computational M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 requirement [7,9]. A brief account of the available literatures has been given in the next section. The present paper makes an attempt to present a sensitivity analysis of the different interphase forces and the turbulence closures within the Eulerian–Eulerian framework of CFD modeling. The CFD study of sensitivity of different interphase forces, like drag, lift and added mass have been performed for the experimental data of Menzel et al. [1], which is an industrially relevant cylindrical column. The understanding was then extended to study the performance of different turbulence closures in simulating the flow development arising out of three different spargers namely a perforated plate (multipoint sparger), a single hole and a porous plate in an laboratory scale bubble column [2,3]. The relative merits of the various force formulations and the turbulence models have been brought out and finally, some recommendations have been made for appropriate selection. 2. Literature review A wide range of CFD studies on bubble column reactors have been conducted in recent years. Mainly two approaches, namely Euler–Euler [7,8,11–13] and Euler–Lagrange [14–19], have been adopted to simulate a gas–liquid system operating under different conditions. The Lagrangian approach gives a direct physical interpretation of the fluid–particle interaction, but it is computationally intensive, and hence cannot be used for simulating large columns, and systems having high dispersed phase volume fraction, as is the case in the present study. Although researchers have carried out both 2D [8,11,12,14,20–26] as well as 3D [10,11,19,24,29–31] Euler–Euler CFD simulations, the results of 3D CFD simulations were found to give good agreements with the experimental results [11,13,20,24,32,33]. Therefore, in this paper we will restrict ourselves in understanding the weakest link in the 3D Euler–Euler CFD simulations by analyzing the published literature, and ultimately make an effort to strengthen it. A brief description of these studies has been given below. Mudde and Simonin [11] have performed transient simulation using standard k–ε turbulence model to simulate the experimental data of Becker et al. [35]. In the interphase force formulation, they had accounted for drag force and virtual mass force, whereas lift force was neglected. The effect of virtual mass force on the simulation was studied, and they concluded that the inclusion of virtual mass force gives better results, which resembles experimental data closely. A low Reynolds number k–ε model was also tested, but it did not lead to any significant improvement in results. Plfeger et al. [24] have carried out transient simulation in a rectangular bubble column with three different types of spargers. Their main focus was on studying the influence of turbulence modelling. Hence, both the laminar and turbulent simulations were performed. A standard k–ε model was used to describe turbulence occurring in the continuous fluid. The results were compared with the experimental measurements carried out using LDA and PIV. The laminar model did not show the harmonic oscillations as observed in experiments. While, the standard k–ε 591 model was found to give good agreement with the experimental liquid velocity and the turbulent kinetic energy profiles. In this work, the interphase momentum transfer was taken into account by considering the drag force, whereas the lift and the added mass forces were neglected. They observed that the inclusion of gas phase dispersion does not have much effect. Deen et al. [34] have simulated the gas–liquid flow in a rectangular bubble column. They have incorporated an expression for bubble induced turbulence along with standard k–ε model. The lift force and virtual mass were neglected, and only drag force was considered. They have studied the bubble plume oscillation, and found that the bubble plume, at lower height of the column, gets positioned permanently in one corner, giving rise to asymmetric velocity profiles. This was attributed to the overestimation of viscosity by k–ε model. The predicted velocity profiles agreed well with the experimental results. Deen et al. [36] have extended the work of Deen et al. [34] and compared the performance of LES and k–ε turbulence models for the bubble column reactor. Unlike their previous study [34], here they have considered both the lift and virtual mass force. They found little influence of the virtual mass force on the results owing to the quasi-stationary state (a k–ε model simulation result that showed low values of acceleration due to unresolved transient details), but the incorporation of lift force was found to improve the results considerably. The LES model was able to capture the transient movement of the bubble plume much better than the k–ε model. The LES results of both velocity and velocity fluctuations show better quantitative agreement with the experiments. Krishna and Van Baten [32] have modeled air–water two phase flows in bubble column as three phase flow by considering two separate bubble class and a liquid phase, but they have not taken the interactions between two bubble classes into account. The turbulence in the liquid phase was described using standard k–ε model, and only the contribution of drag force in the interphase momentum transfer was considered. The lift force was not considered due to the uncertainty in assigning the value. Transient simulations were performed and a good agreement has been shown between the predicted gas hold-up profile and the experimental results of Hills [37]. They have also shown favourable agreement of the average liquid velocity profile and average hold-up with their own experimental results. Pfleger and Becker [27] have performed both experimental and numerical studies in a cylindrical bubble column operating in the homogeneous regime. They have carried out dynamic simulation using standard k–ε model for describing turbulence in the liquid phase. Lift force and the added mass force were neglected, and a constant drag co-efficient of 0.44 was taken to describe the interphase momentum transfer. Their main aim was to study the effect of bubble induced turbulence on the correctness of turbulence results. The simulation results were found to over predict the overall gas hold-up, which they attributed to the grid coarseness and simplified gas inlet modeling. The radial profile of axial velocity and gas hold-up agree well with their experimental results. They concluded that the incorporation of bubble induced turbulence enhances the predictability of radial profile of axial velocity, and further they emphasized the 592 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 need of carrying out dynamic simulations for obtaining good results. Ranade and Tayalia [38] investigated the influence of two different ring spargers on the hydrodynamics of a bubble column using standard k–ε model. Drag formulation proposed by Schwarz and Turner [39] was used and all the other forces were neglected. They observed that the relative mixing performance of single- and double-ring spargers predicted using three-dimensional geometry was much closer to the experimental observations than that predicted by axis-symmetric, two-dimensional geometry. Their computations over predict the overall gas hold-up by 90%. This was attributed to the use of Schwarz’s simplified drag law which does not consider the liquid wake effect on effective slip velocity. Hence, this results in the under-prediction of slip velocity and therefore, an overprediction of gas volume fractions. This work brings out the need to select a proper drag model, and emphasizes on carrying out three-dimensional simulations to get an accurate result. Buwa and Ranade [28] have studied the effect of gas velocity, sparger design on the gas–liquid flow in a rectangular bubble column. The turbulence in the liquid phase was modeled using standard k–ε model, and the bubble induced turbulence was neglected owing to the fact that the extra dissipation due to the small-scale structure compensates the extra generation of turbulence due to large bubble. Both lift and added mass force were considered. They have compared their simulation results of plume oscillation period and radial profile of axial liquid velocity with the experimental results of Pfleger et al. [24] and the gas hold-up profile was compared with the data of Buwa et al. [40]. The asymmetry in experimental and predicted profiles of liquid velocity and gas hold-up was attributed to the insufficient time averaging. Further, they have studied the influence of bubble size on the plume oscillation by considering two limiting value of bubble size (0.5 mm and 20 mm) and found higher oscillation period for larger bubble size. They suggested the use of multi-group CFD model for better representation of different sparger. Lakehal et al. [41] were the first ones to employ LES model for a bubbly flow. Their system involved a vertical bubbly shear flow at a very low void fraction (1.9%). From their study, they suggested that for obtaining better results the optimum cut off filter should be 1.5 times the bubble diameter. They observed that the dynamic approach of Germano does not perform better than Smagorinsky model, which they felt could be due to the inadequate dimensions of their computational domain. This study suggested that the LES approach can be promising for predicting the phase velocities and the void fraction, and it emphasized on the need to carry out more LES simulations on complex systems such as buoyancy driven flow operating at higher void fractions. In this context, Deen et al. [36] were the first one to apply the LES model to simulate a bubble column, and they reported observing a better resolved flow using LES than k–ε model, as discussed previously. Bove et al. [42] used the same geometry as Deen et al. [36] to study the influence of different numerical schemes, of different drag models and of initial flow conditions on LES performance. They observed that the use of a second order FCT scheme for LES simulation enhanced its performance. But they suggested the need to carry out more tests on LES with higher order schemes and finer grid resolutions, before identifying the best numerical scheme for LES simulations. It was seen that the LES results were very sensitive to initial boundary conditions. In this work it seems that the sparger (individual holes in it) was not modeled due to difficulties in adapting the mesh grid to the geometry. The modeling of holes as source points could have made this work more interesting. This work also suggested that there is a need to pay attention towards near wall region as the sub-grid scale models do not account for near wall region processes which can lead to erroneous prediction of frictional stresses at the wall. Bombardelli et al. [43] used a finite element based LES code for simulating and analyzing the phenomena of wandering bubble plumes as in experiments conducted by Becker et al. [35]. They observed that LES was been able to capture the instability associated with wandering more vividly as compared to k–ε simulation. They analyzed the plume for coherent structures arising out of interplay of eddies and bubbles by studying the vorticity contours. Their work showed the evolution of number of eddies and their interrelation with the bubble plume. In this work, the velocity conditions at the wall were set so as to enforce the law of the wall using the LES-Near Wall Modeling (NWM) approach, but the effect this approach has on the simulation results, as compared to the conventional LES-Near Wall Resolution (NWR) approach is not yet known. Zhang et al. [44] took forward the work of Deen et al. [36] by investigating the effect of Smagorinsky constant and the interfacial closures for drag, lift and virtual mass force in two columns having different aspect ratio (H/D = 3.6). They have considered two set of interfacial closures: (1) as proposed by Tomiyama [45], and (2) the Ishii and Zuber [46] drag model, along with lift and added mass forces having a constant coefficient value of 0.5. They observed that, by increasing the value of Smagorinsky constant, the bubble plume dynamics dampens, and consequently a steep mean velocity profile is seen. They obtained good agreement with experimental results when the value of Smagorinsky constant was used in the range of 0.08–0.10. They also suggested that, for taller columns (H/D = 6), interfacial force closures proposed by Tomiyama [45] gave better agreements with experiments. As a future work, they have suggested the need for simulating the dynamic free liquid surface at the top and observing its effect on the velocity profiles in the top region. Kulkarni et al. [3] have presented a comprehensive study of flow profile development both experimentally and numerically for two different spargers. They have solved steady state momentum balance equation with standard k–ε model. Both the lift force and added mass force were considered, and the dispersion of gas was taken into account through a dispersion coefficient. They have captured the development of flow profile, and compared the experimental data of radial profile of axial liquid velocity, gas hold-up turbulent kinetic energy, eddy viscosity and eddy dissipation for single point and a multipoint sparger with the simulation at different H/D. The simulated results compares well with the experimental value at higher H/D. They concluded that the prediction of hydrodynamics near sparger requires further M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 investigation and also pointed out the need of sensitivity analysis of different interfacial forces. Dhotre et al. [47] have assessed the performance of two subgrid scale LES models. They observed that the Smagorinsky model with model constant of Cs = 0.12 performed quite well and gave results identical to those given by the dynamic procedure of Germano. The conceptual advantage that lies in using Germano model is that it estimates Cs value which is not known a priori. Further, it was observed that in comparison with the modified k–ε model (involving extra source terms for turbulent kinetic energy and dissipation rate), the standard k–ε model gives reasonably good agreement with the mean experimental data except for the radial and axial distribution of the fluctuating liquid velocity and turbulent kinetic energy close to the wall. A summary of the published work has been given in Table 1. From Table 1 and the foregoing discussion, it is clear that further investigations are required for understanding the following points: (1) Researchers have used a wide range of formulations for the drag coefficient in order to satisfy the experimental condition under consideration. In many of these studies, there is no valid justification given for the selection of the drag law and the bubble diameter. The drag law is used to simulate the slip velocity, which is dependent upon the bubble size in the system. But, many times, the mean bubble diameter is not selected as per the physical reality, and still one may be able to simulate the slip velocity by choosing a different drag law. Hence, the issue of proper selection of bubble diameter and drag law needs to be brought out. (2) Use of lift force is still not well understood. Several researchers have neglected the lift force due to the lack of understanding and several other researchers have used a constant value of lift coefficient, but no systematic detailed effort has been made to find out how CFD results behave with different value of lift coefficient having different signs. (3) Similarly, for the case of turbulent dispersion force also, no sensitivity analysis has been performed. (4) The performance of different turbulence models like k–ε and LES has been tested and compared. However, the relative merits are far from total clarity. Further, the quantitative performance of the RSM model needs to be included in the overall analysis. (5) Further, despite several notable works involving LES on bubble column, we still do not have enough indication on predictive performance of LES at higher superficial gas velocity and higher void fractions for cylindrical columns. Most of these simulations work have been carried out on rectangular bubble column, and industrially relevant cylindrical columns have not been considered. Also, most of the work is concerned with time averaged profiles, and one needs to look at the use of instantaneous profiles for study of flow structures and flow patterns using LES. In spite of several publications describing the 3D CFD simulation of bubble column in an Euler–Euler framework, a systematic analysis of the different interphase forces and tur- 593 bulence closure in an industrially relevant column is scarce. Therefore, in this paper, an attempt has been made to study the sensitivity of different interphase forces (drag, lift, virtual mass and turbulent dispersion) on simulation results from a column of typical industrial size [1] and operating under heterogeneous regime. Further, the comparative performance of three different turbulence closure models (k–ε, RSM and LES) in capturing the effect of three different spargers (perforated plate, single hole and sintered plate) on hydrodynamics of bubble column has been presented. 3. Mathematical model The numerical simulations presented are based on the twofluid model with the Euler–Euler approach. Here, each fluid (or phase) is treated as a continuum in any size of the domain under consideration. The phases share this domain and interpenetrate as they move within it. The Eulerian modelling framework is based on ensemble-averaged mass and momentum transport equations for each of these phases. These transport equations without mass transfer can be written as: Continuity equation ∂ (ρk ǫk ) + ∇(ρk ǫk uk ) = 0 ∂t (1) Momentum transfer equations ∂ (ρk ǫk uk ) + ∇(ρk ǫk uk uk ) ∂t = −∇(ǫk τk ) − ǫk ∇P + ǫk ρk g + MI,k (2) The terms on the right hand side of Eq. (2) are, respectively, representing the stress, the pressure gradient, gravity and the ensemble-averaged momentum exchange between the phases, due to interface forces. The pressure is shared by both the phases. The stress term of phase k is described as follows:   2 T τk = −µeff,k ∇uk + (∇uk ) − I(∇uk ) (3) 3 where µeff,k is the effective viscosity. The effective viscosity of the liquid phase is composed of three contributions: the molecular viscosity, the turbulent viscosity and an extra term due to bubble induced turbulence. µeff,L = µL + µT,L + µBI,L (4) The calculation of the effective gas viscosity is based on the effective liquid viscosity [21] as follows: µeff,G = ρG µeff,L ρL (5) There are several models to take account of the turbulence induced by the movement of the bubbles. In this study, the model proposed by Sato and Sekoguchi [48] was used. µBI,L = ρL Cµ,BI ǫG dB |uG − uL | with a model constant C␮,BI = 0.6. (6) 594 Table 1 Summary of previous work Geometry Operating details Turbulence model Drag Lift VMF Grid size (∆t) Mudde and Simonin [11] Rectangular column (Becker et al. [35]): W = 0.5 m; D = 0.08 m; H = 1.5 m. Sparger: do = 0.04 m (100 mm from the left wall) Rectangular column: W = 0.2 m; D = 0.05 m; H = 0.45 m. Sparger: set of 8 holes in rectangular configuration:(1) 1 set located in the center, (2) 1 set located off center, (3) 3 sets located centrally Rectangular column: W = 0.15 m, D = 0.15 m, H = 1.0 m. Sparger: perforated plate do = 0.001 m, 49 holes Rectangular column: W = 0.15 m, D = 0.15 m, H = 1.0 m. Sparger: perforated plate do = 0.001 m, 49 holes Cylindrical column: D = 0.38 m, H = 3 m. Sparger: sintered bronze plate with avg. pore size of 50 ␮m Semi batch, VG = 0.0011 m/s Std. k–ε, low Reynolds number k–ε (a) NC 0.5 Semi batch, VG = 0.0013 m/s Std. k–ε Cd = 0.66 NC NC Number of grids (W × D × H) 27 × 10 × 52 cells 38 × 18 × 52 cells Equal size grid of 0.005 m, t = 0.1 s Semi batch, VG = 0.005 m/s Std. k–ε (b) NC NC Semi batch, VG = 0.005 m/s Std. k–ε LES Cs = 0.10 (b) 0.5 0.5 Semi batch, VG = 0.023 m/s Std. k–ε (c) NC NC Cylindrical column: D = 0.288 m, H = 2.6 m. Sparger: (1) perforated plate (do = 0.0007 m, 21 holes) (2) ring (do = 0.0007 m, 20 holes) Cylindrical column: D = 1 m, H = 2 m. Sparger: (1) single ring (ring dia = 0.45 m), (2) double ring (ring dia = 0.45, 0.78 m) Rectangular column: W = 0.2 m, D = 0.05 m, H = 1.2 m. Sparger: (1) sintered sparger, (2) four multipoint ones having 8 holes (do = 0.0008–0.002 m) Convergent channel is divided at bottom with a splitter plate, and each side is supplied independently with mixture of bubbles (dB = 3 mm) and water. W = 0.30 m H = 0.60 m D = 0.04 m (truncated) Std. k–ε Semi batch, VG = 0.0015, 0.005, 0.01 and 0.02 m/s Semi batch, Std. k–ε VG = 0.01, 0.02 and 0.03 m/s Std. k–ε Semi batch, VG = 0.0016–0.0083 m/s Cd = 0.44 NC NC (d) NC NC (e) 0.5 0.5 (b) 0.5 0–0.5 Pfleger et al. [24] Deen et al. [34] Deen et al. [36] Krishna and Van Baten [32] Pfleger and Becker [27] Ranade and Tayalia [38] Buwa and Ranade [28] Lakehal et al. [41] Inlet 0.22 m/s on one side, 0.54 m/s in other side. ∈g,in = 0.019 LES Tqo models. (a) Cs = 0.12, (b) DSM Model Number of grids (W × D × H) 15 × 50 × 160 cells, t = 0.01 s x = y = z = 0.010 m, t = 0.01/0.00 s (r, z, ␪), 30 × 160 × 20, t, 0.0005 s (100 steps) 0.001 s (100 steps) 0.005 s (19800 steps) Number of Grids 60 × 5 × 9, t = 0.1 s 2-grids, r = 23 and 46 cells, z = 24 and 48 cells, θ = 44 and 88 cells, t = 0.0001 s Number of grids (W × D × H) 7 × 7 × 25, 32 × 11 × 47, 61 × 19 × 92, t = 0.001 s, 0.01 s Meshsize varied as function of bubble diameter. x = y = z = = 0.0042 m (1.4 × dB ) = 0.0045 m (1.5 × dB ) = 0.0048 m (1.6 × dB ), t = 0.001 s Bove et al. [42] Rectangular column: W = 0.15 m, D = 0.15 m, H = 1.0 m. Sparger: perforated plate, do = 0.001 m, 49 holes Semi batch, VG = 0.005 m/s VLES Cs = 0.12 (b), (f) Bombardelli et al. [43] Rectangular column: W = 0.5 m, D = 0.15 m, H = 0.08 m Semi batch, VG = 0.0016 m/s LES-NWM In a mixture model framework Forces get cancelled under the assumption of dilute plume hypothesis and use of mixture equation. 0.5 0.5 3 different mesh sizes x = y = z, 0.010, 0.015 and 0.025 m, t = 0.005 s (100 s) x = y = z = 0.010 m, ∆t = 0.1 s M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Author 595 The total interfacial force acting between the two phases may arise from several independent physical effects: MI,L = −MI,G = MD,L + ML,L + MVM,L + MTD,L CD = max  24 0.687 ), 8 Re (1 + 0.15 Re 3  Eo Eo+4  (Tsuchiya et al. [61]); (f) CD = 4 ρL −ρG 1 3 ρL gdB V 2 , T where VT is the terminal rise velocity as defined by Tomiyama et al. [45].; (g) CD = ρL − ρG g. Vs where Vb is the rise velocity of bubble; (d) MD,L = −5 × 104 ǫG ǫL (uG − uL ); (e) if 24 0.687 )ǫ−1.7 , l Re (1 + 0.15 Re (a) CD = Dhotre et al. [47] Zhang et al. [44] Kulkarni et al.[3] Cylindrical column: D = 0.15 m, H = 1 m. Sparger: (1) single point (do = 0.00317 m), (2) multipoint spargers (d0 = 0.00196 m, 26 holes) Rectangular column: W = 0.15 m, D = 0.15 m, H = 1.0 m. Sparger: perforated plate do = 0.001 m, 49 holes Rectangular column: W = 0.15 m, D = 0.15 m, H = 1.0 m. Sparger: perforated plate do = 0.001 m, 49 holes 1 Re < 1000; (b) CD = 23 Eo 2 (Ishii and Zuber [46]); (c) CD = Semi batch, VG = 0.005 m/s LES Smagorinsky Cs varied as 0.08, 0.10, 0.15 and 0.2 LES Smagorinsky Cs = 0.12 and Dynamic Smagorinsky LES Semi batch, VG = 0.005 m/s Semi batch steady state, VG = 0.02 m/s Std. k–ε b 1 4 ρL −ρG 3 ρL gdB V 2 , x = y = z = 0.01 m. t = 0.005 s. Total time = 150 s x = y = z = 0.01 m. t = 0.005 s. Total time = 100 s 0.5 and Tomiyama [45] 0.5 (b) and Tomiyama [45] (b) 0.5 and 0.68 Tomiyama [45] 0.5 0.5 (g) 0.5 r = 22, θ = 42, z = 62 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 (7) The forces indicated above represent the interphase drag force, lift force, virtual mass force and turbulent dispersion force, respectively. In the present hydrodynamic model all forces except the virtual mass force has been used. According to Hunt et al. [49], the contribution of virtual mass force becomes negligible for column diameter greater that 0.15 m. Thakre and Joshi [50], Deen et al. [36] and Sokolichin and Eigenberger [10] have also pointed out the negligible effect of virtual mass force. A brief description of each interfacial force component is presented below. The origin of the drag force is due to the resistance experienced by a body moving in the liquid. Viscous stress creates skin drag and pressure distribution around the moving body creates form drag. The later mechanism is due to inertia and becomes significant as the particle Reynolds number becomes larger. The interphase momentum transfer between gas and liquid due to drag force is given by: 3 CD MD,L = − ǫG ρL |uG − uL |(uG − uL ) 4 dB (8) where CD is the drag coefficient taking into account the character of the flow around the bubble, and dB is the bubble diameter. The lift force arises from the net effect of pressure and stress acting on the surface of a bubble. The lift force in terms of the slip velocity and the curl of the liquid phase velocity can be described as: ML,L = CL ǫG ρL (uG − uL ) × ∇ × uL (9) where CL is the lift coefficient. The sign of this force depends on the orientation of slip velocity with respect to the gravity vector. The turbulent dispersion force, derived by Lopez de Bertodano [51], is based on the analogy with molecular movement. It approximates a turbulent diffusion of the bubbles by the liquid eddies. It is formulated as: MTD,L = −MTD,G = −CTD ρL k∇ǫL (10) where k is the liquid turbulent kinetic energy per unit of mass. CTD is the turbulent dispersion coefficient and its recommended range is between 0.1 and 0.5 [51]. For the case of LES, the Favr-averaged turbulent dispersion force has been used. 3.1. Turbulence equations 3.1.1. k–ε turbulence model When k–ε model is used, the turbulent eddy viscosity is formulated as follow µT,L = ρL Cµ k2 ε (11) 596 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 The turbulent kinetic energy (k) and its energy dissipation rate (ε) are calculated from their governing equations: ∂ (ρL ǫL k) + ∇(ρL ǫL uL k) ∂t   µeff,L ∇k + ǫL (G − ρL ε) = −∇ ǫL σk (12) ∂ (ρL ǫL ε) + ∇(ρL ǫL uL ε) ∂t   µL,eff ε = −∇ ǫL ∇ε + ǫL (Cε1 G − Cε2 ρL ε) σε k (13) The model constants are C␮ = 0.09; σ k = 1.00; σ ␧ = 1.00; Cε1 = 1.44, Cε2 = 1.92. The term G in above equation is the production of turbulent kinetic energy and described by: G = τL : ∇uL (14) 3.1.2. Reynolds Stress Modeling (RSM) turbulence model In the RSM model, individual Reynolds stresses u′i u′j are computed via a differential transport equation. The exact form of Reynolds stress transport equations is derived by taking moments of exact momentum equation. Thus, the RSM model solves six Reynolds stress transport equations. Along with these, an equation for dissipation rate is also solved. The exact transport equations for the transport of Reynolds stresses ρu′i u′j are given by: ∂ ∂ (ǫL ρL u′i u′j ) + (ǫL ρL uk u′i u′j ) ∂t ∂xk    ′ ′ ∂ k2 ∂ui uj 2 ′ = ǫL Pij + ǫL φij + ǫL µL + Cs′ ρ ∂xk 3 ε ∂xk 2 − δij ∈ L ρL ε 3 ∂ ∂ (ρL ǫL ε) + (ρL ǫL ui ε) ∂t ∂xi     ∂ µT,L ∂ε = ǫL µ L + ∂xj σε ∂xj   ∂ui ε ε2 − Cε2 ρL ǫL (17) +ǫL Cε1 ρL u′i u′k ∂xk k k k is calculated from the solved values of normal stress using the Reynolds stress transport equation, as ⎞ ⎛ 1 ⎝  ′ ′⎠ k= (18) ui ui 2 i=1,2,3 3.1.3. Large Eddy Simulation turbulence model Equations for LES are derived by applying filtering operation to the Navier–Stokes equations. The filtered equations are used to compute the dynamics of the large-scale structures, while the effect of the small-scale turbulence is modeled using a SGS model. Thus, the entire flow field is decomposed into a large-scale or resolved component and a small-scale or subgridscale component. The most commonly used SGS models are Smagorinsky model [52] and Dynamic Smagorinsky model, DSM, by Germano et al. [53]. In this work, the Smagorinsky model has been used. In case of LES, the velocities (u) in continuity equations and momentum equations represent the resolved velocities or grid scale velocities. (19) u = uinst − uSGS In case of LES Smagorinsky model, the turbulent eddy viscosity is formulated as follows: µT,L = ρL (Cs ∆)2 |S|2 (15) where φij is the pressure–strain correlation, and P′ , the exact production term, is given by: P ′ = −ρ u′i u′j (∇u)T + (∇u) u′i u′j transport equation: (16) As the turbulence dissipation appears in the individual stress equations, an equation for ε is computed with the model (20) where Cs is the Smagorinsky constant and S is the strain rate. The value of Smagorinsky constant used in this work is 0.10. 4. Numerical details The 3D transient simulations of flow pattern in two different cylindrical bubble columns were carried out using the commercial software ANSYS-CFX-10.0. The details of geometry, operating conditions, and mesh size used are given in Table 2. Table 2 Summary of geometries used for numerical simulations Author Menzel et al. [1] Bhole et al. [2] Geometry details Operating details, VG (m/s) D = 0.60 m H = 5.44 m 0.012 0.024 0.048 0.096 D = 0.15 m H=1m 0.02 Number of nodes used for simulations k–ε RSM 90000 90000 36000 36000 LES – 150000 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 597 The conservation equations were discretized using the control volume technique. The velocity–pressure linkage was handled through the SIMPLE procedure. The hybrid-upwind discretization scheme was used for the convective terms. A time step of 0.05 s was used for the k–ε and RSM simulations, whereas 0.005 s was used for LES simulations. The flow was simulated for 150 s and the data was time averaged over last 130 s. 4.1. Material balance The drift-flux model of Zuber and Findlay [54] is given by the following equation: VG /ǭG = C0 (VG + VL ) + C1 (21) where C0 = ǫG (VG + VL ) ǫG VG + VL  (22) ǫG ǫL VS  ǫG  (23) and C1 = The parameters C0 and C1 are the drift-flux constants. C0 represents the hold-up profile and C1 the bubble rise velocity. The most fortunate characteristic feature of this model is that the values of C0 and C1 are practically independent of the column diameter (of course when D > 150 mm and the sparger region is exceeded). Therefore, for a given gas–liquid system, a few measurements of ǭG with respect to VG and VL (over the range of interest) in a small diameter column (∼150 mm) enable the estimation of C0 and C1 . It is important that Eq. (21) holds for an extreme case of homogeneous regime having C0 = 1. 4.2. Energy balance Fig. 1. A typical radial grid layout for a geometry: (A) Menzel et al. [1] [61000 node], (B) Bhole et al. [2]. Side view layout [37000 node]. A typical grid layout is shown in Fig. 1. The grid independence study has been carried out using two grid resolutions for the std k–ε and Reynolds stress transport models. The grid sizes used for the independence system in simulating Menzel’s experiments were 90,000 and 120,000. For LES, a maximum mesh size of around 100 times the Kolmogorov length scale was used, which worked out to be about 3 mm at the centre. The Kolmogorov length scale was calculated using integral energy dissipation results of the k–ε model. While near the wall, the mesh size is around 0.1 mm. The gas inlet through the spargers was incorporated by creating mass source points at the specified position to mimic the exact sparger. Based on the superficial gas velocity (0.02 m/s), the mass flow rate was specified for each source point. Along the walls, no-slip boundary conditions were adopted. At the outlet of the column, the atmospheric pressure was specified as boundary condition. All the predicted flow patterns must satisfy the energy balance. The rate of energy supply from the gas phase to the liquid phase is given by the following Eq. (12): EI = π 2 D (ρL − ρG )gHD ǫL [VG + (CB − 1)ǭG VS ] 4 (24) When bubbles rise, the pressure energy is converted into turbulent kinetic energy. A fraction of CB is considered to get transferred to the liquid phase; the rate of energy given by Eq. (24) is finally dissipated in the turbulent liquid motion. From the turbulence models used for the prediction of flow pattern, we get the radial and axial variation of ε (energy dissipation rate per unit mass). From this ε field, the total energy dissipation rate can be calculated by suitable volume integration. The total energy dissipation rate must equal the energy-input rate given by Eq. (24). The pertinent detailed discussion has been provided by [12]. 5. Results and discussion For accurate prediction of local hydrodynamics, it is extremely important to properly select the simulation parameters of the interfacial forces like lift force, drag force, dispersion 598 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Table 3 Drag laws considered Author Model Fig. 2. Variation of slip velocity with bubble size for different drag law. () Schiller and Naumaan [62]; () Dalle Ville [63]; () Grace et al. [64]; (×) Ishii and Zuber [46]; () Ma and Ahmadi [56]; () Grevskott et al. [65]; (+) Zhang and Vanderheyden [57]. force and virtual mass force. These selections should always be made based on the considerations of actual physics. For example, a proper choice of drag law and bubble diameter is needed to accurately predict the slip velocity. So, it becomes extremely important that one understands the interrelation between drag force, bubble size and slip velocity. An attempt has been made here to elucidate this point. The slip velocity (VS ), which can be considered as the signature of the multiphase system under a given flow condition, is given as:    4dB ρL − ρG g (25) VS = 3CD ρL From Eq. (25), it is clearly seen that, for a given value of drag force, slip velocity changes with bubble size. To understand this interrelation between drag force, bubble size and slip velocity, slip velocity has been plotted as a function of bubble size for different drag laws (Fig. 2). The list of drag laws considered and its expressions have been given in Table 3. From Fig. 2, it can be seen that a single value of slip can be obtained from several combinations of drag law and bubble size. Similarly, for a particular bubble size, one can obtain several values of slip depending on the drag law. It is therefore of prime importance to choose the correct combination of drag law and bubble size to model the gas–liquid flow in a bubble column. Furthermore, with changes in superficial gas velocity, sparger design, and the change in nature of gas liquid systems, the average bubble diameter changes, and hence, the value of slip velocity also changes. Therefore, for different superficial velocities, either one has to 24 (1 + 0.15 Re0.687 ), if ReB < 1000 B ReB CD = 0.44 if ReB > 1000 2 4.8 = 0.63 + √ ReB 4 gdB (ρL − ρG ) = 3 VT2 ρL 2 0.5 = Eo 3 24 = (1 + 0.1 Re0.75 B ) ReB 5.645 = Eo−1 + 2.835 6 24 + = 0.44 + √ ReB 1 + ReB CD = Schiller and Naumaan [62] Dalle Ville [63] CD Grace et al. [64] CD Ishii and Zube [46] CD Ma and Ahmadi [56] CD Grevskott et al. [65] CD Zhang and Vanderheyde [57] CD VT is the terminal velocity as defined by Grace et al. [63]. take different bubble sizes, if drag law is constant, or select different drag laws, while keeping the bubble size constant. First option is more realistic and represents the actual physical picture. Hence, it is always advisable to change the bubble size with change in superficial gas velocity or change in sparger design. This elucidates the issue of proper choice of combination of drag law and bubble size as far as CFD simulation of bubble column hydrodynamics is concerned. Now, the drag force alone will not be sufficient to predict the local hydrodynamics correctly. Proper description of other interface forces like, lift force, dispersion force, are also important. Therefore, to understand the effect of different interphase forces, a series of simulations have been carried out. For this purpose, the experimental data reported by Menzel et al. [1] have been considered, and the data reported for lowest VG (0.012 m/s) and highest VG (0.096 m/s) have been taken into account for all the simulations. 5.1. Effect of drag law Simulations with all the drag laws given in Table 3 have been carried out. The bubble sizes for all simulations have been chosen in such a way so as to satisfy the average gas hold-up. The lift coefficient for the corresponding bubble size has been taken into account as suggested by Kulkarni [55]. The value reported by Kulkarni [55] follows the same trend as reported by Tomiyama [45], but the absolute value is less in the bubble size range of 6–8 mm. Standard parameter settings are given in Table 4. Simulation results for all the drag force have been given in Table 5. Table 4 Standard parameter setting for drag law (Zhang and Vanderheyden [56]) sensitivity run VG [exp.] (m/s) Bubble diameter, dB (mm) Lift force coefficient (CL) Turbulent dispersion force coefficient (CTD ) Added mass force coefficient 0.012 0.096 6 9 −0.06 −0.2 0.2 0.2 – – 599 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Table 5 Effect of drag law Author VG (m/s) ǭG vL (centerline) (m/s) Material balancea Energy balanceb LHS RHS LHS RHS 0.013 (0.012) 0.107 (0.096) 0.024 (0.029) 0.16 (0.128) 0.33 (0.29) 1.01 (0.86) 0.41 0.75 0.44 0.63 0.025 0.34 0.031 0.43 0.013 (0.012) 0.11 (0.096) 0.026 (0.029) 0.165 (0.128) 0.36 (0.29) 1.22 (0.86) 0.41 0.75 0.44 0.60 0.025 0.28 0.027 0.42 Grace et al. [64] 0.013 (0.012) 0.11 (0.096) 0.036 (0.029) 0.16 (0.128) 0.32 (0.29) 1.21 (0.86) 0.41 0.75 0.32 0.62 0.027 0.3 0.039 0.46 Ishii and Zuber [46] 0.012 (0.012) 0.11 (0.096) 0.034 (0.029) 0.15 (0.128) 0.33 (0.29) 0.928 (0.86) 0.41 0.75 0.35 0.68 0.027 0.33 0.041 0.44 Ma and Ahmadi [56] 0.013 (0.012) 0.102 (0.096) 0.024 (0.029) 0.11 (0.128) 0.33 (0.29) 1.03 (0.86) 0.41 0.75 0.47 0.81 0.025 0.3 0.03 0.42 Grevskott et al. [65] 0.013 (0.012) 0.01 (0.096) 0.035 (0.029) 0.16 (0.128) 0.22 (0.29) 0.92 (0.86) 0.41 0.75 0.29 0.65 0.024 0.32 0.04 0.48 Zhang and Vanderheyden [57] 0.013 (0.012) 0.105 (0.096) 0.028 (0.029) 0.14 (0.128) 0.29 (0.29) 0.88 (0.86) 0.41 0.75 0.38 0.65 0.028 0.036 0.036 0.46 Schiller and Naumaan [62] Dalle Ville [63] The bracketed numbers indicate the experimental values. a LHS = V /ǭ ; RHS = C (V̄ + V̄ ) + C . G G 0 G L 1 b LHS = volume integral of energy dissipation, RHS = energy input. It can be seen from Table 5 that though all the drag laws have been able to predict the average gas hold-up and the centerline velocity quite well, the drag law reported by Ishii and Zuber [46], Ma and Ahmadi [56] and Zhang and Vanderheyden [57] is more closer to the experimental values. Further, axial liquid velocity profile and gas hold-up profile for these three drag laws have been compared with the experimental profile and the results are shown in Fig. 3. For low VG , all the three drag laws give good agreement with the experimental value, but at the higher value of VG , the drag law of Zhang and Vanderheyden [57] was found to give better agreement with experimental value. The drag law of Ishii and Zuber [46] over predicts the hold-up profile. This can be attributed to the fact in this drag law that the slip velocity remains constant for the entire range of bubble size (Fig. 2), which ultimately leads to the over-prediction of hold-up profile at higher VG . The slip velocity calculated from the drag law of Ma and Ahmadi [56] is greater than the slip velocities calculated from all the other drag law in the range of bubble size greater than 5 mm (Fig. 2). For drag law of Ma and Ahmadi [56], the bubble size of 5 mm and 6 mm was found to give average gas hold-up close to experimental value for superficial gas velocity of 0.012 m/s and 0.096 m/s, respectively. The higher slip velocity value in this range therefore justifies the greater axial liquid velocity, which also leads to the steeper hold-up profile. Therefore, the drag law of Zhang and Vanderheyden [57], which gives better prediction for both the VG , has been used for further simulations. 5.2. Effect of lift force To study the effect of lift force, two more values of CL were chosen for both the cases of VG along with the value of lift coefficients given in Table 4. Simulations have been carried out with the value of CL = 0 and the corresponding positive val- Fig. 3. Effect of drag law on: (A) average liquid velocity; (B) gas hold-up. () Experimental data of Menzel et al. [1] [VG = 0.012 m/s], () experimental data of Menzel et al. [1] [VG = 0.096 m/s]; (1) Ishii and Zuber [46], (2), Ma and Ahmadi [56], (3) Zhang and Vanderheyden [57]. 600 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Fig. 4. Effect of lift force, (A) average liquid velocity; (B) gas hold-up. () Experimental data of Menzel et al. [1] [VG = 0.012 m/s], () experimental data of Menzel et al. [1] [VG = 0.096 m/s]; (1) CL = 0 (2) CL = −ve (3) CL = +ve. ues (CL = 0.06 and CL = 0.2), while other parameter values were same as that reported in Table 4. The results for axial liquid velocity profile and gas hold-up profile have been given in Fig. 4. For all values of lift force coefficients, it is seen that deviation in results at lower VG is not that significant as compared to that observed at higher VG . The positive value of CL makes the bubbles move outwards towards the column wall, which leads to a flatter hold-up profile and lower centerline velocity. Therefore, the value of lift coefficient should be chosen depending on the bubble size [unlike several researchers, who have taken constant CL (see Table 1)] and the values reported by Kulkarni [55] were found to give good agreement. 5.3. Effect of turbulent dispersion force Simulations were carried out with three values of CTD (0, 0.2 and 0.5), and all the other parameters were kept same as reported in Table 4. The results of average axial liquid velocity profile and gas hold-up profile for superficial gas velocity of 0.012 m/s and 0.096 m/s have been given in Fig. 5. It can be seen from Fig. 5 that at low VG , the effect of CTD is not that significant, but at higher VG , increase in value of CTD makes the hold-up profile comparatively flatter. Although, CTD value of 0.2 was found Fig. 5. Effect of turbulent dispersion force, (A) average liquid velocity; (B) gas hold-up. () Experimental data of Menzel et al. [1] [VG = 0.012 m/s], () experimental data of Menzel et al. [1] [VG = 0.096 m/s]; (1) CTD = 0 (2) CTD = 0.2 (3) CTD = 0.5. to give a good agreement, but we feel the choice of the value CTD is intuitive, and any value between 0 and 0.5 can be taken [10], depending on the system under consideration, which can predict the hold-up and axial liquid velocity profile closer to the experimental value. 5.4. Effect of added mass force The added mass force was neglected in all the previous simulations in accordance with the observation made by Hunt et al. [49], Thakre and Joshi [50], Deen et al. [36] and Sokolichin et al. [8]. This observation gets further strengthened by two simulations carried out at CVM = 0.5 and CVM = 0. The results are given in Fig. 6. True to the previous observation, it can be seen that, the addition of added mass force has no significant effect on results. With this understanding, simulations have been carried out for the experimental data reported by Menzel et al. [1] for four different VG (Table 2), so as to validate the CFD model for a wide range of superficial gas velocity. The numerical details and the results of the simulations are given in Table 6. The comparison of predicted average liquid velocity and gas hold-up profile with 601 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Fig. 6. Effect of added mass force, (A) average liquid velocity; (B) gas hold-up. () Experimental data of Menzel et al. [1] [VG = 0.012 m/s], () experimental data of Menzel et al. [1] [VG = 0.096 m/s]; (1) CVM = 0 (2) CVM = 0.5. experimental value are shown in Fig. 7. The CFD prediction shows a good agreement with the experimental results. 5.5. Comparison of turbulence models After understanding the sensitivity of interphase forces, we decided to consider the effect of turbulence closer on CFD predictions. Amongst various closure models, the CFD simulations with k–ε turbulence model are plenty, and recently attention has also been given to LES turbulence model (see Table 1). On the other hand, no results have been available with RSM tur- Fig. 7. Comparison between the CFD prediction and experimental data of Menzel et al. [1]. (A) average liquid velocity; (B) gas hold-up. () VG = 0.012 m/s; () VG = 0.024 m/s; () VG = 0.048 m/s; () VG = 0.096 m/s. bulence model. To understand the relative behavior of different turbulence models, the predictive performance of k–ε, RSM and LES turbulence models have been compared. These models are applied to simulate the hydrodynamics (global as well as local) of a cylindrical bubble column arising due to the use of three different spargers (a sieve plate sparger, a single hole and a sintered plate sparger) at the superficial gas velocity of 20 mm/s. The predictions of axial velocity, gas hold-up, turbulent kinetic energy and Reynolds stress profiles have been compared with the experimental data of Bhole et al. [2] and kulkarni et al. [3]. Table 6 Numerical details and comparison of CFD predictions with the data of Menzel et al. [1] Numerical details of Menzel simulation VG [exp] (m/s) 0.012 0.024 0.048 0.096 Simulation results observed Bubble diameter, dB (mm) Lift 6 7 8 9 −0.06 −0.09 −0.15 −0.2 Turbulent dispersion force VG (m/s) 0.2 0.2 0.2 0.2 0.013 0.027 0.05 0.105 The bracketed numbers indicate the experimental values. a LHS = V /ǭ ; RHS = C (V̄ + V̄ ) + C . G G 0 G L 1 b LHS = volume integral of energy dissipation, RHS = energy input. Hold-Up 0.028 (0.029) 0.050 (0.048) 0.104 (0.103) 0.14 (0.128) Material balancea Energy balanceb LHS RHS LHS RHS 0.41 0.50 0.47 0.75 0.38 0.49 0.44 0.65 0.028 0.092 0.185 0.36 0.036 0.11 0.208 0.46 602 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Fig. 8. Comparison between the simulated and experimental profiles of axial liquid velocity at different axial positions in a 150 mm (i.d.) bubble column with sieve plate sparger at VG = 20 mm/s (A) H/D = 1; (B) H/D = 2; (C) H/D = 3; (D) H/D = 4. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. The comparison is presented at four axial locations corresponding to height to diameter ratio of 1 to 4 (z = 0.15 m, 0.3 m, 0.45 m and 0.6 m). Fig. 8 shows the comparison for radial profile of axial velocity for sieve plate bubble column. All the models show that axial flow is distinctly upward in the central region with higher gas hold-up, while a downward counter flow is observed in the near wall region with low gas hold-up. This hold-up gradient is creating the density difference for liquid circulation to take place. The point of flow reversal is clearly seen at a radial location of around r/R = 0.6–0.8 m. It may be noted that near the sparger, none of the models have been able to capture the flow pattern as observed in experiments. At sparger, there is a transient variation of bubbling area, which leads to oscillations and instabilities, and hence near sparger the experimentally observed flow profile is flatter and does not satisfy the material balance at that cross-section. The models have not been able to capture the effect of random variation of bubble formation across the 25 hole of sieve plate; hence we observe that these turbulence models give average axial velocity profiles that nearly satisfy the cross-sectional material balance near sparger. The predictive performance of these turbulence models improves at higher axial location. Fig. 9 shows the comparison for radial profile of axial velocity for sintered plate bubble column. The turbulence models are able to simulate the nearly homogeneous conditions observed in sintered columns. The flat hold-up profiles and flat axial velocity profiles are obtained, with flow reversals occurring at region greater than r/R = 0.8. Fig. 10 shows the comparison for radial profile of axial velocity for single hole sparger bub- ble column. LES is showing better agreement with experiments as compared to the other two RANS based models in all the three cases. This could be because of its ability in capturing the transient central bubble plume movement. Though, RSM is expected to give better results for anisotropic flows involving swirls, acceleration and deceleration, and buoyancy, as the case is in bubble column, the profiles show that RSM is marginally better in predicting the averaged axial liquid velocity profiles than k–ε. Figs. 11–13 show the comparison of the turbulence models in simulating the radial profile of gas hold-up (ǫG ) for sieve plate bubble column, sintered plate bubble column and single hole bubble column, respectively. For sintered bubble column, the simulation gives a flatter and higher averaged gas hold-up profile as compared to that in sieve plate and single hole bubble column. Further, at any radial location, the value of ǫG is higher for the sintered plate than the sieve plate, and further higher for the sieve plate than the single point sparger. All these features get fairly well predicted by all the three models. Figs. 14–16 show the comparison for radial profile of turbulent kinetic energy for sieve, sintered plate sparger and single hole bubble column, respectively. In all the three cases, we observe that k–ε model is not able to predict the turbulent kinetic energy anywhere near to the experimentally obtained values, and the RSM and the LES show comparatively better agreements with the experimental data, with slight over-prediction observed with LES. These deviations could be explained on the basis of how a particular model tries to capture the energetic interactions between the mean flow and the large scale, and the subsequent M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 603 Fig. 9. Comparison between the simulated and experimental profiles of axial liquid velocity at different axial positions in a 150 mm (i.d.) bubble column with sintered plate sparger at VG = 20 mm/s (A) H/D = 1; (B) H/D = 2; (C) H/D = 3; (D) H/D = 4. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. Fig. 10. Comparison between the simulated and experimental profiles of axial liquid velocity at different axial positions in a 150 mm (i.d.) bubble column with single hole sparger at VG = 20 mm/s (A) H/D = 0.5; (B) H/D = 2.5; (C) H/D = 3.4; (D) H/D = 4.6. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. 604 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Fig. 11. Comparison between the simulated and experimental profiles of gas hold-up at different axial positions in a 150 mm (i.d.) bubble column with sieve plate sparger at VG = 20 mm/s (A) H/D = 1; (B) H/D = 2; (C) H/D = 3; (D) H/D = 4. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. energy cascade process from large scale to small scale. In RSM, there is the pressure strain mechanism which does not add or destruct any turbulent kinetic energy but helps to redistribute the energy between the normal components ensuring accuracy. While in k–ε model, which is based on isotropic assumption, the normal stresses get poorly represented, as the turbulent kinetic energy (k) formulation is constructed with a limitation that all the normal components of stresses are equal to each other. This Fig. 12. Comparison between the simulated and experimental profiles of average gas hold-up at different axial positions in a 150 mm (i.d.) bubble column with sintered plate sparger at VG = 20 mm/s (A) H/D = 1; (B) H/D = 2; (C) H/D = 3; (D) H/D = 4. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 605 Fig. 13. Comparison between the simulated and experimental profiles of average gas hold-up at different axial positions in a 150 mm (i.d.) bubble column with single hole sparger at VG = 20 mm/s (A) H/D = 0.5; (B) H/D = 2.5; (C) H/D = 3.4; (D) H/D = 4.6. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. Fig. 14. Comparison between the simulated and experimental profiles of turbulent kinetic energy at different axial positions in a 150 mm (i.d.) bubble column with sieve plate sparger at VG = 20 mm/s (A) H/D = 1; (B) H/D = 2; (C) H/D = 3; (D) H/D = 4. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. 606 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Fig. 15. Comparison between the simulated and experimental profiles of turbulent kinetic energy at different axial positions in a 150 mm (i.d.) bubble column with sintered plate sparger at VG = 20 mm/s (A) H/D = 1; (B) H/D = 2; (C) H/D = 3; (D) H/D = 4. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. results in k getting equally distributed, thus leading to inaccurate kinetic energy profiles. In the case of LES, most of energy containing scales and the associated mechanism of energy transfer are resolved through the use of finer mesh size and the transfer of energy to unresolved scales is done by using sub-grid stresses. Figs. 17 and 18 show the comparison for radial profile of Reynolds stress (axial-radial) for sieve and sintered plate sparger bubble column, respectively. For single hole bubble column, experimental details on Reynolds stress is not available. In case of sintered plate, the Reynolds stress profiles are seen to be flatter and lower in magnitude as compared to sieve plate profile. For sieve plate bubble column the magnitude of Reynolds stress increase from centre to reach a maxima and then decreases towards the wall. The observed values of shear stress can be explained on the basis of presence of dynamic structures [58]. In case of sintered plate bubble columns, the scales are of the size limited to bubble-bubble inter-spacing or bubble wakes. These structures are much smaller than the equipment size and therefore do not have any preferential location. Hence, for sintered plate bubble column, where large-scale dynamic structures are absent and nearly homogeneous flow profile exists, we observe very low mean shear stress values. In the case of sieve plate bubble column, flow structures of the size 5–10 mm are observed in LES simulations. Though these structures are not of the largescale highly dynamic kind, but they do meander slightly and contribute to mean shear stress. Thus, we observe higher mean shear stress values as compared to that in sintered plate. The axial radial stress profiles as obtained for sintered and sieve plate shows that the RSM has been able to show better agreement with the experiments as compared to the LES, which can be seen to overpredict. It is well known that the LES (with a single Smagorinsky constant) is generally unable to represent consistently the correct sub grid-scale stress in various flow situations [53]. This could probably be the reason as to why the LES model slightly over predicts the mean shear stress profile and even the kinetic energy profiles at some places. Despite this limitation with the modeled part, the advantage of LES has been its ability to resolve the flow structures greater than the mesh size used. The information obtained from LES has been used to gain insights into flow profiles, and is discussed below. 5.6. Flow information from instantaneous profiles of large eddy simulation Chen et al. [59] had suggested that the flow structures present within the various flow regions are instantaneous phenomena, which do not get captured when averaging procedures are used to quantify the flow properties. In this context, LES can be an effective tool to study the instantaneous flow profiles, but until now, practically no attempt has been made to harness the capacity of LES to understand the flow structures. Therefore, it was thought desirable to undertake a systematic attempt in this work. Chen et al. [59] had used the PIV technique and observed three flow regimes (dispersed bubble, vortical–spiral flow and tur- M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 607 Fig. 16. Comparison between the simulated and experimental profiles of turbulent kinetic energy at different axial positions in a 150 mm (i.d.) bubble column with single hole sparger at VG = 20 mm/s (A) H/D = 0.5; (B) H/D = 2.5; (C) H/D = 3.4; (D) H/D = 4.6. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. Fig. 17. Comparison between the simulated and experimental profiles of Reynolds stress at different axial positions in a 150 mm (i.d.) bubble column with sieve plate sparger at VG = 20 mm/s (A) H/D = 1; (B) H/D = 2; (C) H/D = 3; (D) H/D = 4. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. 608 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Fig. 18. Comparison between the simulated and experimental profiles of Reynolds stress at different axial positions in a 150 mm (i.d.) bubble column with sintered plate sparger at VG = 20 mm/s (A) H/D = 1; (B) H/D = 2; (C) H/D = 3; (D) H/D = 4. () Experimental; (1) LES model; (2) RSM model and (3) k–ε model. bulent flow regime) for various operating conditions in a 3D bubble column (Fig. 19). For the case of our sieve plate bubble column operating at a superficial gas velocity of 20 mm/sec, a vortical-spiral flow regime is expected. We observe that the instantaneous flow profiles obtained from LES, as shown in Figs. 20 and 21, has indeed been able to nearly capture the vortical-spiral flow regime and the four flow regions (descending flow region, vortical-spiral flow region, fast bubble region and the central plume region) within this regime. Fig. 20(A) has three snapshots taken at three different times, and each is showing instantaneous velocity vectors superimposed with grey scale contours based on instantaneous values of gas hold-up. Near the wall, the downward axial velocity profiles shows the descending flow region, and the higher gas hold-up at the centre shows the existence of the central bubble stream. Fig. 20(B) Fig. 19. Flow structure in a three-dimensional bubble column obtained from Fig. 6 of Chen et al. [59]. shows the snapshot zoomed at a given H/D for better view. Now, in between the descending flow region and the central bubble stream region lies, the vortical-spiral flow region that LES has been able to capture very well. Fig. 21 shows the top view of the flow pattern over a cross-sectional plane at the mid section of the column, where the vortical-spiral regime is clearly seen. Nearer and in-between these vortices, higher gas hold-up regions are observed. The dynamic interactions between the vortices and meandering central bubble stream can be observed by change in their positions with time. The simulation result shows that most of the dynamic vortices (eddies) are seen to have sizes of around 5–10 mm (size equivalent of about 1/2 to 1/3 rd of a H/D). The central bubble stream comprises of fast bubble region and the central plume region. But, LES could not distinguish between these two regions. The fact that whether these two regions are distinct, or whether they have merged to form a central bubbles stream, could not be ascertained with certainty. This is because we are working in the Eulerian LES framework without employing the bubble coalescence and break up models. Chan et al. [59] have opined that, a fast bubble region is characterized by dynamic bubble-bubble interactions with coalescence and break-up, while a central bubble plume region has uniform bubble size distribution. Hence, the knowledge of bubble size distribution, and developing an effective coalescence and break up model that works within the LES frame work could perhaps solve this problem. In spite of this limitation, the LES has been able to capture the vortical-spiral regime very effectively. In the case of sieve plate bubble column, Harteveld et al. [60] observed that bubble columns with uniform aeration give M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 609 Fig. 20. Snapshot of LES results showing instantaneous velocity vectors superimposed with grey scale contours based on instantaneous gas hold-up in a sieve plate bubble column at a plane in center of the column (A) at three different time steps of 12, 30 and 65 s. (B) Snapshot zoomed to give better view of flow for a region with height equivalent to a H/D. Fig. 21. Top view of flow pattern in a cross-sectional area at the middle of the column as obtained for a sieve plate bubble column using LES simulation. 610 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Fig. 22. Snapshot of LES results showing instantaneous velocity vectors superimposed with grey scale contours based on instantaneous gas hold-up in a sintered plate bubble column at a plane in center of the column (A) at three different time steps of 2.5, 30 and 65 s. (B) Snapshot zoomed to give better view of flow for a region with height equivalent to a H/D. very uniform gas distribution, and at a low superficial gas velocity, the familiar large-scale circulation and structures are absent. For our sintered bubble column, having geometry and operating flow conditions similar to Hartveld et al. [60], Fig. 22(A) shows the instantaneous velocity vectors superimposed with grey scale contours based on instantaneous gas hold-up, while Fig. 22(B) shows the snapshot zoomed at a given H/D for better view. At 2.5 s, the flow profile is not fully developed and a complete homogeneous dispersion of gas is not achieved. After some time, as seen at time of 30 s and 60 s, the instantaneous gas hold shows nearly homogeneous behavior, and dynamic large-scale coherent structures or eddies are not seen. Fig. 23 shows the top view of the flow pattern over a plane at the mid section of the column, where nearly homogeneous dispersal of gas is seen as compared to that in sieve plate as shown in Fig. 22, and also, though one can observe one or two eddies in the flow patterns, any distinct vortical-spiral region is absent. The LES has been able to simulate the flow conditions very well, probably because under nearly homogeneous conditions, bubble-bubble coalescence and break up effects are negligible. M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 611 Fig. 23. Top view of flow pattern in a cross-sectional area at the middle of the column as obtained for a sintered plate bubble column using LES simulation. Fig. 24. Snapshot of LES results showing instantaneous velocity vectors superimposed with grey scale contours based on instantaneous gas hold-up in a single hole bubble column at a plane in center of the column (A) at three different time steps of 10, 21 and 60 s. (B) Snapshot zoomed to give better view of flow for a region with height equivalent to a H/D. 612 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 Fig. 25. Top view of flow pattern in a cross-sectional area at the middle of the column as obtained for a single hole bubble column using LES simulation. In the case of single hole bubble column, the LES has been able to capture the bubble plume dynamics very well. The oscillations of central bubble plume, and the vortices surrounding the plume region are visible from Figs. 24 and 25. Fig. 24(A) shows the plume oscillations as simulated by LES at three different time steps (10 s, 21 s and 60 s). While, Fig. 24(B) shows the plume behaviour at a given H/D region. The top view of the flow pattern (in Fig. 25) also shows the eddies and the central bubble plume at the mid-plane of the column. The LES has been successful in capturing the instantaneous flow patterns arising owing to the variation of sparger design. Thus, it can be viewed as an important tool for study of coherent structures. 6. Conclusions CFD simulations have been performed to study the sensitivity of different interphase forces (drag, lift, turbulent dispersion, added mass) and different turbulence models (k–ε, RSM and LES). Conclusions of this study are as follows. The importance of proper choice of the combination of bubble size and drag law has been brought out. Seven different drag laws were studied, almost all the drag laws were able to predict the global hydrodynamics within certain limits. It is the prediction of local hydrodynamics where the drag law of Zhang and Vanderheyden [57] scores over the others. The effect of lift force was studied by choosing different value of CL , and the importance of choice of CL value according to bubble size has been highlighted. The CL value as a function of bubble size reported by Kulkarni et al. [55] was found to give better predictions. The effect of turbulent dispersion force on the simulation results has been studied. Although, we feel the choice of accurate value of CTD is intuitive, this activity has brought out the importance of the effect of CTD on the local hydrodynamics of bubble columns. As CTD increases, the hold-up profile becomes flatter. No significant contribution of added mass force on the simulation of local hydrodynamics of bubble column is seen, which is in accordance with the observation made earlier by Hunt et al. [49], Thakre and Joshi [50], Deen et al. [36] and Sokolichin and Eigenberger [8]. This is mainly because the acceleration and deceleration effects are restricted to small end regions of the column. The performance of three CFD models, viz k–ε, RSM and LES, has been compared with the experimental data of Bhole et al. [2]. Near the sparger, none of the turbulence model has been able to predict the axial velocity profiles and the gas holdup profiles as reported experimentally. The predictive capability improves at higher axial locations, where the flow gets developed. Contrary to the expectations, the RSM has not been able to show better predictive performance than k–ε in predicting the average axial velocity profiles. The profiles of Reynolds stress and kinetic energy for all the three turbulence models are in partial agreement with some deviations. In predicting turbulent kinetic energy, RSM has done a better job than k–ε, which could be due to its intrinsic ability in capturing the anisotropic energy transfer mechanism. LES has been successful in capturing averaged behavior of the flow. Though, at some locations, it slightly over predicts the kinetic energy and stress profiles, which could be due to the fact that with a single Smagorinsky constant model, the LES is generally unable to represent consistently the correct sub grid-scale stress for various flow situations. We feel that LES needs to be used in a more effective way as a tool to study the instantaneous flow to capture the coherent structures and to develop better turbulence models. On the whole, with RSM and LES simulations, there is very little gain in information at the cost of higher computational resources. Hence, the k–ε model can be preferred over the RSM and LES models for simulating 3D bubble column for getting average information. LES has been able to capture the instantaneous phenomena quite well. The LES shows the flow structures and the flow regions of the vortical-spiral flow regime in the sieve plate column, and the bubble plume dynamics along with vortical structures for the single hole sparger. Perhaps, the incorporation of bubble coalescence and break up model in LES could help to distinguish the two regions (fast bubble region and central bubble plume region) in the vortical-spiral regime, and thus help in better prediction of regions. For sintered columns, LES gives nearly homogeneous flow conditions with absence of dynamic eddies. LES can thus be effectively used for study of coherent structures and instantaneous flow profiles. Acknowledgement One of the authors would like to acknowledge the fellowship support given by the University Grants Commission (UGC), Government of India. References [1] T. Menzel, T. Weide, O. Staudacher, U. Onken, Reynolds stress model for bubble column reactor, Ind. Eng. Chem. Res. 29 (1990) 988–994. [2] M.R. Bhole, S. Roy, J.B. Joshi, Laser doppler anemometer measurements in bubble column: effect of sparger, Ind. Eng. Chem. Res. 45 (26) (2006) 9201–9207. [3] A.A. Kulkarni, K. Ekambara, J.B. Joshi, On the development of flow pattern in a bubble column reactor: experiments and CFD, Chem. Eng. Sci. 62 (2007) 1049–1061. [4] Y.T. Shah, B.G. Kelkar, S.P. Godbole, W.D. Deckwer, Design parameters estimations for bubble column reactors, AIChE J. 28 (3) (1982) 353–379. [5] L.K. Doraiswamy, M.M. Sharma, Heterogeneous Reactions, vol. 2, J. Wiley, New York, 1986. [6] W.D. Deckwer, Bubble Column Reactors, John Wiley and Sons, New York, 1992. M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 [7] J.B. Joshi, Computational flow modeling and design of bubble column reactors, Chem. Eng. Sci. 56 (2001) 5893–5933. [8] A. Sokolichin, G. Eigenberger, Simulation of buoyancy driven bubbly flow: established simplifications and open questions, AIChE J. 50 (2004) 24–45. [9] M. Rafique, P. Chen, M.P. Dudukovicı̌, Computational modeling of gas–liquid flow in bubble columns, Rev. Chem. Eng. 20 (2004) 225–375. [10] A. Sokolichin, G. Eigenberger, Gas–liquid flow in bubble columns and loop reactors. Part-I. Detailed modelling and numerical simulation, Chem. Eng. Sci. 49 (1994) 5735–5746. [11] R.F. Mudde, O. Simonin, Two and three dimensional simulation of a bubble plume using a two fluid model, Chem. Eng. Sci. 54 (1999) 5061–5071. [12] M.T. Dhotre, J.B. Joshi, Two-dimensional CFD model for prediction of pressure drop and heat transfer coefficient in bubble column reactors, Trans. Inst. Chem. Eng. 78 (2004) 689–707. [13] K. Ekambara, M.T. Dhotre, J.B. Joshi, CFD simulations of bubble column reactors: 1D, 2D and 3D approach, Chem. Eng. Sci. 60 (2005) 6733–6746. [14] A. Lapin, A. Lubbert, Numerical simulations of the dynamics of two-phase gas–liquid flows in bubble columns, Chem. Eng. Sci. 49 (1994) 3661–3674. [15] N. Devanathan, M.P. Dudukovic, A. Lapin, A. Lubbert, Chaotic flow in bubble column reactors, Chem. Eng. Sci. 50 (16) (1995) 2661–2667. [16] E. Delnoij, F.A. Lammers, J.A.M. Kuipers, W.P.M. Van Swaaij, Dynamic simulation of dispersed gas–liquid two-phase flow using a discrete bubble model, Chem. Eng. Sci. 52 (1997) 1429–1458. [17] E. Delnoij, J.A.M. Kuipers, W.P.M. Van Swaaij, A three-dimensional CFD model for gas–liquid bubble columns, Chem. Eng. Sci. 54 (1999) 2217–2226. [18] S. Lain, D. Broder, M. Sommerfeld, Experimental and numerical studies of the hydrodynamics in a bubble column, Chem. Eng. Sci. 54 (1999) 4913–4920. [19] V.V. Buwa, D.S. Deo, V.V. Ranade, Eulerian–Lagrangian simulations of unsteady gas–liquid flows in bubble columns, Int. J. Multiphase Flow 32 (7) (2006) 864–885. [20] A. Sokolichin, G. Eigenberger, Applicability of the standard turbulence model to the dynamic simulation of bubble columns. Part-I. Detailed numerical simulations, Chem. Eng. Sci. 54 (1999) 2273–2284. [21] H.A. Jakobsen, B.H. Sannæs, S. Grevskott, H.F. Svendsen, Modeling of vertical bubble-driven flows, Ind. Eng. Chem. Res. 36 (10) (1997) 4052–4074. [22] V.V. Ranade, Modeling of turbulent flow in a bubble column reactor, Chem. Eng. Res. Des. 75 (1997) 14–23. [23] O. Borchers, C. Busch, A. Sokolichin, G. Eigenberger, 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 (1999) 5927–5935. [24] D. Pfleger, S. Gomos, N. Gilbert, H.G. Wagner, Hydrodynamic simulation of laboratory scale bubble columns: fundamental studies of Eulerian– Eulerian modeling approach, Chem. Eng. Sci. 54 (1999) 5091–5099. [25] J. Sanyal, S. Vasquez, S. Roy, M.P. Dudukovic, Numerical simulation of gas–liquid dynamics in cylindrical bubble column reactors, Chem. Eng. Sci. 54 (1999) 5071–5083. [26] Y. Pan, M.P. Dudukovic, M. Chang, Dynamic simulation of bubble flow in bubble columns, Chem. Eng. Sci. 54 (1999) 2481–2489. [27] D. Pfleger, S. Becker, Modelling and simulation of the dynamic flow behaviour in a bubble column, Chem. Eng. Sci. 56 (2001) 1737–1745. [28] V.V. Buwa, V.V. Ranade, Dynamics of gas–liquid flows in rectangular bubble columns: experiments and single/multi-group simulations, Chem. Eng. Sci. 57 (2002) 4715–4736. [29] A. Lapin, T. Paaschen, K. Junghans, A. Lubbert, Bubble column fluid dynamics, flow structures in slender columns with large-diameter ringspargers, Chem. Eng. Sci. 57 (8) (2002) 1419–1424. [30] R.S. Oey, R.F. Mudde, H.E.A. Van den Akker, Sensitivity study on interfacial closure laws in two-fluid bubbly flow simulations, AIChE J. 49 (2003) 1621–1636. [31] S.M. Monahan, V.S. Vitankar, R.O. Fox, CFD predictions for flow-regime transitions in bubble columns, AIChE J. 51 (2005) 1897–1923. [32] R. Krishna, J.M. Van Baten, Eulerian simulations of bubble columns operating at elevated pressures in the churn turbulent flow regime, Chem. Eng. Sci. 56 (2001) 6249–6258. 613 [33] F. Bertola, M. Vanni, G. Baldi, Application of computational fluid dynamics to multiphase flow in bubble column, Int. J. Chem. Reactor Eng. 1 (2003) 1–13. [34] N.G. Deen, T. Solberg, B.H. Hjertager, Numerical simulation of gas–liquid flow in a square cross section bubble column, in: Presented at 14th International congress of chemical and process engineering (CHISA), Praha, Czech Republic, 2000. [35] S. Becker, A. Sokolichin, G. Eigenberger, Gas–liquid flow in bubble columns and loop reactors. Part II. Comparison of detailed experiments and flow simulations, Chem. Eng. Sci. 49 (1994) 5747–5762. [36] N.G. Deen, T. Solberg, B.H. Hjertager, Large eddy simulation of gas–liquid flow in a square cross-sectioned bubble column, Chem. Eng. Sci. 56 (2001) 6341–6349. [37] J.H. Hills, Radial non-uniformity of velocity and voidage in a bubble column, Trans. Inst. Chem. Eng. 52 (1974) 1–13. [38] V.V. Ranade, Y. Tayalia, Modelling of fluid dynamics and mixing in shallow bubble column reactors: influence of sparger design, Chem. Eng. Sci. 56 (4) (2001) 1667–1675. [39] M.P. Schwarz, W.J. Turner, Applicability of the standard k–ε turbulence model to gas-stirred baths, Appl. Math. Model. 12 (3) (1988) 273–279. [40] V.V. Buwa, V.V. Ranade, Dynamics of gas–liquid flow in 2D bubble columns: influence of sparger design, column diameter and physical properties, in: Presented at international conference of multiphase flows (ICMF)-2001, New Orleans, USA, 2001. [41] D. Lakehal, L.B. Smith, M. Milelli, Large-eddy simulation of bubbly turbulent shear flows, J. Turbulence 3 (2002) 1–21. [42] S. Bove, T. Solberg, B.H. Hjertager, Numerical aspects of bubble column simulations, Int. J. Chem. Reactor Eng. 2 A1 (2004) 1–22. [43] F.A. Bombardelli, G.C. Buscaglia, M.H. Garcı́a, E.A. Dari, Simulation of bubble-plume wandering phenomena in a bubble plume using a k–ε model and a large-eddy-simulation (LES) approach, Mecanica Computational XXIII (2006) 1969–1994. [44] D. Zhang, N.G. Deen, J.A.M. Kuipers, Numerical simulation of dynamic flow behavior in a bubble column: a study of closures for turbulence and interface forces, Chem. Eng. Sci. 61 (2006) 7593–7608. [45] A. Tomiyama, Drag lift and virtual mass forces acting on a single bubble, in: Third International Symposium on Two phase flow modelling and experimentation, Pisa, Italy, 22–24 September, 2004. [46] M. Ishii, N. Zuber, Drag coefficient and relative velocity in bubbly, droplet or particulate flows, AIChE J. 25 (1979) 843–855. [47] M.T. Dhotre, B. Niceno, B.L. Smith, Large eddy simulation of a bubble column using dynamic sub-grid scale model, Chem. Eng. J. 136 (2008) 337–348. [48] Y. Sato, K. Sekoguchi, Liquid velocity distribution in two-phase bubbly flow, Int. J. Multiphase Flow 2 (1975) 79–95. [49] J.C.R. Hunt, T.R. Auton, K. Sene, N.H. Thomas, R. Kowe, ICHMT International seminar on transient phenomena in multiphase flow, Dubrovnik, Yugoslavia, 1987, pp. 103–125. [50] S.S. Thakre, J.B. Joshi, CFD simulation of bubble column reactors: importance of drag force formulation, Chem. Eng. Sci. 54 (1999) 5055–5060. [51] M. Lopez de Bertodano, Turbulent Bubble flow in a triangular duct, Ph.D. Thesis, Rensselaer Polytechnic Institute, Troy New York (1991). [52] J. Smagorinsky, General circulation experiments with the primitive equations. 1. The basic experiment, Monthly Weather Rev. 91 (1963) 99– 164. [53] M. Germano, U. Piomelli, P. Moin, W.H. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A7 (1991) 1760–1765. [54] N. Zuber, J.A. Findlay, Average volumetric concentration in two phase flow systems, J. Heat Transfer 87 (1969) 453–468. [55] A.A. Kulkarni, Transport phenomena and nonlinear dynamics in multiphase systems, PhD Thesis, UICT, Mumbai (2003). [56] D. Ma, G.J. Ahmadi, A thermodynamical formulation for dispersed turbulent flows-1: basic theory, Int. J. Multiphase Flow 16 (1990) 323–340. [57] D.Z. Zhang, W.B. Vanderheyden, The effects of mesoscale structures on the disperse two-phase flows and their closures for dilute suspensions, Int. J. Multiphase Flow 28 (2002) 805–822. [58] R.F. Mudde, J.S. Groen, H.E.A. Van den Akker, Liquid velocity field in a bubble column: LDA experiments, Chem. Eng. Sci. 52 (1997) 4217–4229. 614 M.V. Tabib et al. / Chemical Engineering Journal 139 (2008) 589–614 [59] R.C. Chen, J. Reese, L.S. Fan, Flow structure in three dimensional bubble column and three phase fluidized bed, AIChE J. 40 (1994) 1093– 1104. [60] W.K. Harteveld, R.F. Mudde, H.E.A. Van Den Akker, Dynamics of a bubble column: influence of gas distribution on coherent structures, Can. J. Chem. Eng. 81 (2003) 389–394. [61] K. Tsuchiya, A. Furumoto, L.-S. Fan, J. Zhang, Suspension viscosity and bubble rise velocity in liquid–solid fluidized beds, Chem. Eng. Sci. 52 (18) (1997) 3053–3066. [62] L.A. Schiller, Z. Naumaan, A drag coefficient correlation, Ver Deutsch. Ing. 77 (1935) 138. [63] J.M. Dalla Ville, Micrometrics, Pitman Publishing Co., New York, 1948. [64] J.R. Grace, T. Wairegi, T.H. Nguyen, Shapes and velocities of single drops and bubbles moving freely through immiscible liquids, Trans. Inst. Chem. Eng. 54 (1976) 167–173. [65] S. Grevskott, B.H. Sannaes, M.P. Dudkovic, K.W. Hjarbo, H.F. Svendsen, Liquid circulation, bubble size distributions and solids movement in twoand thee-phase bubble columns, Chem. Eng. Sci. 51 (1996) 1703–1713.