Modeling of Bubble Column Reactors: Progress and Limitations

Ind. Eng. Chem. Res. 2005, 44, 5107-5151

5107

Modeling of Bubble Column Reactors: Progress and Limitations Hugo A. Jakobsen,* Ha˚ vard Lindborg, and Carlos A. Dorao Department of Chemical Engineering, Norwegian University of Science and Technology, NTNU, Sem Sælands vei 4, NO-7491 Trondheim, Norway

Bubble columns are widely used for carrying out gas-liquid and gas-liquid-solid reactions in a variety of industrial applications. The dispersion and interfacial heat- and mass-transfer fluxes, which often limit the overall chemical reaction rates, are closely related to the fluid dynamics of the system through the liquid-gas contact area and the turbulence properties of the flow. There is thus considerable interest, within both academia and industry, to improve the limited understanding of the complex multiphase flow phenomena involved, which is preventing optimal scale-up and design of these reactors. In this paper, the progress reported in the literature during the past decade regarding the use of averaged Eulerian multifluid models and computational fluid dynamics (CFD) to model vertical bubble-driven flows is reviewed. The limiting steps in the model derivation are the formulation of proper boundary conditions, closure laws determining turbulent effects, interfacial transfer fluxes, and the bubble coalescence and breakage processes. Examples of both classical and more recent modeling approaches are described, evaluated, and discussed. Physical mechanisms and numerical modes creating bubble movement in the radial direction are outlined. Special emphasis is placed on the population balance modeling of the bubble coalescence and breakage processes in two-phase bubble column reactors. The constitutive relations used to describe the bubble-bubble and bubble-turbulence interactions, the bubble coalescence and breakage criteria, and the daughter size distribution models are discussed with a focus on model limitations. The demand for amplified modeling, more accurate and stable numerical algorithms, and experimental analysis providing data for proper model validation is stressed. 1. Introduction Chemical reaction engineering (CRE) has emerged as a methodology that quantifies the interplay between transport phenomena and kinetics on a variety of scales and allows the formulation of quantitative models for various measures of reactor performance. The ability to establish such quantitative links between measures of reactor performance and input and operating variables is essential in optimizing the operating conditions in manufacturing, in determining proper reactor design and scale-up, and in correctly interpreting data in research and pilot-plant work. A starting point for reaction engineers is the formulation of a reactor model for which the basis is the microscale species mass and enthalpy balances. For practical applications, the direct solution of these equations is excessively costly, and simplifications or average representations are usually introduced. The choice of averages (e.g., global reactor volume, cross-sectional area, or length) over which the balance equations are integrated (averaged) determines the level of sophistication of the reactor model. It is very common in tubular reactors to have flow predominantly in one spatial direction, say, z. The major gradients then occur in that direction. For many cases, then, the crosssectional average values of concentration and temperature are used instead of the local values. In this way, a one-dimensional dispersion model is obtained. If the convective transport is completely dominant over the * To whom correspondence should be addressed. E-mail: [email protected]. Tel.: + 47 73594132. Fax: + 47 73594080.

diffusive transport, the diffusive term can be neglected. The resulting equations are denoted the ideal plug-flow reactor (PFR) model. When the entire reactor can be considered to be uniform in both concentration and temperature (i.e., because of very large dispersion coefficients), one can neglect gradients in all spatial directions and integrate the equations globally over all spatial dimensions (assuming convective flows at the boundaries), leading to the ideal reactor model of the continuous stirred tank reactor (CSTR). For more complex flow patterns, more elaborated and complete models are required where the flow fields are described via the solution of the Navier-Stokes equations. The understanding of the complex flow phenomena involved as well as the solution of these vector equations make the problem much more difficult to analyze within the constraint of reasonable costs and efforts. In these cases, the full set of governing equations can be averaged over local but finite spatial and temporal scales to obtain formulations that are solved with feasible space and time resolutions. The latter type of models can also be obtained using other local averaging procedures as well. Two-phase and slurry bubble columns are widely used in the chemical and biochemical industries for carrying out gas-liquid and gas-liquid-solid (catalytic) reaction processes.1-6 The historical development of bubble column modeling was discussed in a recent paper by Dudukovic.4 In the past, the axial dispersion model was commonly applied for both the gas and liquid (slurry) phases in bubble columns.2,7 The balance equations determining the liquid- (slurry-) phase axial dispersion model can be written as

10.1021/ie049447x CCC: $30.25 © 2005 American Chemical Society Published on Web 01/20/2005

5108

Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005

S d(FLvz,L ωi,L)

dz

) LFLDz,L

d2ωi,L dz2

+

/ kLa(FL,i

- FL,i) +

L

∑r γi,rRL,r

(1)

and S FLCP,Lvz,L

dTL dz

) Lλz,L

d2TL dz2

+ hWaW(TW - TL) + L

∑r (-∆Hr)Rr,L

(2)

In these balance equations, all terms should be described at the same level of accuracy. It certainly does not pay to have the finest description of one term in the balance equations if the others can only be very crudely described. Current demands for increased selectivity and volumetric productivity require more precise reactor models and also force reactor operation to churn turbulent flow, which, to a great extent, is uncharted territory. An improvement in accuracy and a more detailed description of the molecular-scale events describing the rate of generation term in the heat- and mass-balance equations has, in turn, pushed forward a need for a more detailed description of the transport terms (i.e., in the convection/advection and dispersion/ conduction terms in the basic mass and heat balances). Experimental evidence shows that the liquid axial velocity is far from being flat and independent of the radial space coordinate,2 and the use of a cross-sectional average velocity variable seems not to be sufficient. The back-mixing induced by the global liquid flow pattern has commonly been taken into account by adjusting the axial dispersion coefficient accordingly. However, even though (slurry) bubble column performance often can be fitted with an axial dispersion model, decades of research have failed to produce a predictive equation for the axial dispersion coefficient. Hence, research is in progress to quantify these parameters based on first principles (e.g., see Jakobsen et al.,8 Joshi,9 Rafique et al.,10 and references therein). However, despite their simple construction, the fluid dynamics observed in these columns is very complex. Even though CFD modeling concepts have been extended over the past two decades in accordance with the rapid progress in computer performance, the model complexity required to resolve all of the important phenomena in these systems is still not feasible within reasonable time limits. The multifluid model is found to represent a tradeoff between accuracy and computational efforts for practical applications. Unfortunately, the present models are still on a level aiming at reasonable solutions with several model parameters tuned to known flow fields. For predictive purposes, these models are hardly able to predict unknown flow fields with a reasonable degree of accuracy. It appears that CFD evaluations of bubble columns by use of multidimensional multifluid models still have very limited inherent capabilities to fully replace the empirical-based analyses (i.e., in the framework of axial dispersion models) in use today.11,12 After two decades of performing fluid dynamic modeling of bubble columns, it has been realized that there is certainly a limit to how accurate one will be able in formulating closure laws using the Eulerian framework.

A severe problem is that, although the Eulerian modeling prospects are not outstanding, no other concepts available are favorable for the purpose of predicting bubble column flows. A more relevant question is thus: What can be achieved by adopting the Eulerian multifluid modeling framework for the purpose of describing chemical processes operated in bubble column reactors? The authors intend the following summary to provide certain guidelines for the discussion but, of course, no complete answer to this question. Fluid Dynamic Modeling. Considering modeling in further detail, the general picture from the literature is that the forces acting on the dispersed phase are8 inertia, gravity, buoyancy, viscous, pressure, lift, wall, turbulent stress, turbulent dispersion, steady-drag, and added-mass forces. In the latest papers (i.e., published after the review by Jakobsen et al.8) performing 3D simulations, the force balances in vertical bubbly flows were found to be determined by only a few of these forces. The axial component of the momentum equation for the gas phase is dominated by the pressure and steady-drag forces only, indicating that algebraic slip models might be sufficient,13-15 whereas most multifluid models also retain the inertia and gravity (buoyancy) terms. The axial momentum balance for the liquid phase considers the inertia, turbulent stress, pressure, steadydrag, and gravity forces. Only pressure and buoyancy forces are acting on a motionless bubble in a liquid at rest. In the radial and azimuthal directions, the force balances generally includes the steady-drag force, i.e., a force that opposes motion, whereas the pertinent forces causing motion are more difficult to define. In a very short inlet zone, the wall friction is likely to induce a radial pressure gradient that pushes the gas bubbles away from the wall, whereas a few column diameters above the inlet, the radial pressure gradient vanishes. It is still an open question whether this pressure gradient is sufficient to determine the phase distributions observed in these systems. It is expected that the presence of the wall induce forces that act on the dispersed particles farther from the inlet, but there is no general acceptance on the physical mechanisms and formulations of these forces. It is also a matter of discussion whether this wall effect should be taken into account indirectly through the liquid wall friction or, in view of the model averaging performed, directly as a force in the gas-phase equations. The bulk lift forces do not induce any radial bubble movement without an initial velocity gradient, so other forces are important in the initial phase of the flow pattern development. In the generalized drag formulation, the interfacial coupling is expressed as a linear sum of independent forces; this point of view is probably not valid when the void fraction exceeds a few percent. Moreover, the parametrizations used for the coefficients occurring in the interfacial closures vary significantly, especially at higher void fractions. Single-particle drag, added-mass, and lift coefficients are most frequently used, whereas swarm corrections have been included in some codes. For higher void fractions, other rather empirical corrections have been introduced as well.16 The flow of the continuous phase is considered to be initiated by a balance between the interfacial particlefluid coupling and wall friction forces, whereas the fluidphase turbulence damps the macroscale dynamics of the bubble swarms, thereby smoothing the velocity and volume fraction fields. Temporal instabilities induced

Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5109

Figure 1. Classification of regions accounting for the macroscopic flow structures: left, 2D bubble column;21 right, 3D bubble column.23

by the fluid inertia terms create inhomogeneities in the force balances. Unfortunately, proper modeling of turbulence is still one of the main open issues in gas-liquid bubbly flows, and this flow property can significantly affect both the stresses and the bubble dispersion.15 It was shown by Svendsen et al.,17 among others, that the time-averaged experimental data on the flow pattern in cylindrical bubble columns apparently becomes close to axisymmetric. Fair agreement between experimental and simulated results are generally obtained for the steady velocity fields in both phases, whereas the steady phase distribution is still a problem. Therefore, it was anticipated that 2D axisymmetric simulations would capture the pertinent time-averaged flow pattern needed for the analysis of many (not all) mechanisms of interest for chemical engineers. Sanyal et al.18 and Krishna and van Baten,19 for example, stated that 2D models provide good engineering descriptions, although they are not able to capture the high-frequency unsteady behavior of the flow, and can be used for approximately predicting the low-frequency time-averaged flow and void patterns in bubble columns. The early 2D steady-state model proposed by Jakobsen20 was able to capture the global flow pattern fairly well, provided that a large negative lift force was included. However, after the first elaborated experimental studies on 2D rectangular bubble columns were published by Tzeng et al.21 and Lin et al.,22 it was commonly accepted that time-average computations cannot provide a rational explanation of the transport processes of mass, momentum, and energy between the bubbles and liquid. The experimental data obtained were analyzed, and sketches of their interpretations of the dynamic flow patterns in both 2D and 3D columns were given, as shown in Figure 1. It was concluded that a proper bubble column model should consider the transient or instantaneous flow behavior. A few years later, Sokolichin and Eigenberger,24 Sokolichin et al.,25 and Sokolichin and Eigenberger26 claimed that dynamic 3D models were needed to provide sufficient representations of the high-frequency unsteady behavior of these flows. Very different dynamic flow patterns can result in quantitatively similar long-time-averaged flow profiles. This limits the use of long-time-averaged flow profiles for validation of bubbly flow models. van den

Figure 2. Lateral movement of the bubble hose in a flat bubble column.28 Reprinted with permission from Elsevier.

Akker,12 among others, questioned the early turbulence modeling performed (or rather the lack of any) in these studies and argued that the apparently realistic simulations of the transient flow characteristics could be numerical modes rather than physical ones (see also Sokolichin et al.15). Insufficiently fine grids might have been used in the simulations, resulting in numerical instabilities that could be erroneously interpreted as physical ones. However, after the observation of Sokolichin and Eigenberger was reported, intensive focus on these instability issues were made, mainly by researchers from the fluid dynamics community. The experimental data provided by Becker et al.27,28 (see Figure 2) still serve as a benchmark test and are often used for validation of dynamic flow models. The numerical investigations were restricted to bubbly flow hydrodynamics (i.e., no reactive systems were analyzed), where additional simplifications were made: isothermal conditions, no interfacial mass transfer, constant liquid density, gas density constant or depending on local pressure as described by the ideal gas law, and no bubble coalescence and breakage. However, after the dynamic flow structures were observed, bubble columns have generally been simu-

5110


lated using either 2D or 3D dynamic models for both cylindrical18,29-36 and rectangular 2D and 3D column geometries.14,15,34,35,37-43 The gas is introduced adopting both uniform and localized feedings at the bottom of the column. The modeling of systems uniformly gassed at the bottom is more difficult than the modeling of partly aerated ones. Simulating systems with continuous liquid flow is also more difficult than keeping the liquid in batch mode. Finally, it is also noted that the different research groups applied both commercial codes (e.g., CFX, FLUENT, PHOENICS, CFDLIB, ASTRID, NPHASE) and several in-house codes in which the inherent choices of numerical methods, discretizations, grid arrangements, and boundary implementations vary quite widely. These numerical differences alter the solutions to some extent, so it should not be expected that the corresponding simulations will provide identical results. An open and unified research code available for all research groups could assist by eliminating any misinterpretations of numerical modes as physical mechanisms and visa versa. Considering the interfacial and turbulent closures for vertical bubble-driven flows, no extensive progress has been observed in the later publications. However, two diverging modeling trends seem to emerge because of the lack of understanding of the phenomena involved and how to deal with these phenomena within an average modeling framework. One group of papers considers only phenomena that can be validated with the existing experimental techniques and thus contains a minimum number of terms and effects. Papers in the other group include a large number of weakly founded theoretical hypothesis and relations intended to resolve the missing mechanisms. Steady-State or Dynamic Simulations, Closures, and Numerical Grid Arrangements. Not only dynamic models have been adopted investigating these phenomena. Lopez de Bertodano,44 for example, used a 3D steady finite-volume method (FVM) (PHOENICS) with a staggered grid arrangement to simulate turbulent bubbly two-phase flow in a triangular duct. In this study, the lift, turbulent dispersion, and steady-drag forces were assumed to be dominant. Standard literature expressions were adopted for the drag and lift forces, whereas a crude model for the turbulent dispersion force was developed. An extended k- model, considering bubble-induced turbulence, was also developed for the liquid-phase turbulence. It was proposed that the shear-induced turbulence and the bubbleinduced turbulence could be superimposed. The lift force was found to be essential to reproducing the experimentally observed wall void peaks satisfactorily. Anglart et al.45 adopted many of the same closures within a 2D steady version of the same code (PHOENICS), predicting low-void bubbly flow between two parallel plates. They found satisfactory agreement with experimental data when applying drag, added-mass, lift, wall (lubrication), and turbulent diffusion forces in their study. The extended k- model of Lopez de Bertodano44 was applied for the liquid-phase turbulence. In a more recent paper, Antal et al.46 adopted very similar closures for 3D steady-state bubble column simulations using the NPHASE code. The NPHASE CMFD code employs an FVM on a collocated grid. A three-field multifluid model formulation was used to simulate two-phase flow in a bubble column operating in the churn-turbulent flow regime. The gas phase was subdivided into two fields

(i.e., small and large bubbles) to more accurately describe the interfacial momentum-transfer fluxes. The third field was used for the liquid phase. The model results were validated against a few time-averaged data sets for the liquid axial velocity and the gas volume fraction. Global flow patterns for all three fields and the overall gas volume fractions were shown. The simulations were in fair agreement with the experimental observations. Using a 2D in-house FVM code with a staggered grid arrangement, Dhotre and Joshi47 predicted the flow pattern, pressure drop, and heat-transfer coefficient in bubble column reactors. The model used contained steady-drag, added-mass, and lift forces, as well as a reduced pressure gradient formulated as an apparent form drag. Turbulent dispersion was taken into account by use of mass-diffusion terms in the continuity equations. A low-Reynolds-number k- model was incorporated, merely constituting a standard k- model with modified treatment of the near-wall region. The turbulence model used contained an additional production term accounting for the large-scale turbulence produced within the liquid flow field because of the movement of the bubbles. A semiempirical mechanical energy balance for the gas-liquid system was imposed. The simulated results were in surprisingly good agreement with experimental literature data on the axial liquid velocity, gas volume fraction, friction multiplier, and heat-transfer coefficient. Deen et al.40 used the lift force in addition to the steady-drag and added-mass forces in their dynamic 3D model to obtain the transverse spreading of the bubble plume that is observed in experiments. A prescribed zero-void wall boundary was used to force the gas to migrate away from the wall. The continuous-phase turbulence was incorporated in two different ways, using either a standard k- or a VLES model. The effective viscosity of the liquid phase was composed of three contributions: the molecular, shear-induced turbulent, and bubble-induced turbulent viscosities. The calculation of the turbulent gas viscosity was based on the turbulent liquid viscosity as proposed by Jakobsen et al.8 These simulations were performed using the commercial code CFX, so an FVM on a collocated grid was employed. Sample results simulating a square 3D column at low void fractions using the 3D VLES model of Deen et al.40 are shown in Figure 3. Krishna and van Baten29,35 used a steady-drag force as the only interfacial momentum coupling in their transient 2D and 3D three-fluid models with an inherent prescribed and fixed bimodal distribution of the gas bubble sizes. The two bubble classes were denoted small and large bubbles. The small bubbles were in the range of 1-6 mm, whereas the large bubbles were typically in the range of 20-80 mm. An unfortunate simplification made in the model is that there are no interactions or exchanges between the small and large bubble populations, as will be discussed later. An FVM on an unstaggered grid was used to discretize the equations (i.e., in CFX). No-slip conditions at the wall were used for both phases. A k- model was applied for the liquid-phase turbulence, whereas no turbulence model was used for the dispersed phases. To prevent a circulation pattern in which the liquid flowed up near the wall and came down in the core, the large bubble gas was injected on the inner 75% of the radius. The time-averaged volume fraction and velocity profiles calculated from the predicted 3D flow field were in reasonable agreement with experimental


Figure 3. Snapshots of the instantaneous isosurfaces of Rd ) 0.04 and liquid velocity fields after 30 and 35 s for the LES model.40 Reprinted with permission from Elsevier.

data. Lehr et al.32 used a similar three-fluid model (implemented in CFX), combined with a simplified population balance model for the bubble size distribution. The simplified population balance relation used contained semiempirical parametrizations for the bubble coalescence and breakage phenomena. It was concluded that the calculated long-time-average volume fractions, velocities, and interfacial area density were in good agreement with the experimental data. Pfleger et al.39 applied a two-fluid model using the same code (CFX) to a 2D rectangular column with localized spargers. It was concluded that a 3D model including the steady-drag force and a standard k- model is sufficient to correctly capture the unsteady behavior of bubbly flow with very low gas void fractions. The dispersed phase was considered laminar. Sanyal et al.18 and Bertola et al.34 used similar two-fluid models and FVMs on collocated grid arrangements (FLUENT and CFDLIB) to simulate both cylindrical and rectangular columns, confirming the results found by Pfleger et al.39 It should be noted, however, that Bertola et al.34 solved the k- turbulence model for both phases, whereas Sanyal et al.18 adopted an approach based on Tchen’s theory48-51 to predict turbulence in the dispersed phase. Mudde and Simonin38 performed both 2D and 3D simulations of a 2D rectangular column using a similar FVM on a collocated grid arrangement (ASTRID). Their two-fluid model contained an extended k- turbulence model formulation for the liquid-phase turbulence, as well as drag and added-mass forces. The dispersed-phase turbulence was assumed to be steady and homogeneous and was described by the extended Tchen’s theory approach mentioned above. This code predicted reasonable highfrequency oscillating flows only when the added-mass force was included; without this force low-frequency,

almost steady flows were obtained. Using an FD scheme, Tomiyama et al.52 reported that the transient transverse migration of bubble plumes in vessels can be well predicted including steady-drag, added-mass, and lift forces in their two-fluid model describing laminar flows containing very low void fractions (i.e., below 0.5%). Later, Tomiyama53 and Tomiyama and Shimada54 used the bubble-induced turbulence model of Sato and Sekoguchi55 and Sato et al.56 and an extended k- model in their work on turbulent flows. In a recent paper, Sokolichin et al.15 concluded that the model by Sato and Sekoguchi55 for bubble-induced turbulence strongly underestimates the turbulence level in a number of test cases. Another disadvantage of this approach consists of their local properties, because it considers the increase of the turbulence intensity only locally in the reactor where the gas phase is actually present. In reality, the turbulence induced by the bubbles at some given point can spread and affect regions farther from the turbulence source. Oey et al.14 applied a 3D twofluid model containing the steady-drag, added-mass, and turbulent dispersion forces together with an extra source term in the k- turbulence model to account for the effects of the interface in their in-house staggered FVM code (ESTEEM). Because of the controversy in the literature regarding dispersed-phase turbulence, both laminar and turbulent gas simulations were performed. In the turbulence case, the extended Tchen’s theory approach was adopted. The liquid tangential velocity components close to the wall were found using wall functions, whereas no wall friction was taken into account for the dispersed phase. They found that, in 3D, the steady-drag force was sufficient to capture the global dynamics of the bubble plume whereas the other forces moreover had secondary effects only. Furthermore, Oey et al.14 also concluded that the discretization schemes adopted for the convection terms in the fluid momentum and dispersed-phase continuity equations had a severe impact on the solutions. This is an important aspect of any multiphase flow modeling: the numerical and modeling issues cannot be investigated (completely) separately, as their interplay is of considerable importance.8,11,57 Most of the studies mentioned above adopted a kind of k- model to describe the liquid turbulence in the system, whereas there is less consensus regarding whether the dispersed phase should be considered turbulent or laminar, or even whether any deviating stress terms should remain in the dispersed-phase equations at all. However, even the k- model predictions are questioned by Deen et al.40 and Bove et al.43 This group showed that only low-frequency unsteady flow is obtained using the k- model because of overestimation of the turbulent viscosity. These model predictions were found not to be in satisfactory agreement with the more high-frequency experimental results. On the other hand, when a 3D Smagorinsky LES (large-eddy simulation) model was used instead, the strong transient movements of the bubble plume observed in experiments were captured. For completeness, it should also be mentioned that, in a few recent papers,58,59 it has been claimed that 2D mixture model formulations can be used to reproduce the time-dependent flow behavior of 2D bubble columns. It has been found that the crucial physics can be captured by 2D models if suitable turbulence models are used. The predictions reported by Bech58 rely on the

5112


inclusion of a mass-diffusion term in the dispersedphase continuity equation, whereas the predictions of Cartland Glover and Generalis59 rely on the inclusion of a Reynolds stress model. However, the great majority of the investigations reported conclude that, for both rectangular and cylindrical columns, the high-frequency instabilities are 3D and have to be resolved through the use of 3D models. An apparent conclusion drawn from these papers (although not explicitly stated), that the fluid dynamics had to be essentially perfect before the chemistryrelated topics could be considered, might have been a severe hindrance in the development of integrated models of interest for chemical engineers. The main limitation is said to be the tremendous CPU demands needed for the high-frequency instabilities (assuming that these are important for the chemical conversion), making such 3D simulations infeasible for most research groups. Unfortunately, there are several indications that these local high-frequency dynamics and coherent structures of the flow are important in determining the conversion of the system.15,37,42,60,61 Mixing times predicted by steady flow codes, for example, are found not to be in accordance with experimental data, whereas mixing times predicted by high-frequency transient flow codes are in fair agreement with the corresponding measurements. The interfacial and turbulence closures suggested in the literature also differ in considering the anticipated importance of the bubble size distributions. It thus seemed obvious for many researchers that further progress on the flow pattern description was difficult to obtain without a proper description of the interfacial coupling terms, and especially of the contact area or projected area for the drag forces. The bubble column research thus turned toward the development of a dynamic multifluid model that is extended with a population balance module for the bubble size distribution. However, the existing models are still restricted in some way or another because of the large CPU demands required by 3D multifluid simulations. Multifluid Models and Bubble Size Distributions. To gain insight into the capability of the present models to capture physical responses to changes in the bubble size distributions, a few preliminary analyses have been performed adopting the multifluid modeling framework. Carrica et al.13 developed the most comprehensive multifluid/population balance model known to date (to the knowledge of the authors) for the description of bubbly two-phase flow around a surface ship. In the multifluid model part, the inertia and shear stress tensors were assumed to be negligible for the gas bubble phases or groups. The interfacial momentum transfer terms included for the different bubble groups are steady-drag, added-mass, lift, and turbulent dispersion forces. Algebraic turbulence models were used both for the liquid-phase contributions and for the bubbleinduced turbulence. The two-fluid model was solved using an FVM on a staggered grid. In the population balance part of the model the intergroup transfer mechanisms included were bubble breakage and coalescence and the dissolution of air into the ocean. Fifteen size groups were used with bubble radii at normal pressure between 10 and 1000 µm. It was found that intergroup transfer is very important in these flows for determining both a reasonable two-phase flow field and

the bubble size distribution. The population balance was discretized using a multigroup approach. It was pointed out that the lack of validated kernels for bubble coalescence and breakage limited the accuracy of the model predictions. Politano et al.62 adapted the population balance model of Carrica et al.13 for the purpose of 3D steady-state simulation of bubble column flows. However, no details regarding the necessary model modifications were provided. In a later study by Politano et al.,63 a 2D steady-state version of this model was applied for the simulation of polydisperse two-phase flow in vertical channels. The two-fluid model was modified using an extended k- model for the description of liquid-phase turbulence. A two-phase logarithmic wall law was developed to improve on the boundary treatment of the k- model. The interfacial momentum-transfer terms included for the different bubble groups were the steadydrag, added-mass, lift, turbulent dispersion, and wall forces. The two-fluid model equations were discretized using an FD method. The bubble mass was discretized in three groups. The effect of bubble size on the radial phase distribution in vertical upward channels was investigated. Comparing the model predictions with experimental data, it was concluded that the model is able to predict the transition from the near-wall gas volume fraction peaking to the core peaking beyond a critical bubble size. Tomiyama53 and Tomiyama and Shimada54 adopted an (N + 1)-fluid model for the prediction of 3D unsteady turbulent bubbly flows with nonuniform bubble sizes. Among the (N + 1) fluids, one fluid corresponds to the liquid phase, and the N fluids correspond to gas bubbles. To demonstrate the potential of the proposed method, unsteady bubble plumes in a water-filled vessel were simulated using both (3 + 1)-fluid and two-fluid models. The gas bubbles were classified and fixed in three groups only, so a (3 + 1)- or four-fluid model was used. The dispersions investigated were very dilute so the bubble coalescence and breakage phenomena were neglected, whereas the inertia terms were retained in the three bubble-phase momentum equations. No population balance model was then needed, and the phase continuity equations were solved for all phases. It was confirmed that the (3 + 1)-fluid model gave better predictions than the two-fluid model for bubble plumes with nonuniform bubble sizes. As mentioned earlier, three-fluid models have also been used by a few groups.32,35 Krishna and van Baten35 solved the Eulerian volume-averaged mass and momentum equations for all three phases. However, no interchange between the small- and large-bubble phases was included so each of the dispersed bubble phases exchanges momentum only with the liquid phase. No population balance model was used as the bubble coalescence and breakage phenomena were neglected. Lehr et al.32 extended the basic three-fluid model by including a simplified population balance model for the bubble size distribution. In most studies reported so far, two-fluid models have been used,13,31,33,41,62,64,65 assuming that all of the particles have the same average velocities. In other words, possible particle collisions due to buoyancy effects are neglected even though these contributions have not been proven insignificant. This means that the two momentum equations for the two phases are solved together with the continuity equation of the liquid phase and the N population balance equations for the dispersed


Figure 4. Simulated steady-state results obtained with (s) and without (- - -) the population balance at the axial level 2 m above the column inlet. Points (‚) are experimental data.76 vsd ) 8 cm/s.

phases.13 An alternative, often used when commercial CFD codes are used because of limited access to the solver routines, is to solve the full two-fluid model in the common way using the dispersed-phase continuity equation together with the two momentum equations and the liquid-phase continuity. Within the IPSA-like calculation loop, the N - 1 population balance equations are solved in another step considering additional transport equations for scalar variables. Unfortunately, using this approach, it can be difficult to ensure mass conservation for the dispersed phase, and/or negative concentrations can be predicted for the last class. From the solution of the size distribution of the dispersed phase, the Sauter mean diameter is calculated. This diameter is then used to compute the contact area, so that the two-way interaction between the flow and the bubble size distribution is established. When dilute dispersions are considered, the interfacial momentumtransfer fluxes due to particle collisions, coalescence, and breakage phenomena are normally neglected. Politano et al.62 are probably the only ones to have considered the interfacial momentum-transfer fluxes between the bubble groups induced by bubble coalescence and breakage, but still, they did not include the bubble collisions resulting in rebound. Several recent studies indicate that algebraic slip models are sufficient for modeling the flow pattern in bubble columns,13-15 as only the pressure and steadydrag forces dominate the axial component of the gas momentum balance. Therefore, the population balance model can be merged with an algebraic slip model to reduce the computational cost required for preliminary analysis.66 Furthermore, by adopting this concept, the restriction used in the two-fluid model of assuming that all of the particles have the same velocity can be avoided. This means that only one set of momentum and continuity equations for the mixture is solved, together with the N population balance equations for the dis-

persed phases. The individual phase and bubble-class velocities are calculated from the mixture velocities using algebraic relations. Even simpler approaches are used in solving a single transport equation for one moment of the population balance only, determining a locally varying mean particle size or the interfacial area density.36,67-69 Luo70 initiated the population balance investigations within our group. The main contribution was a closure model for binary breakage of fluid particles in fully developed turbulence flows based on isotropic turbulence and probability theories.71 The authors also claimed that this model contains no adjustable parameters, although a better phrase might be “no additional adjustable parameters”, as both the isotropic turbulence and probability theories involved contain adjustable parameters and distribution functions. Hagesaether et al.72-75 continued the population balance model development within the framework of an idealized plug-flow model, whereas Bertola et al.66 combined the extended population balance module with a 2D algebraic slip mixture model for the flow pattern. Bertola et al.66 studied the effect of the bubble size distribution on the flow fields in bubble columns. They used an extended k- model to describe the turbulence of the mixture flow. Two sets of simulations were performed, i.e., both with and without the population balance. Four different superficial gas velocities, i.e., 2, 4, 6, and 8 cm/s were used, and the superficial liquid velocity was set to 1 cm/s in all cases. The population balance contained six prescribed bubble classes with diameters set to d1 ) 0.0038 m, d2 ) 0.0048 m, d3 ) 0.0060 m, d4 ) 0.0076 m, d5 ) 0.0095 m, and d6 ) 0.0120 m. Figures 4 and 5 show simulated and experimental results for the superficial gas velocity vsg ) 8 cm/s. Figure 4 shows the axial and radial liquid velocity

5114


Figure 5. Calculated steady-state bubble number density at the axial level 2 m above the column inlet. Points (‚) are experimental data.76 vsd ) 8 cm/s.

components, the axial gas velocity component, and the gas fraction 2.0 m above the column inlet. Figure 5 shows the number density in each class 2.0 m above the inlet. The results from the two simulations (i.e., with and without the population balance) are nearly identical. In both cases, the simulated results are in fair agreement with the experimental data, but in the center of the reactor, the deviations between the simulated and experimental velocity and void fraction profiles are rather large. The number density 2.0 m above the inlet is shown in Figure 5. In general, it was found that the initial bubble size was not stable and was further determined by break-up and coalescence mechanisms. The simulation provides results in fair agreement with the experimental data for classes 3-6, for which the bubble number densities are of the same order of magnitude as the experimental data. In bubble classes 1 and 2, the experimental bubble number densities are considerably underestimated in the simulations. However, in other cases, the model predictions deviate much more from each other and were in poor agreement the experimental data on measurable quantities such as phase velocities, gas volume fractions, and bubble size distributions. An obvious reason for this discrepancy is that the breakage and coalescence kernels rely on ad hoc empiricism determining the particle-particle and particle-turbulence interaction phenomena. The existing parametrizations developed for turbulent flows are high-order functions of the local turbulent energy dissipation rate that is often determined by the k- turbulence model. This approach is not accurate enough,32,33,77 as this model variable (i.e., ) merely represents a closure for the turbulence integral length scale with model parameters fitted to experimental data for idealized single-phase flows. The population balance kernels are also difficult to validate on the mesoscale level, as the physical mechanisms involved (e.g., con-

sidering eddies and eddy-particle interactions) are vague and not clearly defined. If possible, the coalescence and breakage closures should be reparametrized in terms of measurable quantities. Well-planed experimental analysis of the mesoscale phenomena are then required to provide data for proper model validation. Initial and Boundary Conditions. Initial and boundary conditions are also very important parts of any model formulation. However, there is still very limited knowledge regarding the formulation of proper boundary conditions even for simple bubble columns.11,12 Inlet boundary conditions for the local gas volume fraction, bubble size distribution, and local gas velocity components are difficult to determine, although a number of more or less well founded suggestions have been reported in the literature.43,66,78 The complex flows in this section have been studied experimentally in several papers,61,79,80 and there is clearly a need for improvements as Harteveld et al.61 found significant discrepancies between their experimental observations and the present model predictions regarding the onset of the significant flow instabilities. Harteveld et al.61 claimed that uniform aeration gives only a very small entrance region and no large-scale circulation or coherent structures. This visual observation clearly contradicts most of the model predictions cited above. The specification of proper outlet boundary conditions is also a problem for these flows, as the recirculating motion of the liquid phase continues as long as gas bubbles are present at sufficiently high fluxes. Consistent formulations of such boundaries have not yet been reported.81 This limitation restricts the use of explicit discretization schemes to situations with low gas fractions or cases where the local recirculating flux contains very little gas. For implicit solution methods, approximate and rather crude outlet conditions are used for the flow variables. The boundary is usually idealized considering the liquid


phase in batch mode only.78 Following the experimental laboratory practice of keeping the height of the gasliquid dispersion lower than the actual column height, the top surface of the column can be modeled as an outlet for both the gas and liquid phases. It is anticipated that the solution of the model equations will determine the actual height of the gas-liquid dispersion and only gas will exit from the column outlet. In adopting this approach, one has to ensure that the modeling closures are well posed and that the discretization scheme applied to the governing equations as well as the iterative solver used is capable of handling the steep gradients in the volume fraction and density profiles (discontinuities) that occur when the continuous phase change from being liquid below the gas-liquid interface to gas above it. Because of the large density difference between these phases, such an attempt very often leads to nonphysical pressure and velocity values close to the interface and encounters severe convergence difficulties. To enhance the convergence behavior in the interface region, empirically adjusted smooth profiles of the continuous-phase density profile and/or the void fraction profile might help to maintain a fairly stable solution.82 For incompressible flows, this numerical approach could enforce mass conservation, whereas for reactive and other density-varying systems, mass conservation can be a severe problem so this boundary treatment should be avoided. As numerical instabilities and inaccurate mass conservation are frequently encountered in solving problems containing sharp interfaces within the calculation domain, the solution domain is often restricted to the height of the gas-liquid dispersion. In this case, the local liquid velocity components normal to the outlet plane are fixed at zero, as there is no net flux through the column outlet cross section. The top surface of the solution domain is assumed to coincide with the free surface of the dispersion, which might or might not be assumed flat. This assumption is rather crude, as an exact value of the gas-liquid dispersion height is not known a priori. To induce apparently physical flow characteristics at the outlet, approximate boundaries are required for the other flow variables. The tangential shear stress and the normal fluxes of all scalar variables are set to zero at the free surface. The gas bubbles are free to escape from the top surface. In commercial codes, this implementation might not be possible, and further approximations have been proposed. In some cases, the top surface of the dispersion is defined as an “inlet” where the normal liquid velocity is set to zero and the normal gas velocity is set to an approximate terminal rise velocity. In other cases, the top surface of the dispersion is modeled as a “no shear wall”. This boundary sets both the gas and liquid velocities to zero. To represent escaping gas bubbles, an artificial gas sink is defined for all of the grid cells attached to the top surface. A similar approach was used by Lehr et al.32 The free surface assumed to be located at the top of the column was replaced by an apparent semipermeable wall. In this way, the gas could leave the system, whereas the liquid surface acted as a frictionless wall for the liquid. The liquid was considered to leave the system through an outflow at the periphery of the column. Bove et al.43 specified an outlet pressure boundary instead and determined the axial liquid velocity components in accordance with a global mass balance. This

approach is strictly valid only when the variations in liquid density due to interfacial mass transfer or temperature changes are negligible, as the local changes are not known a priori. Furthermore, in many industrial systems, the liquid phase is not operated in batch mode; rather, a continuous flow of the liquid phase has to be allowed. However, because of numerical problems, most reports on bubble column modeling introduce the simplifying assumption that the continuous phase is operated in batch mode. Further work is needed on the continuous-mode boundary conditions. The assumption of cylindrical axisymmetry in the computations prevents lateral motion of the dispersed gas phase and leads to an unrealistic radial phase distribution that has also been reported by other authors.35 Krishna and van Baten35 obtained better agreement with experiments when a 3D model was applied. However, experience shows that it is very difficult to obtain reasonable time-averaged radial void profiles even in 3D simulations, even though the predicted radial velocity profiles seem reasonable. At the wall boundary, both the 2D and 3D turbulentviscosity-based models rely on the assumption that the single-phase logarithmic law of the wall is still valid for bubble-driven upward flows in bubble columns. Troshko and Hassan83,84 claimed that this assumption is reasonable for dilute and downward bubbly flows, but not recommended for upward bubble-driven flows as found in bubble columns. A corresponding two-phase logarithmic wall law was derived with the intention that it would be in better agreement with experimental data on homogeneous bubbly flows. Numerical Schemes and Algorithms. In terms of progress on the numerical side, several novel schemes and algorithms for solving the fluid dynamic part of the model have been published (the solution methods intended for the population balance equation will be discussed later). This work has concentrated on several items. Most important, as mentioned earlier, one avoids using the very diffusive first-order upwind schemes while discretizing the convective terms in the multifluid transport equations. Instead, higher-order schemes, e.g., second- or third-order TVD (total variation diminishing) schemes, which are much more accurate, have been implemented in the codes.8,14,25,26,57 The numerical truncation errors induced by the discretization scheme adopted for the convective terms can severely alter the numerical solution, and this can destroy the physics reflected by the model equations. It is also well-known that the strong coupling between the phasic equations prevents efficient and robust convergence for the implicit iteration process.85 Exchanging the well-known partial elimination algorithm (PEA) of Spalding86,87 and reducing the interaction between the phasic velocities in the drag terms of the momentum equations with a coupled solver46,88-93 that simultaneously iterates on one velocity component of all phases seems to improve the numerical stability and the overall convergence rate.94 In addition, the discretization schemes together with the solution algorithms lead to large sparse linear systems of algebraic equations that need to be solved. Previously, the TDMA algorithm was applied to multidimensional problems to determine a linewise Gauss-Seidel approach. During the past decade, there have been developments in full-field solvers along the lines of Krylov subspace methods and in the field of multigrid

5116


schemes.90,91,94-96 These methods make effective use of sparsity and are efficient methods for the solution of large linear systems.96 Furthermore, the original interphase-slip algorithm (IPSA) of Spalding86 was developed to introduce an implicit coupling between pressure and volume fractions. The algorithm contains an attempt to approximate the simultaneous change of volume fractions and velocity with pressure. However, most versions of the IPSA algorithm were merely extensions of the single-phase SIMPLE approach; thus, the pressure was computed by assuming that all of the velocity components, but not the volume fraction variable, depend on pressure changes. The pressure-volume fraction relationship was not considered in a satisfactory implicit manner. Therefore, the extension of the singlephase pressure-velocity coupling to multifluid models led to low convergence rates and pure robustness of the iterating procedure.20 Lately, numerical methods originally intended for multiphase models, rather than being extended single-phase approaches, have been investigated.46,89,91-93 A fully implicit coupling of the phasic continuity and compatibility equations within the framework of pressure-volume fraction-velocity correction schemes seems to have potential if the resulting set of algebraic equations can be solved by an efficient and robust parallel solver.93 However, so far, severe stability problems have been identified within the iterative solution process. The numerical properties of the resulting set of algebraic equations are not optimized for robust solutions. In summary, significant improvements of the numerics have been obtained during the past decade, but the present algorithms are still far from being sufficiently robust and efficient. Further work on the numerical solution methods in the framework of FVMs should proceed along the paths sketched above. Very few attempts exploring the capabilities of alternative methods such as FEMs97 and fully spectral methods have been reported. Chemical Reaction Engineering. As a consequence of the large number of modeling limitations discussed above, caution concerning the predictive power of the multifluid model applied to bubbly flows is certainly justified. However, even though there are obviously many open questions and shortcomings related to the fluid dynamic modeling of bubble columns, preliminary attempts have been performed in predicting chemical reactive systems. For example, the early CFD models that were developed in our group to describe global steady flow pattern were tested aiming at predicting chemical conversion of a reactive system operated in a bubble column.20,98 The system investigated was CO2 absorption in a methyldiethanolamine (MDEA) solution. The starting point for the numerical investigation was a steady two-fluid flow model tuned to the air/ water system. This air/water model was then applied to the reactive system without retuning of any model parameters but updating the physical properties in accordance with the reactive system. It was found that, with this model, the global flow pattern, the interfacial mass-transfer fluxes, and the conversion were all still difficult to predict because of the limited accuracy reflected by the interfacial coupling models and especially the relations used for the contact area (and the projected area). Several semiempirical models for the locally varying mean bubble size and contact areas have been suggested20,99 with limited success. The accuracy

of the experimental data used for model validation was also questioned. According to Dudukovic et al.,5 still no fundamental models for the interfacial heat- and mass-transfer fluxes have been coupled successfully to the flow models, and reliable reactor performance predictions based on these models are not imminent. The mechanisms of coalescence and breakage are far from being sufficiently understood yet. The physicochemical hydrodynamics determining the bubble coalescence and breakage phenomena cannot be captured by any continuum model formulation. Therefore, the fluid dynamic models described above do not significantly improve on the prediction of the interfacial transfer fluxes and thus the chemical conversion, as the pertinent physics on bubble, interface, and molecular scales still have to be considered using empirical parametrizations. A multiscale reactor model with an inherent two-way coupling between mechanisms on scales ranging from discrete molecular ones (molecular dynamics) to continuum macroscale reactor scales should be elucidated.100 The weakest link in modeling reactive systems within bubble columns is thus still the fluid dynamic part, considering multiphase turbulence modeling, interfacial closures, and especially the impact and descriptions of bubble size and shape distributions. For reactive systems, the estimates of the contact areas and thus the interfacial mass-transfer rates are likely to contain large uncertainties. The preliminary simulations performed to date clearly indicate that future work should also consider the possibility of developing more efficient and stable numerical models for the integrated multiphase/ population balance models. Bertola et al.66 found that the CPU demand was increased by about 15% for every bubble class added to a simulation. Outline of the Paper. In the following sections, a fairly rigorous modeling framework is outlined with a focus on population balance modeling. Emphasis is placed on the closures for the coalescence and breakage terms. Bottlenecks that need further considerations are identified. The numerical methods considered good candidates in solving the population balance model for bubbly flows are the moment, volume (i.e., the class and multigroup), and spectral methods. The optimal choice of numerical solution methods is discussed. Finally, concluding remarks are given. 2. Modeling Framework Average multifluid models with an integrated population balance module have been found to represent a tradeoff between accuracy and computational efforts for practical applications. If bubbles of different mass have to be considered, separate continuity and momentum balances are required for each bubble size and for the continuous liquid phase in a rigorous model formulation. The multifluid model equations are listed in the following sections, together with the interfacial closures that are adopted in most gas-liquid (two-fluid) analyses. Multifluid Formulation. The multifluid model represents a direct extension of the well-known two-fluid model and is described in detail in Carrica et al.,13 Pfleger et al.,30 Pfleger et al.,39 Tomiyama and Shimada,54 and Fan et al.101 The governing set of equations consists of the continuity and momentum equations for


N + 1 phases; one phase corresponds to the liquid phase, and the remaining N phases are gas bubble phases. The fundamental form of the multifluid continuity equation for phase k reads

∂ ∂t

N

(RkFk) + ∇‚(RkFkvk) )

Γk,l ∑ l)1

The drag coefficient can be estimated using the relation suggested by Tomiyama et al.102

{ [

CD,d ) max min

(3)

where Fk is the density and Rk is the volume fraction of phase k. The first term on the left-hand side denotes the transient change of mass within a control volume, the second term denotes the convective flux of mass through the surfaces of the control volume, and the term on the right-hand side describes the net mass-transfer flux to phase k from all other phases l. The phasic volume fractions also satisfy the compatibility condition

0.15ReB,d0.687),

E0,d )

∂

j k) + (RkFkvk) + ∇‚(RkFkvkvk) ) -Rk∇p - ∇‚(Rkσ

gz|Fc - Fd|ds,d2 σI

ReB,d )

Mk,l ∑ l)1

(5)

The terms on the left-hand side of eq 5 denote the inertia force, whereas the terms on the right-hand side denote all of the additional forces acting on phase k. These are the pressure force, the deviating normal and shear stresses, the gravitational force, and the interfacial momentum-transfer terms accounting for all momentumtransfer fluxes between phase k and the other N phases. When sufficiently dilute dispersions are considered, only particle-fluid interactions are significant, and the two-fluid closures can be adopted. For the gas bubble phases (d), the interaction with the continuous liquid phase (c) and the wall (w) through the last term on the right-hand side of eq 5 is expressed as N+1

Md,l ≈ Md,c + Md,w ) RcFd ) RcFD,d + RcFV,d + ∑ l)1

RcFL,d + RcFTD,d + RcFd,w (6)

i.e., the sum of steady-drag, added-mass, lift, turbulent diffusion, and wall forces, respectively. All of the force terms are multiplied by the liquid fraction because of the reduced liquid volume available at considerable gas loadings. The net interfacial momentum-transfer term for the liquid phase (i.e., excluding the wall forces) can be written as N

∑

N

Mc,d ) -

d)1

∑ Md,c

(7)

d)1

The drag force is given by

FD,d ) -FD,d(vd - vc) ) -

3 R F C |v - vc|(vd 4ds,d d c D,d d vc) (8)

(11)

The acceleration of the liquid in the wake of the bubble can be taken into account through the added-mass force given by Ishii and Mishima103

N

RkFkg +

(10)

Fcds,d|vd - vc| µc

FV,d ) -RdFcCV,d

∂t

(9)

and ReB,d is the particle Reynolds number

(4)

In a consistent manner, the momentum balance for phase k yields

}

]

3E0,d 3A , ReB,d 8(E0,d + 4)

For pure systems, the parameter A ) 16, whereas for contaminated systems A ) 24. The Eo¨tvo¨s number, E0,d, is given as

N+1

∑ Rk ) 1 k)1

A (1 + ReB,d

(

)

Dvd Dvc Dt Dt

(12)

where CV,d ) 0.5 is derived for potential flow. The lift force on the dispersed phase due to shear in the liquid phase is expressed as104

FL,d ) -RdFcCL,d(vd - vc) × (∇ × vc)

(13)

where CL,d ) 0.5 is derived for potential flow. The dispersion of bubbles in turbulent liquid flow can be modeled as suggested by Carrica et al.13

FTD,d ) -

νc,t F ∇R RdSct,d D,d d

(14)

where the Schmidt number is defined as Sct,d ) νc,t/νd,t. An additional lift force that pushes the dispersed phase away from the wall was suggested by Antal et al.105 to be given by

(

FW,d ) max 0, Cw1 + Cw2

)

ds,d |vd - vc|2 RdFc nw (15) y ds,d

where Cw1 ) -0.1 and Cw1 ) 0.35. This force, eq 15, represents an extension of the original model of Antal et al.105 A defect in the original model, namely, that a bubble located far from the wall is attracted to the wall, has been removed.63 Turbulence Closures. It is expected that future multifluid models will adopt VLES turbulence closures.40 However, as noted in the Introduction, so far, the standard (single-phase) k- model is usually adopted as the basis for describing liquid-phase turbulence in Eulerian multiphase simulations. The structure of the turbulence model equations thus corresponds to the well-known generalized scalar transport equation for multiphase flow.39 A few extensions of the standard k- model have also been used, accounting for bubbleinduced turbulence. The turbulence models adopted in the sample simulations presented in this paper represent three versions

5118


of the standard k- model (or somehow including the standard k- model), as will be defined in further detail later. Initial and Boundary Conditions. The initial and boundary conditions used in the sample simulations presented in this paper are sketched below. Initially, there is no gas and no flow in the column (Rc ) 1, vk ) 0). The turbulent energy and dissipation rates are low but nonzero. When the gas starts entering the column, uniform void and velocity profiles for both phases are specified at the inlet. Neuman conditions are used at the other boundaries, except for that at the outlet, where a prescribed atmospheric-pressure boundary is applied. In addition, no inflow of gas is allowed. The inlet conditions for Rk and vk,z are determined using a 1D mass balance over the inlet plane, requiring that the relative velocity inside the column be determined by prescribed terminal velocities. Superficial velocities for the liquid and dispersed phases are prescribed. This approach neglects the complex bubble trajectories, added-mass and turbulent production mechanisms in the 2-3-cm inlet zone. The wall boundary conditions are based on the wall function approach consistent with the single-phase k- model. Assuming that the liquid velocity profile near the wall is similar to that of single-phase flow, the skin friction coefficient can be obtained from the logarithmic part of the “law of the wall”. At a distance y+ from the wall, the liquid wall shear stress can then be written as

τw ) -FcCf|vc|vc,z

(16)

The skin friction coefficient can be obtained from the “log-law”

(

)

〈vc,z〉y 1 1 ) ln ExCf ) ln(ExCf Rey) (17) κ ν κ c xCf 1

where y+ ) v*y/νc ) xCf〈vc,z〉y/νc ) xCf Rey and v* ) xCf〈vc,z〉. Here, κ ) 0.4 and E ) 9.973. As observed by Marie´ et al.,106 the wall exerts a friction force on the bubbles. A wall boundary condition for the gas momentum equations has been formulated by Lopez de Bertodano44 as

RdFc |v |v A Fd,w ) -Cwb,d 100ds,d c d,z w

(18)

∫V

(

)

R/d Fnd - Rnd Fnd dv ) ∆t

∫V∇‚(Rnd Fnd vnd) dv

The continuous-phase volume fraction is found by use of eq 4. 2. Temporary velocities are then found by solving the momentum balances without the interfacial coupling terms using explicit discretizations of all variables except for the volume fractions. The temporary volume fractions obtained in step 1 are used in all terms except for the convective terms

∫V

(

)

R/k Fnk v/k - Rnk Fnk vnk dv ) ∆t -

∫V∇‚(Rnk Fnk vnk vnk) dv + ∫V∇‚(R/k σj nk) dv ∫VR/k∇pn dv + ∫VR/k Fnkg dv

The governing two-fluid equations used in this contribution are discretized by a finite-volume algorithm. A staggered grid arrangement is used. A fractional-step scheme similar to those used by Tomiyama and Shimada54 and Lathouwers91 is developed on the basis of the single-phase algorithm reported by Jakobsen et al.107 and Lindborg et al.96 The fractional steps defining the algorithm are sketched in the following discussion. 1. Temporary volume fractions of the dispersed phase are first calculated by explicit discretization of the continuity equation

(20)

3. All momentum coupling terms, except for the lift force, are solved using semiimplicit time discretizations to avoid severe restrictions on the time step

∫V

(

)

/ n / R/k Fnk v// k - Rk Fk vk dv ) ∆t

and

∫V

(

)

/ dv ∫VFL,k

(21)

/ n // R/k Fnk v/// k - Rk Fk vk /// /// dv ) V(R/c FD,k + R/c FV,k + ∆t / /// / /// Rc FTD,k + Rc FW,k) dv (22)

∫

4. An equation for the pressure update was found by introducing a time-independent volume fraction (and density) step

∫V

(

)

R/k Fnk vn+1 - R/k Fnk v/// k k dv ) ∆t

∫VR/k∇φn+1 dv

(23)

where φn+1 ) pn+1 - pn. Applying the divergence operator to the terms in this relation and recognizing that the first term in the resulting relation represents an extended estimate of the transient term in the continuity equation (defined in step 5), yields

∫V

[

]

/ n /// ∇‚(R/k Fnk vn+1 k ) - ∇‚(Rk Fk vk ) dv ∆t ∂ - (RkFnk ) - ∇‚(R/k Fnk v/// k ) ∂t ) V dv ) ∆t

∫

[

-

where Cwb,d ) 1.0 is a common value. 3. Numerical Solution Algorithm

(19)

]

(24)

∫V∇‚(R/k∇φn+1) dv

Time discretization, followed by volume integration of the transient term and finally multiplication of the whole relation by ∆t/Fnk , yields

Rnk - Rn+1 k 1 ∆V - n ∆t F k

∫ ∇‚(R/k Fnk v/// k ) dv ) -

∆t Fnk

∫V∇‚(R/k∇φn+1) dv

(25)

Summing the corresponding equations for both phases and requiring that eq 4 be fulfilled both at time n and at n + 1, the equation for the pressure update is obtained as


-

∫

∫

1 1 / n /// ∇‚(R/d Fnd v/// d ) dv - n ∇‚(Rc Fc vc ) dv n Fd Fc )-

(26)

∫V∇‚(R/d∇φn+1) dv - ∆t ∫ ∇‚(R/c∇φn+1) dv Fn V

∆t Fnd

c

5. The void fractions are updated by solving the continuity equation for the dispersed phase using the newest estimate of the Rk and vk variables

∫V

(

)

Rn+1 Fnd - Rnd Fnd d dv ) ∆t

∫V∇‚(R/d Fnd vn+1 d ) dv

(27)

For variable-density flows, a suitable equation of state (EOS) is used to calculate the variations in the densities. The fractional-step concept applied consists of successive applications of the predefined operators determining parts of the transport equation. The convective and diffusive terms are further split into their components in the various coordinate directions. The timetruncation error in the splitting scheme is of first order. The convective terms are solved using a second-order TVD scheme in space and a first-order explicit Euler scheme in time. The TVD scheme applied was constructed by combining the central difference scheme and the classical upwind scheme by adopting the “smoothness monitor” of van Leer108 and the monotonic centered limiter.109 In the turbulence model, the diffusive terms were discretized with a second-order central difference scheme in space and a first-order semiimplicit Euler scheme in time. 4. Sample Results and Discussion In this section, the capabilities and limitations of 2D axisymmetric two-fluid models for simulating cylindrical bubble column reactor flows are discussed. Model limitations are explained, and sample results are given to illustrate the properties of the interfacial, turbulence, and wall interaction closures in use today. The aim is to evaluate whether a 2D model can be sufficient in providing reasonable predictions when descriptions of variable bubble size distributions, chemical reactive systems, continuous flow of the liquid phase, compressible effects, etc., are to be included in our fluid dynamic analysis. In this work, simulations were performed with three different turbulence closures for the liquid phase. A standard k- model was used to examine the effect of shear-induced turbulence (case a). In the first alternative, case b, both shear- and bubble-induced turbulence are taken into account by linearly superposing the turbulent viscosities obtained from the k- model and the model of Sato and Sekoguchi.55 A third approach, case c, is similar to case b in that both shear- and bubble-induced turbulence contributions are considered. However, in this model formulation, case c, the bubbleinduced turbulence contribution is included through an extra source term in the turbulence model equations.8 As mentioned earlier, it has been shown by Svendsen et al.,17 among others, that the time-averaged experimental data on the flow pattern in cylindrical bubble columns apparently becomes close to axisymmetric. Therefore, it was anticipated that steady-state 2D axisymmetric simulations could capture the pertinent time-averaged flow pattern in these columns and that a time-averaged flow description was sufficient for the

analysis of interest for chemical reaction engineers. Because 2D computations are much faster to perform, it is also important to learn as much as possible from such simulations before one performs dynamic 3D computations that actually might not provide any further information from a chemical engineering point of view. Therefore, in line with the early CFD analyses of bubble column flows, in this work, we are merely interested in the time-averaged behavior of the flow system, comparing the model results with time-averaged experimental data. As we proceed in discussing our results, we intend to reveal the pertinent reasons why this type of bubble column simulation has been performed for more than two decades with a rather limited degree of success. First, the lack of proper boundary conditions, as discussed in the Introduction, is further elucidated. At the column outlet, numerical problems arise for higher superficial gas velocities, because of the intense recirculating flow pattern produced by the free surface boundary that redistributes the liquid flow. Further numerical problems are observed for the case of operating the bubble column in a continuous liquid mode. The assumption of cylindrical axisymmetry in the computations prevents lateral motion of the dispersed gas phase and leads to an unrealistic radial phase distribution wherein the maximum void is away from the centerline (Figures 6 and 7), as also reported by other authors.35 Krishna and van Baten35 obtained better agreement with experiments by applying a 3D model. Even though the time-dependent terms in our 2D model were retained merely to ease the numerical solution of the 2D pseudo-steady-state problem, so that the resulting transients predicted by our 2D model are considered numerical modes, it is interesting to note that the flow development during the first 3-5 s of the simulations determines whether the fully developed steady-state flow pattern will result in wall, intermediate, or center void peaking (apparently multiple steady states). If the gas phase moves toward the wall or the center, there are no 3D secondary flows available to allow for any redistribution of the phases. Comparing our 2D axisymmetric model results at steady state with experimental data on the axial velocities and void fraction profiles at the axial level 0.3 m above the inlet, as in Figures 7 and 8, it can be observed that the magnitude of the simulated liquid axial velocity component is overestimated at the centerline in cases a and c, whereas the higher effective viscosity level in case b yields an estimation of the liquid axial velocity closer to the experimental results. The simulated axial velocities in the two phases are more tightly coupled than in the experiment, resulting in overestimation of the axial gas velocity component at the centerline and underestimation toward the wall. The corresponding void profile is in fair agreement with experiment, but with an unphysical local minimum at the center that has also been observed by other investigators. The cross-sectionalaveraged void, however, is overestimated by about 2535%, so the corresponding cross-sectional averaged gas velocity is underestimated. This is an indication that 3D effects are important at least in the inlet section. There is only one way to improve on this axisymmetric model limitation: a full 3D model is needed. However, studies reported in the literature show that it is very difficult to obtain reasonable time-averaged radial void profiles even in 3D simulations, even though the pre-

5120


Figure 6. Axial liquid velocity, axial gas velocity, gas voidage, and viscosity profiles as a function of column radius at z ) 2.0 m after 80 s (steady-state) using only drag and added mass as interfacial forces: ×, experimental data;76 s, case a; ‚‚‚, case b; - - -, case c. Grid resolution, 20 × 72; time resolution, 2 × 10-4 s.

Figure 7. Axial liquid velocity, axial gas velocity, gas voidage, and viscosity profiles as a function of column radius at z ) 0.3 m after 80 s (steady-state) using only drag and added mass as interfacial forces: ×, experimental data;76 s, case a; ‚‚‚, case b; - - -, case c. Grid resolution, 20 × 72; time resolution, 2 × 10-4 s.

dicted radial velocity profiles seem reasonable. Alternatively, or rather in addition, the deviation between the experimental data and the simulated results might also indicate that the inlet boundary conditions, phase distribution, and interfacial coupling are inaccurate and not in agreement with the real system, as will be discussed later. The present simulations indicate that the most important mechanism for radial movement of the gas, as the gas starts entering the bottom of the column, is the liquid wall friction (Figure 8). In our 2D pseudo-

transient simulations, the wall friction creates a significant pressure gradient at the inlet that pushes the gas toward the center line. An effect of the axisymmetric boundary condition is that it reflects most of the momentum in the fluid streams. Thus, if the effective viscosity in the bulk is low, only a small amount of the kinetic energy of the liquid will dissipate into heat, and the reflected fluxes of liquid and gas might have significant radial velocity components. In contrast, at the wall, the kinetic energy seems to be dissipated to a larger extent because of the wall friction. The radial


Figure 8. Gas and liquid velocity components and gas voidage profiles as a function of column radius showing center peaking in the bubble column after 4 s (left) and 8 0s (right). Grid resolution, 20 × 72; time resolution, 2 × 10-4 s.

Figure 9. Gas and liquid velocity components and gas voidage profiles as a function of column radius showing wall peaking in the bubble column after 4 s (left) and 8 0s (right). Grid resolution, 20 × 72; time resolution, 2 × 10-4 s.

boundary conditions thus make the radial flow unphysically sensitive to the turbulence level in the column. If the resulting radial velocity away from the center is high enough, a large fraction of the fluids seems to reach the wall, and the buoyancy force will ensure that the gas flows upward along the wall. This inlet flow pattern results in wall peaking, as shown in Figure 9. Alternatively, if the resulting radial velocity away from the center become lower because of a higher effective viscosity in the fluid phase, only negligible fluid streams will reach the wall. The liquid will then entrain the gas downward along the wall as the liquid continuity enforces a circulation cell, and a void center peaking profile occurs as shown in Figure 8. Apparently, similar numerical modes can be produced in dynamic 3D simulations with insufficient interfacial and turbulence closures. The limited accuracy reflected by the inlet boundary condition and the prediction of bubble size distribution phenomena is considered next. The specifications of proper inlet conditions for the phase distribution, local velocities, bubble size distributions, turbulence quantities are not trivial. The approach used here relies on the possibility of estimating an average bubble size or bubble size distribution at the inlet zone and a net terminal velocity and neglects the effects of the local

bubble jets found at the locations of the holes in the gas distributor. This approach thus induces a diffusive effect distributing the momentum uniformly over the column cross section. A better approach could be to resolve the flow through the gas distributor, as an integrated part of the model. Bertola et al.66 studied the effect of the bubble size distribution on the flow fields in bubble columns using a 2D axisymmetric mixture model with an inherent population balance. The initial bubble size was found not to be stable and was further determined by breakup and coalescence mechanisms that were to a large extent influenced by the local turbulence level. Even for the pure flow calculations, the common simplifying assumption of a uniform phase distribution at the inlet was found questionable, as the breakage and coalescence kernels rely on ad hoc empiricism determining the particle-particle and particle-turbulence interaction phenomena. For reactive systems, the estimates of the contact areas and thus the interfacial mass-transfer rate at the inlet are likely to contain large uncertainties. Turbulence model limitations are difficult to avoid even for single-phase flows. To enable multiphase simulations with the k- model, a few model adjustments were required to avoid divergence of the numer-

5122


ical algorithms. Accounting for shear-induced turbulence only, case a, a minimum turbulent energy level, kmin ) 5 × 10-3 m2/s2, was used in the first few seconds of the pseudo-dynamic simulation to ensure a sufficiently high viscosity level at the inlet to avoid having the turbulence level die out before the gas entered into the column. For the same reason, the contribution of the bubble-induced turbulence in case b was enforced by a factor of 10 during the first 3 s of the simulation until the k--model started to produce turbulent energy. In the case of introducing the bubble-induced turbulence through a source term in the turbulence model, case c, the start-up of the simulation was performed in the same way as for case a. Comparisons between the steady-state profiles at axial levels z ) 2.0 m and z ) 0.3 m for cases a-c are given in Figures 6 and 7, respectively. For z ) 2.0 m, all of the models give reasonable velocities compared to the experimental data. The turbulent viscosity of case b is approximately 20% higher than that of case a, which results in lower velocities that are closer to the experimental data. The turbulent viscosity profile in case c lies between the two other profiles. However, the velocities in case c are about equal to the velocities in case b. For z ) 0.3 m, there are larger discrepancies between the simulated results and the experimental data, as discussed earlier. The voidage profiles observed in the experiment are lower than at z ) 2.0 m, whereas the profiles from the simulations are approximately the same. The discrepancy indicates that the physics at the bottom of the column is not sufficiently captured by the 2D model. The interfacial coupling and wall forces are difficult to parametrize with sufficient accuracy. However, the steady-drag force is found to be the most important interfacial force in the system, as expected from literature observations. Reasonable agreement with experimental data is achieved by using this interfacial force only. A variety of different formulations of this force, and especially of the drag coefficient, can be found in the literature, valid for deformable and spherical particles and extended for different flow regimes, swarm effects, and pure and contaminated fluid systems. Kurul and Podowski111 compared several expressions for the drag coefficient in a boiling channel and showed that the different correlations for bubbly flows appeared to have minor effects on the predicted void fraction profiles. The different drag coefficients applied by Deen et al.,40 Krishna and van Baten,35 and Tomiyama53 deviate by less than 15% in our simulations, indicating the accuracy of this parameter. However, in several papers on bubble column modeling,14,34 drag coefficient correlations valid for nondeformable spheres,112,113 which give rise to drag force coefficients being about one-fourth the values of the other coefficients mentioned above for high Reynolds numbers, were adopted. When these spherical particle coefficient correlations were applied in our simulations, all of the profiles became flat. The gas is thus not able to entrain the liquid to the center in the bottom of the column to create circulation cells. The literature profiles, on the other hand, seem very similar to the experimental data. It is speculated that this deviation can be explained by the different grid arrangements used in the different codes. It is believed, however, that the zero void fraction occurring at the wall should be considered

a net effect of the forces acting on the bubbles, thus prescribing the wall void fraction results in an overconstrained set of equations, and should be avoided.105 Further examination of these issues is required to avoid any misinterpretations of the physical phenomena involved. Implementing the added-mass force has barely any influence on the steady-state solution; this conclusion is in agreement with the observations reported by Deen et al.,40 among others. Deen et al.40 explained this negligible effect by the fact that the simulations yield a quasi-stationary state in which there is only a little acceleration. The bubble jets observed close to the distributor plate are then not considered. However, the convergence rate and thus the computational costs are significantly improved when this force is implemented. The net effect of the transverse lift force is to push the gas radially inward or outward in the column depending on the sign of the lift coefficient. The force is induced by a velocity gradient and can therefore only reinforce or smooth out gradients in the flow fields. In our 2D simulations, this force alone cannot change a center peak flow pattern to a wall peak pattern or vice versa. This observation is quite surprising, as the flow pattern produced by the early steady-state model of Jakobsen20 was totally dominated by this force. A reasonable explanation could be that the simulations by Jakobsen20 were strongly influenced by numerical diffusion, as the resolution in the axial direction was limited. Effects similar to those of the lift force are observed when a turbulent dispersion force that is proportional to the gradient of the gas voidage is implemented. The only effect seen is that it smoothes out sharp gradients in the domain until the profiles become completely flat, as the models tested seemed to overestimate the diffusion effect. This trend is probably related to the low accuracy reflected by the effective viscosity variable determined by the turbulence model. The modeling of this force should also be seen in connection with the inherent numerical diffusion caused by the numerical approximation of the convective terms. The radial wall lift force proposed by Antal et al.105 requires a certain minimum resolution of the grid close to the wall. The choice of model coefficients is crucial for the development of the profiles as these parameters determine the magnitude of the force as well as the operative distance from the wall. The magnitude of the force is highest for larger bubbles, and it seems to be very sensitive to the bubble size. However, in our 2D model simulations, the force could not be used to alter the flow development significantly except very close to the wall. There was no mechanism available for transporting the gas further toward the center. Lopez de Bertodano44 correctly suggested that several unresolved phenomena close to the wall influence the phase distribution. A net axial friction force acting on the gas in the vicinity of the wall was derived. In the present work, the net effect of this force is to reduce the time scale for the flow development from the gas starts entering the column. However, this force has a significant effect on the solution when the drag coefficient for nondeformable spheres is used, as it causes center peaking, whereas the profiles become flat without this force as mentioned above. This indicates that a noslip condition specified for the gas phase as well as a fixed zero void fraction on the wall, which is used


as boundary conditions in some numerical codes, will give reasonable results even with this drag coefficient. Another limiting model simplification adopted in most reports on bubble column modeling is the assumption of a constant gas density. Simulating columns that are several meters high, adopting an ambient pressure boundary at the outlet, cannot be consistent with this simplifying assumption, as density gradient effects become important. Similar inconsistency problems arise when reactive systems are simulated, as the gas density has to be allowed to vary in accordance with the chemical process behavior. In addition, in our experience, interfacial mass-transfer fluxes and chemical reactions often induce numerical stability problems and the lack of proper convergence. Finally, multiphase flow simulations are numerically unstable considering both 2D and 3D models because of the applications of low-accuracy discretizations, large density differences between the phases, ill-conditioned implementations of interfacial closures in the limit of zero void or zero liquid fractions, and large gradients or discontinuities in the flow fields that can occur in these flow systems. Therefore, further studies on these issue are highly recommended. Preliminary results using our 2D model show that accurate discretization schemes for the dispersed-phase continuity equation and for the convection terms within the fluid momentum equations are crucial for accurate predictions of the flow development. Large discretization errors can change the flow development from center to wall peaking or visa versa. Similar observations were reported by Oey et al.14 regarding the simulation of 2D rectangular columns in 3D. The physical interpretation and relevance of dynamic 2D analysis of bubble column flows have been matters of debate for many years. Judging from the experimental analyses discussed in the Introduction, it is evident that bubble column flows should be considered dynamic and 3D in any rigorous fluid dynamic analysis. Turbulence is a dynamic 3D phenomenon. However, direct numerical simulations (DNS) of these systems is still not feasible, so some kind of filtering or averaging is needed. Dynamic 3D simulations by use of filtered VLES models are computationally intensive and perform accurately only when the cutoff turbulence length scale (eddy size) is placed at a spectral gap in the turbulence energy spectrum, thereby ensuring that the energy flux between the resolved and subgrid scales is very small. Accurate modeling of these fluxes itself is thus necessary only when the subgrid-scale fluxes get large in comparison with the resolved fluxes (as in most VLES simulations114). Unfortunately, in engineering practice, the cutoff length scale is often apparently chosen arbitrarily as no energy gap is found in the energy spectrum, requiring more sophisticated subgrid closures. This might be one of the reasons why most dynamic 3D simulations of bubble column flows have been performed using Reynolds-averaged models. Among the Reynolds-averaged turbulence models available, the simple k- model is adopted in almost all of the investigations reported in the literature. There is experimental evidence, however, that the flow in bubble columns is highly anisotropic.115 It is thus an open question whether a 3D simulation using a k- model (predicting isotropic turbulence) provides the necessary 3D effects or whether perhaps Reynolds stress or VLES

models are required to represents the dynamics of these flows. However, because of the impact of phenomena that cannot be resolved by continuum models, it is not expected that even the phase distributions predicted by the use of VLES models will always be physically realistic. Furthermore, when unsteady turbulent flow simulations using Reynolds-averaged models are performed, providing physical interpretations of the model results is not trivial. Schlichting and Gersten116 suggested that these flow variables could be interpreted as timeaverage quantities whereby the time interval used in the averaging operator had been chosen large enough to include all turbulent fluctuations but still small enough to contain no effect from the transient, or periodic, parts. This concept is based on the assumption that the transient process occurs very slowly or that the frequency of the oscillation is very small and lies outside the turbulent spectrum. Telionis117 discussed an extension of the standard Reynolds averaging approach: a triple decomposition, where the instantaneous quantities are decomposed into three parts. Recall that, in the standard Reynolds decomposition, we split the instantaneous quantities into a time-independent time average (mean) and a fluctuation, whereas in the triple decomposition approach, the mean motion is also dependent on time. The mean is thus generally made up of a time-independent part and a time-dependent part. In other cases, the application of the standard Reynolds averaging concept has been further extended to simulations of transients that are within the turbulent spectrum. To compare the simulated results obtained with this type of model with experimental data that are averaged over a long time period to give steady-state data (representing the whole spectrum of turbulence), both the modeled and the resolved scales have to be considered in much the same manner as for VLES.118 For bubbly flows, this model interpretation induces another severe problem related to the bubble-induced turbulence phenomena. For steady-state simulations using the standard Reynolds-averaged models, the turbulent kinetic energy variable, k, denotes a measure of the mean energy considering all time scales within the flow (i.e., for the whole turbulence energy spectrum). In a typical energy spectrum, most of the energy is accumulated on the larger scales of turbulence and very little on the smaller scales. This k quantity is thus sometimes used as a measure of the energy level on the integral scales (i.e., the larger energy containing scales of turbulence). Considering the energy spectrum in terms of turbulence length scales, the energy-containing scales represented by the k- model are thus normally much larger than the bubble size. The inclusion of turbulence production due to the bubbles’ relative motion is therefore based on the assumption of an inverse cascade of turbulence.8 For dynamic simulations, the k quantity can represent scales less than or at the same order of magnitude as the particle size. In the cases where the time averaging operator is chosen in such a way that the dispersed phase is resolved, no bubble-induced turbulence mechanisms should be included in the turbulence model. In summary, one can certainly conclude that dynamic simulations of bubble column flows are not a trivial task. Physical interpretations of the simulated 3D dynamics should be performed with caution when any kind of turbulence model is adopted.

5124


Considering the much simpler approach used in our simulations, it was found that the 2D steady-state bubble column model has very limited inherent capabilities of predicting bubble column flows with sufficient accuracy. However, after tedious adjustments of model parameters, boundary conditions, time and space resolutions, numerical implementations, and solution algorithms, it is often possible to reproduce known flow fields to a certain extent. The global liquid flow pattern is usually captured fairly well, whereas the corresponding gas fields are more questionable and difficult to evaluate because of the lack of reliable experimental data. Unfortunately, the possibility of achieving a completely erroneous steady-state flow pattern where the liquid flows up along the wall and down in the core of the column seems ever-present. The almost definitive rejection of this reactor model formulation is related to the fact that it is often difficult to reproduce the gas volume fraction fields with sufficient accuracy. The latter limitation would severely affect the contact areas, the interfacial heat- and mass-transfer fluxes, and the projected areas, as well as the interfacial momentumtransfer fluxes, making reliable predictions of the performance of chemical processes impossible. It is speculated that the more specific model limitations that have to be eliminated to enable reliable model predictions can be found among the topics outlined in the following. On the physical side, the model has to be made 3D and dynamic, proper initial and boundary conditions have to be formulated, the high-frequency turbulence mechanisms and lower-frequency coherent structures have to be resolved (requiring sufficiently accurate turbulence closures), the fluid particleturbulence interaction mechanisms as well as interparticle interactions are probably important, the description of the interactions between the dispersed phase and the turbulent wall boundary layer structure should be extended, the pertinent phase distribution mechanisms have to be identified and described with sufficient accuracy, bubble size and shape distributions have to be considered (i.e., the microscopic bubble coalescence and breakage processes have to be described), the interfacial momentum closures have to be described with higher accuracy (especially the transverse forces), the free surface at the top of the column might have to be resolved, the mechanical and turbulence energy balances in the system have to be fulfilled (for validation of the interfacial transfer fluxes and the energy dissipation rates), phasic density variations should be allowed, and the impact of heat- and masstransfer mechanisms and chemical conversion might need to be investigated in further detail. On the numerical side, any algorithms, grid arrangements, unphysical source implementations (e.g., unphysical model extensions implemented to avoid ill-conditioned implementations of interfacial and turbulence closures in the limit of zero void fractions), or discretization truncation errors producing solutions merely determined by numerical modes rather than physical mechanisms should be removed (if possible). Finally, proper convergence criteria should be formulated and fulfilled for every time step in the simulation. In view of the comprehensive list of limitations provided above, it is concluded that a truly predictive model is not imminent. However, the more pressing improvements required soon are those that will enable a sufficiently accurate description of the time-averaged

phase distribution in the system. These improvements should also ensure that one would not achieve, on a time-averaged basis, downward flow of liquid in the core of the column and upward flow of liquid close to the wall that completely disagree with the available experimental observations. The authors believe that, in this context, one of the most important model extensions required is the inclusion of a proper description of the bubble size and shape distributions. Further work in our group continues to elucidate these issues along the lines described in the remaining sections of this paper. 5. Population Balance Framework The chemical engineering community began the first efforts of association with the concepts of the population balance in the 1960s. Two fundamental modeling frameworks emerge for formulating the early population balances, in quite the same way as the kinetic theory of gases and the continuum theory were proposed deriving the governing conservation equations in fluid mechanics (e.g., appendixes in Williams119). A third less rigorous approach is also used in formulating the population balance directly on the volume-averaging scale considering the internal coordinates (i.e., particle size), an analogue to multiphase mixture models. Considering dispersed two-phase flows, a few research groups have thus preferred to describe the dispersed phase on the basis of a statistical Boltzmann-type equation determining the temporal and spatial rates of change of a suitably defined distribution function. Performing “Maxwellian” averaging, integrating all terms over the whole velocity space to eliminate the velocity dependence, one obtains the generic population balance equation. Reliable closures are required for the unknown terms resulting from the averaging process. However, this theory provides rational means of understanding the one-way coupling between discrete particle physics and the average continuum properties. The procedure sketched above has much in common with the granular theory (i.e., that could be defined as the flow of powder in a vacuum) of solid particles.120 However, the great majority of research groups within the chemical engineering community adopted an approach based on an extension of the classical continuum theory instead developing the population balance equation. The main reason for this choice was probably that the basic concept is familiar to most chemical engineers from basic courses in fluid mechanics, whereas the principles of kinetic theory might be less accessible from a pedagogic point of view. On the other hand, the continuum theory gives only an average representation of the dispersed phase (i.e., on a macroscopic level relative to the particle scale) and does not provide any information on the unresolved mechanisms formulating the population balance closures. In most cases, the balance principle is applied by formulating the particle number balance on the integral form. Extended versions of the Leibnitz and Gauss theorems are then required to transfer the integral balances to the differential form. For turbulent flows, the governing microscopic continuum equations on the differential form are then time(Reynolds-) averaged after being volume averaged to obtain trackable volume and time resolutions, giving rise to additional closure requirements.121 The third group of population balances are thus formulated directly on the averaging scales.The Leibnitz and Gauss


theorems are used to cast all of the terms in the equation into a volume integral. The governing integral equation is thereafter converted to the differential form. In this formulation, the closures are purely empirical parametrizations based on intuitive relationships rather than sound scientific principles. Even though the kinetic theory approach is considered more general, the basic balance formulation might need to be reformulated for every novel system considered to ensure that the pertinent physics on the discrete particle level is consistent with the problem in question. An optimized choice of distribution function definition might be necessary, and novel closures might be required, making the approach rather demanding theoretically. The population balance concept was first presented by Hulburt and Katz.122 Rather than adopting the standard continuum mechanical framework,123 the model derivation was based on the alternative continuum approach based on particle statistics familiar from classical statistical mechanics. The main problems investigated stem from solid particle nucleation, growth, and agglomeration. Randolph124 and Randolph and Larson,125 on the other hand, formulated a generic population balance model based on the extended continuum mechanical framework. Their main concerns were solid particle crystallization, nucleation, growth, agglomeration/ aggregation, and breakage. Ramkrishna126,127 adopted the concepts of Randolph and Larson to explore biological populations. An outline of the population balance model derivation from the continuum mechanics point of view was discussed. Similar approaches are also frequently used in the theory of aerosols in which the gas is the continuous phase121,128 and in chemical, mechanical, and nuclear engineering describing multiphase droplet flow dynamics;129,130 such an approach was also used for incompressible bubbly two-phase flows by Guido Lavalle et al.131 Carrica et al.13 developed the population balance from kinetic theory using the particle mass as the internal variable while investigating compressible bubbly twophase flow around a surface ship, whereas most earlier work on solid particle crystallization used particle volume (or diameter). The two formulations are completely equivalent in the case of incompressible gases. However, in flows where compressibility effects in the gas are important (as in the case of experimental bubble columns operated at atmospheric conditions), the use of mass as an internal variable was found to be advantageous because it is conserved under pressure changes. Lehr et al.,32 Millies and Mewes,68 Lehr and Mewes,36 and Pilon et al.132 sketched a possible alternative formulation using particle volume (diameter) as the internal coordinate. In their approach, several growth terms have to be considered to express the effects of gas expansion due to changes in gas density in accordance with a suitable EOS. The use of a mass density form of the population balance derived according to the kinetic theory approach (i.e., instead of the more common number density) has also been discussed, having several advantages in reactor technology.133,134 In chemical engineering, Coulaloglou and Tavlarides135 were among the first to introduce the simpler macroscopic formulation, describing the interaction processes in agitated liquid-liquid dispersions. Drifting

from the fundamental equations, the closures became an integrated part of the discrete numerical discretization scheme adopted (i.e., there is no clear split between the numerical scheme and the closure laws). Lee et al.136 and Prince and Blanch77 adopted the basic ideas of Coulaloglou and Tavlarides135 to formulate the population balance source terms directly on the averaging scales to perform an analysis of bubble breakage and coalescence in turbulent gas-liquid dispersions. The source term closures were completely integrated parts of the discrete numerical scheme adopted. The number densities of the bubbles were thus defined as the number of bubbles per unit mixture volume and not as a probability density in accordance with the kinetic theory of gases. Luo and Svendsen71 extended the work of Coulaloglou and Tavlarides,135 Lee et al.,136 and Prince and Blanch,77 formulating the population balance directly on the macroscopic scales where the closure laws for the source terms were integrated parts of the discrete numerical scheme used to solve the model equations. The works of Millies and Mewes,68 Lehr and Mewes,36 Lehr et al.,32 Hagesaether et al.,73-75 and Wang et al.137 all adopted the governing population balance of Luo and Svendsen.71 The extension of these modeling concepts to bubbledriven gas-liquid flows is not as straightforward as one might think from either a physical or a computational point of view. Many physical or physicochemical mechanisms are still not understood or possibly even not discovered yet, and as the pertinent interfacial phenomena involved are expected to interfere at the molecular scale, tremendous computational efforts are required. Further challenges are foreseen considering reactive systems in industrial chemical reactors, as the changes in species compositions, heats of reaction, and surface processes are strongly linked. However, preliminary attempts have already been made in the modeling of bubble coalescence and breakage in gas-liquid systems, adopting the population balance framework and mathematical modeling tools proposed by the above-mentioned pioneers. Integrated multiphase flow/population balance models have been applied to bubble column reactors and two-phase stirred vessels with limited success.31,32,41,65,66,138,139 The simulations reveal that these models are not able to predict nonequilibrium bubble size distributions sufficiently accurately as a function of space and time because of the very restricted predictive capabilities reflected by the present population balance kernels for both bubble coalescence and breakage. In the following subsections, two forms of the population balance are formulated in accordance with the continuum theory. First, a macroscopic balance is formulated directly on the averaging scales in terms of number density functions. A corresponding set of source term closures is presented as well. Second, a more fundamental microscopic population balance as well as fairly general expressions for the source terms are formulated in terms of number probability densities. This balance is then averaged using the popular timeafter-volume averaging procedure. The fairly general form of the constitutive relations adopted for the bubble coalescence and breakage phenomena are formulated purely on the basis of physical reasoning and intuitive interpretations of the mechanisms involved. Further links to the discrete particle-scale phenomena are usually expressed on the averaging scales by extrapo-

5126


lating relationships and concepts from the kinetic theory of gases. However, the basic kinetic theory of gases does not provide any coupling between the particles and the ambient fluid. In a more generalized framework, however, the kinetic theory concept allows for the introduction of external forces determining the interaction between the fluid and the dispersed particles. The approach described by Ramkrishna127 is outlined. Macroscopic Population Balance Equation. In this subsection, the macroscopic population balance formulation of Luo70 is outlined. In the work of Luo,70 no growth terms were considered; the balance equation thus contained a transient term, a convection term, and four source terms due to binary bubble coalescence and breakage. The population equation was expressed as

∂ni + ∇‚(vini) ) BB,i - DB,i + BC,i - DC,i ∂t

( )

1 s m3 (28)

where ni is the number density with units of 1/m3 and N vi is the mass-average velocity vector, vi ) (∑i)1 N nivi,pFgVi)/(∑i)1niFgVi). [Note that, alternatively, the N N nivi,p)/(∑i)1 ni); particle number-average, vi,ni ) (∑i)1 N N size-average, vi,di ) (∑i)1nivi,pdi)/(∑i)1nidi); surface-avN N nivi,pai)/(∑i)1 niai); or volume-average, erage, vi,ai ) (∑i)1 N N vi,Vi ) (∑i)1nivi,pVi)/(∑i)1niVi) velocity vector could be used as well.140] The source terms express the bubble number birth and death rates per unit dispersion volume for bubbles of size di at time t due to coalescence and breakage. The source terms are assumed to be functions of bubble size di, bubble number ni, and time t. The birth of bubbles of size di due to coalescence stems from the coalescence between all bubbles of size smaller than di. Hence, the birth rate for bubbles of size di, BC,i, can be obtained by summing all coalescence events that form a bubble of size di. This gives di/2

BC,i )

∑ d )d j

ΩC(dj:di - dj)

min

( ) 1

s m3

(29)

where dmin is the minimum bubble size and depends on the minimum eddy size in the system. The source term definition implies that bubbles of size dj coalesce with bubbles of size (di - dj) to form bubbles of size di. The upper limit of the sum stems from symmetry considerations or avoidance of counting the coalescence between the same pair of bubble sizes twice. Similarly, the death of bubbles of size di due to coalescence stems from coalescence between two bubbles in class di or between one bubble in class di and other bubbles. Hence, the bubble death rate for bubbles of size di, DC,i, can be calculated as dmax-di

DC,i )

∑ d )d j

min

ΩC(dj:di)

( ) 1

s m3

(30)

where dmax is the maximum bubble size in the system. The upper limit indicates that the bubble size formed by coalescence must not exceed dmax. The birth of bubbles of size di due to breakage stems from the breakage of all bubbles larger than di. The breakage birth rate, BB,i, can be obtained by summing all breakage events that form bubbles of size di

dmax

BB,i )

∑ ΩB(dj:di) d )d j

i

( ) 1

s m3

(31)

The death of bubbles of size di due to breakage stems from breakage of bubbles within this class; thus

DB,i ) ΩB(di)

( ) 1 s m3

(32)

The local gas volume fraction in any grid cell can be N calculated as Rg ) ∑i)1 ni(π/6)di3. The macroscopic model formulation has numerical properties similar to those of the discrete discretization scheme (discussed later) and can be solved directly without further considerations. To ensure number and mass conservation, Hagesaether et al.73-75 extended this model by adopting a numerical procedure to redistribute the bubbles on pivot points in accordance with the discrete solution method. This modification of the solution strategy was introduced because of an inherent model limitation, i.e., the breakage closure of Luo and Svendsen71 was not necessarily conservative, as will be explained later. Macroscopic Source Term Closures. In accordance with the work of Coulaloglou and Tavlarides135 and Prince and Blanch,77 Luo70 assumed that all macroscopic source terms determining the death and birth rates could be defined as the product of a collision density and a probability. Thus, modeling of bubble coalescence means modeling of a bubble-bubble collision density and a coalescence probability, whereas modeling of bubble breakage means modeling of an eddy-bubble collision density and a breakage probability. Models for the collision densities were derived assuming that the mechanisms of the bubble-bubble and eddy-bubble collisions are analogous to those of collisions between molecules as in the kinetic theory of gases.135 Models for the Binary Bubble Coalescence Rate, ΩC(di:dj). For coalescence between bubbles of class, di, and bubbles of class, dj, the coalescence rate is expressed as

ΩC(di:dj) ) ωC(di:dj) pC(di:dj)

( ) 1 s m3

(33)

Models for the Bubble-Bubble Collision Density, ωC(di:dj). Although not necessarily valid for bubble columns, the dispersions are considered sufficiently dilute that only binary collisions need to be considered. A collision of two bubbles can occur when the bubbles are brought together by the surrounding liquid flow or by body forces such as gravity. At least three sources of relative motion can be distinguished: motion induced by turbulence in the continuous phase; motion induced by mean velocity gradients; and motion induced by buoyancy (or, more generally, body forces) arising from different bubble slip velocities, wake interactions, or helical/zigzag trajectories.140 However, most studies on bubble columns have been restricted to considering only models for the contribution of turbulence to coalescence; the contributions from mean velocity gradients and body forces are generally neglected without proper validation. These collision mechanisms are generally both signifi-


Figure 10. Sketch of a collision tube of an entering bubble moving through the tube with a velocity vp. The bubbles within the tube are assumed to be “frozen” or stationary.

cant, which greatly complicates the construction and validation of coalescence models. The parametrization of the collision densities was obviously influenced by the more fundamental studies on formulating the source terms at the microscopic level followed by some kind of averaging and a discrete numerical discretization scheme, as will be further outlined in a later section. At this point, it is simply observed that the units of the volume-average coalescence frequency correspond to the units of the effective swept volume rate in the kinetic theory of gases. Therefore, rather than the actual collision density, the volume swept by a moving particle (i.e., bubble or eddy) multiplied by the particle density squared is used to define the collision densities.77,135,138 To derive an expression for the effective swept volume rate, one usually considers a particle as it travels in a straight path from one collision to the next in a monodisperse dispersion, as sketched in Figure 10. The particle’s speed and direction of motion changes with each collision. Further imagine that, at a given instant, all particles except for the one in question are frozen in position and this particle moves with an average speed, vp. At the instant of collision, the center-to-center distance of the two particles is d. The collision cross section of the target area of the particle is σA ) 1/4πd2. It should be noted that the collision cross-sectional area is defined in several ways. In one approach, inspired by kinetic theory, the moving particle is approximated as a point-like particle. Thus, the cross section of the moving particle is not considered in calculating the effective cross-sectional area; hence, σA ) 1/4πd2. In the case of small particles, this approximation might be sufficient. However, an improved estimate of the effective cross-sectional area is achieved by taking into account the size of the moving particle as well; thus σA ) πd2, as sketched in Figure 11. In time ∆t, the moving particle sweeps out a cylindrical volume of length vp∆t and cross section 1/4πd2. Any particle whose center is in this cylinder will be struck by the moving particle. The number of collisions in the time ∆t is f1(1/4πd2vp∆t), where f1 is assumed to be locally uniformly distributed in space. The effective swept volume rate hC(d) is given by

1 hC(d) ) σAvp ) πd2vp 4

( ) m3 s

(34)

where σA ) 1/4πd2 is the cross-sectional area of the socalled “collision tube”.

Figure 11. Sketch of a collision tube of a bubble, defining the collision diameter and the effective cross-sectional area.

The collision density of a single particle is defined as the number of collisions per unit time and length (diameter)

[ ]

1 ωf1(d) ) f1(d) hC(d) ) f1(d,t) πd2vp 4

1 s (m)

(35)

This relation represents a very crude collision model for a dispersion containing only one type of particle. Sometimes, the collision density representing the number of collisions per unit mixture volume and per unit time and length (diameter) is a more convenient quantity. Multiplying the collision density of a single particle by the particle number density, one obtains the modified collision density

[

1 1 ωf1f1(d) ) hC(d) f1f1 ) f1f1 πd2vp 2 4

]

1 m s (m) (m) (36) 3

The factor of 1/2 appears because h(d) represents twice the number of collisions. In accordance with the pioneering work of Smoluchowski,141 similar considerations can be repeated for dispersions containing two types of particles having diameters d, d′ and particle densities f1, f′1. The resulting collision density is given by

ωC(d,d′) ) hC(d,d′) f1f′1 ) π d + d′ 2 vrel f1f′1 4 2

(

)

[

]

1 (37) m s (m) (m) 3

where the collision frequency is calculated using the mean collision diameter d ≡ (d + d′)/2, as sketched in Figure 12. Note that, in this case, the 1/2 factor is not included. Kolev130 discussed the validity of these relations for fluid particle collisions considering the obvious discrepancies resulting from the different nature of the fluid particle collisions compared with the random molecular collisions. The basic assumptions in kinetic theory that the molecules are hard spheres and that the collisions are perfectly elastic and obey the classical conservation laws do not hold for real fluid particles because these particles are deformable and elastic and can agglomerate or even coalesce after random collisions. The collision density is thus not truly an independent function of the coalescence probability.

5128


subrange of turbulence, then δv2(d) ) C(d)2/3. A relationship between the constant C and the Kolmogorov constant Ck (i.e., the parameter in the energy spectrum function for the inertial subrange) can be deduced from the power-law spectrum and the second-order structure function, C ) (27/55)Γ(1/3)Ck.144,145 Setting Ck ) 1.5 in accordance with experimental data yields C ≈ 2.0. Interpreting the second-order structure function (i.e., because it is a two-point correlation function) as a relative particle velocity (i.e., actually a mean velocity obtained using one-point statistics seems more appropriate), the collision density suggested by Prince and Blanch77 becomes

ωTC(di:dj) ≈ ninj Figure 12. Sketch of the mean collision diameter and the effective cross-sectional area.

As mentioned above, bubble-bubble collisions can occur as a result of a variety of mechanisms. Prince and Blanch77 modeled bubble coalescence in bubble columns considering bubble collisions due to turbulence, buoyancy, and laminar shear through analysis of the coalescence probability (efficiency) of collisions. It was assumed that the collisions from the various mechanisms were cumulative. The collision density resulting from turbulent motion was expressed as a function of bubble size, concentration, and velocity in accordance with the work of Smoluchowski141

(

)

π d + d′ 2 ωT(d,d′) ≈ f1f′1 (vj t,d2 + vj t,d′2)1/2 4 2

[

]

1 (38) m s (m) (m) 3

Estimates of the turbulent velocities have been obtained using certain relations developed in the classical theory on isotropic turbulence due to Kolmogorov. 142 At this point, the population balance closures seem to become rather vague, as the application of the classical turbulence relations developed for fluid velocity fluctuations or vortices to describe flows of discrete fluid particles are not always in accordance with the basic turbulence theory restrictions. The theory of Kolmogorov142 states that, if the distance λ between two points in the flow field is much smaller than the turbulence macroscale L but much larger than the Kolmogorov microscale λd, then the second-order velocity structure function is a function of only the turbulent energy dissipation rate and the magnitude of the length λ, δv2(λ) ) C(λ)2/3. This two-point velocity structure function should not be confused with the normal component of the Reynolds stresses, which is a one-point velocity correlation function. A summary of the theory involved is given by Pope143 (Chapter 6 and Appendix G). By definition, the secondorder velocity structure function is the covariance of the difference in velocity between two points in physical space. Considering the second-order structure function defined by Kolmogorov142 for the continuous fluid flow turbulence, it can be expressed at a separation equal to the bubble diameter d as δv2(d) ) [vz(z+d)-vz(z+d)]2. If the magnitude of diameter d lies within the inertial

(

)

π di + dj 2 (vj t,di + vj t,dj)1/2 ≈ 4 2

0.089πninj(di + dj)21/3(di2/3 + dj2/3)1/2

( )

1 (39) s m3

(Note that the collision density was originally expressed in accordance with the macroscopic population balance.70 If otherwise acceptable, it seems that the coalescence closures formulated directly within the macroscopic framework can be transformed and expressed in terms of probability densities and used within the average microscopic balance framework without further considerations if one adopts a discrete numerical discretization scheme.) Note that the parameter value used by Luo70 and Luo and Svendsen71 is a factor of 4 larger than this, one as they used a different definition of the effective crosssectional area. In this interpretation, two bubbles being dispersed in a large-scale fluid eddy (i.e., larger than the bubble scale, d) are assumed to be advected together with the fluid without approaching each other. However, the validity of the Kolmogorov hypothesis describing bubbly flows is not yet strictly verified. For convenience, the relationship for turbulent fluid velocity fluctuations has thus been interpreted and extrapolated in many ways. There seems to be some confusion in the literature regarding different interpretations of the Kolmogorov velocity scale and the second-order structure function.143 The second-order structure function, δv2(d), has thus (apparently erroneously) been used as an estimate for an imaginary absolute turbulent bubble velocity in the inertial subrange of isotropic turbulence, in other cases considered a relative bubble velocity in a turbulent flows,77 and as a third possibility been used a measure of intrabubble oscillations. The buoyancy collision density ωBC(di:dj) is expressed as (Prince and Blanch;77 Williams and Loyalka,128 p 164)

ωBC(di:dj) ≈ ninj

(

)

π di + d j 2 |vj r,di - vj r,dj| 4 2

( ) 1 s m3

(40) where vjr,di is the rise velocity of the particle. (Again, note that the collision density was originally expressed in accordance with the macroscopic population balance formulation.70) The functional form of the collision rate due to laminar shear is given by (Friedlander,121 p 200; Williams and Loyalka,128 p 170)


( ) ( )

dvc 4 3 ωLS C (di:dj) ≈ ninj (di + dj) 3 dR

1 s m3

(41)

where vc is the continuous-phase circulation velocity and R is the radial coordinate of the column. The term (dvc/dR) is the average shear rate. [Again, the collision density was originally expressed in accordance with the macroscopic population balance formulation (e.g, Luo70).] The net coalescence frequency of bubbles of diameters di and dj was then calculated by superposing the different bubble collisions mechanisms as a linear sum of contributions multiplied by a common efficiency (a similar closure can be obtained for the average microscopic formulation if one adopts a discrete numerical discretization scheme, as will be described later)77

ΩC(di:dj) ) [ωTC(di:dj) + ωBC(di:dj) + ωLS C (di:dj)]pC(di:dj)

( )

1 (42) s m3

Hagesaether146 adopted this approach when modeling bubble column dispersions and found that, with the choice of parameter values used in his analysis, the turbulent contribution dominated the collision rate for the bubbles in the system. However, the most important model parameter as the turbulent energy dissipation rate is difficult to compute, and only crude estimates were obtained.32,70,77 In addition, Kolev130 argued that the frequency of coalescence of a single bubble should rather be determined for individual efficiencies

ΩC(di:dj) ) ωTC(di:dj) pTC(di:dj) + ωBC(di:dj) pBC(di:dj) + 1 LS (43) ωLS C (di:dj) pC (di:dj), s m3

( )

Kolev130 further assumed that, for the coalescence processes induced by buoyancy and nonuniform velocity fields, the forces leading to collisions inevitably act toward coalescence for contracting particle free paths. Therefore, the probabilities must have the following ranges: 0 E pTC(di:dj) E 1 and pBC(di:dj) ) pLS C (di:dj) ≈ 1. However, no firm conclusions on the relative importance of the various contributions have yet been drawn, as the turbulent particle velocity closures are at best inaccurate. Models for the Probability of Coalescence, pC(di:dj). The probability or efficiency of oscillatory bubble coalescence, e.g., induced by turbulent fluctuations, is expected to be determined by physical mechanisms on various scales. The coalescence process is assumed to occur in three consecutive stages. First, bubbles collide, trapping a small amount of liquid between them under the action of the continuous phase. Second, this liquid then drains over a period of time from its initial thickness until the liquid film separating the bubbles reaches a critical thickness, under the action of the film hydrodynamics. The hydrodynamics of the film depends on whether the film surface is mobile or immobile, and the mobility, in turn, depends on whether the continuous phase is pure or a solution. Third, at this point, film rupture occurs because of film instability, resulting in instantaneous coalescence (i.e., usually not modeled in further detail). The probability of coalescence is thus defined as a function of the ratio of the time interval within which the bubbles are touching each other, called

the collision time interval, ∆tcol, and the time interval required to push out the surrounding liquid and to overcome the strength of the capillary microlayer between the two bubbles, called the coalescence time interval, ∆tcoal130

pC ∝ f(∆tcoal/∆tcol)

(44)

The coalescence efficiency thus represents the fraction of particles that coalesce out of the total number of particles that have been colliding. The functional relationship given by Coulaloglou and Tavlarides,135 i.e.

pC ≈ exp(-∆tcoal/∆tcol)

(45)

is often used for the modeling of bubble coalescence probabilities for bubble column flows. The duration of such interactions is limited, and coalescence will occur only if the intervening film can drain to a sufficiently small thickness to rupture in the time available. Coalescence is obviously possible if there is enough contact time for completing the coalescence

∆tcol > ∆tcoal

(46)

Luo,70 Luo and Svendsen,71 and Hagesaether et al.73-75 adopted this approach. The complexity of the film draining phenomena involved is a severe problem and might be best illustrated by briefly discussing an early modeling attempt. Oolman and Blanch147 derived an expression for the coalescence times in stagnant fluids by examining the time required for the liquid film between bubbles to thin from an initial thickness to a critical value where rupture occurs. The original model considers the flow rate of fluid from the liquid film by capillary pressure, augmented by the Hamaker contribution (reflecting the mutual attraction of fluid molecules on opposite sides of the liquid film) at very low film thicknesses, bubble deformation, and the changes in the concentration of surfactant species

-

{ [ ( ) (

8 -4c dσI dh ) dt Rd2Fc RT dc

2

+ h2

)]}

2σI A + rb 6πh3

1/2

(47)

where h is the film thickness, Rd is the radius of the liquid disk between the coalescing bubbles, R is the gas constant, T is the temperature, A is the Hamaker constant, σI is the surface tension, and c is the concentration of a surfactant species. To solve this equation, sufficient initial conditions are required. The rather arbitrary initial thickness used for the film in air-water systems was set to h0 ) 1 × 10-4 m, whereas the final thickness was set at hf ) 1 × 10-8 m. However, in practice, several of the effects included in the model are very cumbersome to describe, and no accurate solution of this equation have yet been found. Several simplifying assumptions are usually adopted (e.g., no Hamaker contributions, no surface impurities, simplified interface geometry, etc.), making analytical solutions possible. A typical coalescence time was found by integration, ∆tij ) (rij3Fc/16σI)1/2 ln(h0/hf). rij denotes the equivalent bubble radius. A similar collision time, ∆tcol, was suggested by Luo70 for the flowing bubble column system, assuming deformable and fully mobile interfaces

5130


∆tcol ) 2tmax ≈ (1 + ξij)

x

Fcdi3 3(1+ ξij2)(1 + ξij3) σI (48) (Fd/Fc + CV)

where the bubble ratio is ξij ) di/dj, CV is the addedmass coefficient, and tmax is the time between the first contact and when the film area between the two colliding bubbles reaches its maximum value. The corresponding coalescence time, ∆tcoal, first suggested by Chester148 and further extended by Luo,70 yields

Fcvj ijdi2 ∆tcoal ≈ 0.5 (1 + ξij)2σI

(49)

where vj ij ) (vj i2 + vj j2)1/2 ) vj i(1 + ξij-2/3)1/2 and vj i is the mean turbulent bubble velocity. To determine the mean approach velocity of bubbles in turbulent collisions, a series of unfortunate assumptions was made by Luo and Svendsen.71 First, in accordance with earlier work on fluid particle coalescence,77,135 the colliding bubbles were assumed to take the velocity of the turbulent fluid eddies having the same size as the bubbles. Luo and Svendsen71 further assumed that the turbulent eddies in liquid flows would have approximately the same velocity as neutrally buoyant droplets in the same flow. Utilizing the experimental results obtained in an investigation on turbulent motion of neutrally buoyant droplets in stirred tanks reported by Kuboi et al.,149,150 the mean square droplet velocity was expressed as vrms2(d) ) v2(d) ) 2.0(d)2/3. The experimental data also indicated that the turbulent velocity component distributions of droplets follow a Maxwellian distribution function (a brief description of the properties of the Maxwellian state is given by Gidaspow120); thus, vj drops ) (8vrms2(d)/3π)1/2 ) x1.70 (d)1/3. It was further noticed that the value of the coefficient was sensitive to the density ratio between the continuous and dispersed phases, Fd/Fc. This approach will be discussed in further detail in the breakage section. Note, however, that Brenn et al.151 investigated unsteady bubbly flow with very low void fractions and concluded that the velocity probability density functions of bubbles and liquid are better described using two superimposed Gaussian functions. Hagesaether et al.72 derived a model for film drainage in turbulent flows and studied the predictive capabilities. In line with all other attempts to estimate these time scales, it was concluded that the present models are not sufficiently accurate and that sufficiently detailed data on bubbly flows are not available for model validation. For droplet flows, it was found that the pure drainage process (without interfacial mass-transfer fluxes) was predicted with fair accuracy, whereas no reliable coalescence criteria was found (see also Klaseboer et al.152,153). Furthermore, it was concluded that a head-on collision is not representative for all possible impact parameters. Orme154 and Havelka et al.,155 among others, noticed that the impact parameter is of great importance for the droplet-droplet collision outcome in gas flows. However, no such collision outcome maps have yet been published for bubble-bubble collisions. Saboni et al.156,157 developed drainage models for partially mobile plane films to describe film drainage and rupture during coalescence in liquid-liquid disper-

sions, taking into account the interfacial tension gradients generated by interfacial mass transfer. The resulting Marangoni forces accelerated the film drainage, which, in general, corresponds to dispersed-tocontinuous phase transfer and diminishes film drainage in the negative case. Similar effects are expected to occur for the gas-liquid systems operated in bubble columns, but no detailed experimental analysis on gasliquid dispersions has yet been reported. Models for the Macroscopic Breakage Rate, ΩB(di,dj). During the past decade, considerable attention has been focused on the macroscale modeling of bubble breakage in gas-liquid dispersions.67,77,137,145,158 A brief outline of the important milestones is given next. Fluid particle breakage controls the maximum bubble size and can be greatly influenced by the continuousphase hydrodynamics and interfacial interactions. Therefore, a generalized breakage mechanism can be expressed as a balance between external stresses, σij + σt,ij, that attempt to disrupt the bubble and the surface stress, σI/d, that resists the particle deformation. Thus, at the point of breakage, these forces must balance, σ ≈ σI/(d/2). This balance leads to the prediction of a critical Weber number, above which the fluid particle is no longer stable. It is defined by159

Wecr )

(σij + σt,ij)dmax 2σI

(50)

where dmax is the maximum stable fluid particle size and σ reflects the hydrodynamic conditions responsible for particle deformation and eventual breakage. In the case of turbulent flow, particle breakage is caused by velocity fluctuations resulting in pressure variation along the particle surface (Wecr ) Fc v′c2dmax/2σI). In laminar flow, viscous shear in the continuous phase will elongate the particle and cause breakage [Wecr ) (µc(∂vz/∂r)dmax)/2σI]. In the absence of net flow of the continuous phase such as rising bubbles in a liquid, the fluid particle breakage is caused by interfacial instabilities due to Raleigh-Taylor and Kelvin-Helmholtz instabilities.160 In most bubble column analyses, the flow is considered turbulent, and the effects of both viscous and interfacial instability are neglected basically without any further validation. The fluid particle fragmentation phenomena in a highly turbulent flow are related to the fact that the velocity in a turbulent stream varies from one point to another (i.e., validated by two-point measurements145). Therefore, different dynamic pressures will be exerted at different points on the surface of the fluid particle. Under certain conditions, this will inevitably lead to deformation and breakage of the fluid particle. According to Kocamustafaogullari and Ishii,67 the force due to the dynamic pressure can develop either through the local relative velocity around the particle, which appears because of inertial effects, or through changes in the eddy velocities over the length of the particle. Most of the published literature on bubble breakage is derived from the theories outlined by Kolmogorov161 and Hinze.159 In one phenomenological interpretation, bubble breakage occurs through bubble interactions with turbulent eddies bombarding the bubble surface. If the energy of the incoming eddy is sufficiently high to overcome the surface energy, deformation of the


surface results, which can finally lead to the formation of two or more daughter bubbles. For bubble breakage to occur, the bombarding eddies must be smaller than or equal in size to the bubble, as larger eddies only transport the bubble. To model the breakup process, the following simplifications are generally made:71 (1) The turbulence is isotropic. (2) Only binary breakage of a bubble is considered. (3) The breakage volume ratio is a stochastic variable. (4) The occurrence of breakage is determined by the energy level of the arriving eddy. (5) Only eddies of a size smaller than or equal to the bubble diameter can cause bubble breakage. The basis of the macroscopic theories is the engineering interpretation considering the velocity fluctuations as imaginary discrete entities denoted fluid slabs or eddies. Following this concept, a large number of eddies exist in the flow having a size (diameter) distribution ranging from the Kolmogorov microscale to the vessel dimensions. The basic ideas for the model development of Luo and Svendsen71 were adopted from an earlier paper by Coulaloglou and Travlarides135 considering droplet breakage in turbulent flows. In the model of Coulaloglou and Travlarides,135 the breakage density was expressed as a product of an integral or average breakage frequency (∝ 1/tb) and an integral (average) breakage efficiency determining the fraction of particles breaking. The breakage time (tb) was determined from isotropic turbulence theory, and the breakage efficiency was determined from a probable fraction of turbulent eddies colliding with droplets that have kinetic energy greater than the droplet surface energy. Prince and Blanch77 further postulated that bubble breakage is a result of collisions between particles and turbulent eddies and that the collision density can be calculated following arguments from the kinetic theory of gases. An integral (average) breakage density ΩB(di:dj), with units 1/(s m3), was thus expressed as a product of an average eddy-bubble collision density ωB(di:dj) and an average breakage efficiency pB(di:dj); thus, ΩB(di:dj) ) ωB(di:dj) pB(di:dj). No explicit breakage frequency function was given. Luo and Svendsen,71 on the other hand, considered the individual eddy-particle collisions and argued that, for a bubble to break, each of the colliding eddies must have sufficient energy to overcome the increase in bubble surface energy, as illustrated in Figure 13, and have a size on the order of the bubble diameter. The differential breakage density was then expressed as a product of an eddy-bubble collision probability density, ωB,λ(di,λ), and a breakage efficiency, pB(di:dj,λ), both of which depended on the eddy size (λ). The total breakage rate for a bubble of size di is then dmax

ΩB(di) )

∑

dj)dmin

ΩB(di:dj)

( ) 1

s m3

(51)

∫λd

min

ωTB(di,λ) pB(di:dj,λ) dλ

( ) 1 s m3

Figure 14. Sketch of a collision tube of an entering eddy moving through the tube with a velocity v j λ. The bubbles within the tube are assumed to be frozen or stationary.

assumption that only eddies of size smaller than or equal to the bubble diameter can cause bubble breakage. A bubble being dispersed in a large-scale fluid eddy (larger than the bubble scale, d) is assumed to be advected together with the eddies in the fluid. Unfortunately, the breakage frequency is found to be highly sensitive to the choice of integration limits.146,158 Note that no explicit breakage frequency was given. Luo and Svendsen71 considered the collision density of eddies with a given velocity, vj λ, bombarding a number, ni, of locally frozen bubbles, as illustrated in Figure 14

π ωTB(di,λ) ) nifλhC(di,λ) ) nifλ (di + λ)2vj λ 4

[

1 3

s m (m)

whereas the individual rate of breaking a parent bubble of size di into the daughter size classes dj is expressed as

ΩB(di:dj) )

Figure 13. Sketch of the bubble breakage surface energy balance. The mean kinetic energy of an eddy of size λ breaking a bubble of size di, ej(di,λ), is assumed to be larger than the increase of the bubble surface energy required breaking the parent bubble di into a daughter bubble dj and a second corresponding daughter bubble, es(di,dj).

(52)

and the simple eddy-bubble collision probability density, ωTB(di,λ), has units 1/[s m3 (m)]. The upper integration limit for the eddy size is based on the model

]

(53)

where fλ is the number of eddies of size between λ and λ + dλ with units of 1/[m3 s (m)] and vj λ is the turbulent velocity of eddies of size λ. Politano et al.162 defined the effective collision cross-sectional area as σA ) (π/4)[(di + λ)/2]2, in accordance with the standard definition in kinetic theory (see also Kolev130). The mean turbulent velocity of eddies with size λ in the inertial subrange of isotropic turbulence was assumed to be equal to the velocity of the neutrally

5132


buoyant droplets measured by Kuboi et al.149,150 Kuboi et al.149,150 found that the turbulent velocity of droplets could be expressed by the Maxwell distribution function, giving the mean eddy velocity of

[ ]

8v2(λ) vj λ ) 3π

1/2

(54)

where v2(λ) ) (8β˜ /2π)(λ)2/3 is in accordance with the classic turbulence theory. The coefficient value is coincidentally the same as that obtained by Kuboi et al.,149,150 (8β˜ /3π) ≈ C ≈ 2.0. However, it was commented by Kuboi et al.149,150 that the comparison of these relations cannot be very decisive, in view of the fact that there is a large difference between the processes and types of fluids used to obtain these relations. Kuboi et al.149,150 also investigated the effect of a difference in particle fluid density producing nonneutrally buoyant particle flows and concluded that the parameter value discussed above is very sensitive to the density ratio. Therefore, the application of the above relation as an approximation for the bubble velocity is highly questionable. To employ eq 53, the number probability density of eddies of a particular size must be determined. Luo and Svendsen71 assumed that the turbulence is isotropic and that the eddy size of interest lies in the inertial subrange. An expression for the number probability density of eddies as a function of wavelength for these conditions was formulated using conceptual ideas from Azbel163 (p 85) and Azbel and Athanasios.164 The turbulent energy spectrum function, E(k), can be interpreted as the kinetic energy contained within eddies of wavenumbers between k and k + dk, or equivalently, of size between λ and λ + dλ, per unit mass. A relationship between fλ and E(k) can thus be obtained by formulating an energy balance for eddies being interpreted both as discrete entities and as a wave function

[

J m (m)

Eeddies(λ) ) Espectra(λ)

3

]

(55)

Thus

fλ

[(

]

( ) [

dk 1 π 3 2 F λ vj λ ) E(k)FL(1 - Rg) 2 L6 dλ

)

]

J m3 (m) (56)

The functional form of the energy spectrum in the inertial subrange of turbulence is defined as143

E(k) ) Ck2/3k-5/3

[

m2 (m) s2

]

(57)

and the relationship between the wavenumber and the size of the eddy (wavelength) is k ) 2π/λ (dk/dλ ) -2πλ-2), vj λ2 ) β(λ)2/3, β ) 8β˜ /3π, and β˜ ) 3/5Γ(1/3)Ck ) 2.41 [Γ(1/3) ≈ 2.6789.] [The value of the latter parameter is difficult to determine as different authors give different values: Luo70 used β˜ ) 3/5Γ(1/3)Ck ) 2.41; Batchelor144 (p 123) defined β˜ ) 9/5Γ(1/3)Ck ) 7.23; Pope143 (p 232) obtained another value, β˜ ≈ 2.0; Martıńez-Bazań et al.165 referred to Batchelor166 and claimed that β˜ ≈ 8.2; Risso and Fabre145 referred to Batchelor167 (p 120) and defined β˜ ) 27/55Γ(1/3)Ck ) 1.97, which is about the

same as the value given by Pope;143 and Lasheras et al.158 reviewed several breakage models and pointed out that the value of this parameter ranges from about β˜ ) 2.045 to β˜ ) 8.2. Because of the disagreement on this value in the literature, we have chosen to proceed with the value of Luo.70) The number density of eddies, fλ, was thus defined as

[

9Ck 1 0.822(1 - Rg) fλ ) (1 - Rg) 5/3 2/3 4 ) β˜ 2 π λ λ4

]

1 m3 (m) (58)

where Rg is the volume fraction of the dispersed phase. However, it is not clear whether Luo and Svendsen71 distinguish between the theoretical parameter value for the turbulent eddy velocities and the empirical parameter for the bubble velocity in the model implementation. In this report, the “theoretical value” is adopted for the turbulent eddy properties. The collision density of eddies with sizes between λ and dλ on particles of size di can be expressed as

(di + λ)2 π ωTB(di,λ) ) (0.822)x2.045(1 - Rg)ni(λ)1/3 4 λ4 1 (59) 3 s m (m)

[

]

In dimensionless variables, this is

ωTB(ξ) ) 0.923(1 - Rg)(di)1/3ni

(1 + ξ)2 di2ξ11/3

(60)

where ξ ) λ/di. To close the collision density model, a breakage probability (efficiency) is needed. For each particular eddy hitting a particle, the probability for particle breakage was assumed to depend not only on the energy contained in the arriving eddy but also on the minimum energy required to overcome the surface area increase due to particle fragmentation. The latter quantity was determined by the number and sizes of the daughter particles formed in the breakage processes. To determine the energy contained in eddies of different scales, a distribution function of the kinetic energy for eddies in turbulent flows is required. A Maxwellian distribution function might be a natural and consistent choice,136 as the eddy velocity is assumed to follow this distribution, but Luo and Svendsen71 preferred an empirical energy distribution density function for fluid particles in liquid developed by Angelidou et al.168 The breakage probability function used yields

pB(di:dj,λ) ) 1 -

∫0e

[

]

es(di,dj) 1 exp de′s ) ej(di,λ) je(di,λ)

s,c

1-

∫0 exp(-χ) dχ′ ) exp(-χc) χc

(61)

(Notice that this efficiency function is not necessary volume-or mass-conserving, as Luo and Svendsen71 considered the breakage efficiency function to be equal to the kinetic energy distribution function. It would probably be better to consider the breakage distribution function as purely proportional to the empirical kinetic energy distribution function and determine the probability constant by requiring bubble volume or mass


conservation within the breakage process.) In eq 61, χ ) es(di,dj)/ej(di,λ) is a dimensionless energy ratio, where ej(di,λ) is the mean kinetic energy of an eddy of size λ and es(di,dj) is the increase in surface energy when a particle of diameter di is broken into two particles of size dj and (di3 - dj3)1/3. A critical dimensionless energy for breakage to occur was thus defined as χc ) [es(di,dj)/ ej(di,λ)]|c ) [-(12CfσI)/βFc2/3di5/3ξ11/3]|c. The mean kinetic energy of an eddy with size λ, ej(di,λ), was expressed as 2 π vj λ 2β 2β ) Fc λ11/32/3 ) Fc ξ11/3di11/32/3 ej(di,λ) ) Fc λ3 6 2 3 3 (62)

whereas the increase in surface energy was given by

es(di,dj) ) πσI[dj2 + (di3 - dj3)1/3 - di2] ) πσIdi2[fVij2/3 + (1 - fVij)2/3 - 1] ) Cf(fVij)πσIdi2 (63) where Cf(fVij) is defined as the increase coefficient of surface area, which depends only on the breakage volume fraction, fVij ) dj3/di3. The es(di,dj) function is symmetrical about the breakage volume fraction, fVij ) 0.5. Combining eqs 60, 61, and 52, the breakage density of one particle of size di that breaks into particles of sizes dj and (di3 - dj3)1/3 is given by

ΩB(di:dj) )

∫λd

min

ωTB(di,λ)pB(di:dj,λ)dλ

) 0.923(1 - Rg)ni

( ) di2

1/3

(

( ) 1 s m3

(64)

)

12CfσI (1 + ξ)2 exp dξ min ξ11/3 βFc2/3di5/3ξ11/3

∫ξ1

where ξmin ) λmin/d, λmin/λd ≈ 11.4 - 31.4, and λd is the Kolmogorov microscale. Recall that Cf is a function of dj, i.e., Cf ) Cf(fVij) ) Cf(dj3/di3). No explicit breakage frequency was determined in this model. A severe limitation of the original breakage density closure of Luo and Svendsen71 is that the model produces too many very small particles regardless of the numerical size resolution used; the model is thus not grid-independent (i.e., in the particle size grid). Several contributions have focused on the modification of the basic model of Luo and Svendsen71 intending to avoid this limitation.32,74,137 Hagesaether et al.74 and Wang et al.137 extended the basic breakage model of Luo and Svendsen71 by assuming that, when an eddy of size λ has kinetic energy e(λ), the daughter bubble size is limited by two constraints: One is that the dynamic pressure of the turbulent eddy 1/ F v 2 2 c j λ must be larger than the capillary pressure σI/d′′, resulting in a minimum breakage fraction fVmin (or bubble size dj,min). d′′ denotes the diameter of the smaller daughter particle (or two times the minimum radius of curvature). The second constraint is in accordance with the work of Luo and Svendsen71 that the eddy kinetic energy must be larger than the increase of the surface energy during the breakage. These constraints result in a maximum breakage fraction fVmax (or bubble size dj,max). Lehr et al.32 derived a similar breakage density function using only the capillary constraint and assumed that the interfacial force on the bubble surface

and the initial force of the colliding eddy balance each other. In contrast, Hagesaether et al.74 and Wang et al.137 assumed that, during breakage, the inertial force of the colliding eddy is often larger than the interfacial force and bubble deformation is strengthened until breakage occurs. Therefore, the force balance of Lehr et al.32 might not be satisfied during bubble breakage. Politano et al.162 studied the equilibrium between coalescence and breakage in homogeneous flows with isotropic turbulence using a population balance model. The population balance was solved using a multigroup approach. The daughter size distribution function was approximated by a uniform function, by a delta function, and by the model proposed by Luo and Svendsen.71 The breakage rate was calculated using either the model proposed by Luo and Svensen71 or the model proposed by Prince and Blance.77 Significant differences in the resulting bubble breakage rate, and therefore in the bubble size distribution, were observed comparing the model performance with experimental data on bubbly flows in an agitated tank. The model of Luo and Svendsen71 was in very good agreement with the experimental data, whereas the model of Prince and Blanch77 appears to overpredict the bubble size. In summary, most of the models discussed so far are based on eddy collision arguments that rely on the interpretation that turbulent flows consist of a collection of fluid slabs or eddies that are treated using relationships developed in the kinetic theory of gases. This engineering eddy concept is impossible to validate regarding the eddy shape, the number densities, and the breakage mechanisms discussed above. Further, the resulting breakage models require the specification of the minimum and maximum eddy sizes that are capable of causing particle breakage, as well as a criterion for the minimum bubble size dmin. Finally, the macroscopic formulation is considered less fundamental and less general than the average microscopic one (discussed in the next section). An unfortunate consequence of this modeling framework is that a discrete numerical scheme for the particle size is embedded within the source term closures; thus, optimizing the numerical solution procedure is not a trivial task. Microscopic Population Balance Formulation. The dispersed-phase system is considered as a population of particles of the dispersed phase distributed not only in physical space (i.e., in the ambient continuous phase) but also in an abstract property space.122,124 In the terminology of Hulburt and Katz,122 one refers to the spatial coordinates as external coordinates and the property coordinates as internal coordinates. The joint space of internal and external coordinates is referred to as the particle state space. The quantity of basic interest is the average number of particles per unit volume of the particle state space. The population balance is thus an equation for the number density and can be regarded as representing a number balance on particles of a particular state. The balance equation basically accounts for the various ways in which particles of a specific state can either form or disappear from the system. When particle states are continuous (i.e., in the internal coordinates as well), then processes, which cause their smooth variation with time, must contribute to the rates of formation and disappearance of specific particle types. Such processes can be viewed as convective processes because they result from convective motion in particle state space.

5134


The number of particles of a particular type can change by processes that create new particles (“birth” processes) and destroy existing particles (“death” processes). In bubble columns, birth of new particles can occur through both breakage and coalescence processes. The bubble breakage and coalescence processes also contribute to death processes, as a particle type that either breaks (into other particles) or coalesces with another particle no longer exists as such following the event. The phenomenological utility of population balance models lies in the convective processes as well as the birth and death processes. Let f1(x,r,t) be the average number of particles per unit volume of the particle state space at time t, at location x ≡ (x1, x2, ..., xd) in property space (d representing the number of different quantities associated with the particle) and r ≡ (r1, r2, r3) in physical space. This average number density function is considered a smooth function of its arguments x, r, and t. Thus, f1(x,r,t) can be differentiated as many times as desired with respect to any of its arguments. The continuous-phase variables are represented by a vector Y(r,t) ≡ (vc, vd, Rc, Rd, p, Fc, Fd, ...) that is calculated from the governing continuum transport equations and the problem-dependent boundary conditions. If dVx and dVr denote infinitesimal volumes in property space and physical space, respectively, located at (x, r), then the average number of particles in dVx dVr is given by f1(x,r,t) dVx dVr. The local (average) number density in physical space, that is, the total number of particles per unit volume of physical space, denoted N(r,t), is given by

N(r,t) )

∫V f1(x,r,t) dVx x

(65)

Whereas the rate of change of external coordinates refers to motion through physical space, that of internal coordinates refers to motion through an abstract property space. Separate velocities can then be defined for the internal coordinates, X4 (x,r,Y,t), and for the external coordinates, R4 (x,r,Y,t). It is thus possible to identify particle (number) fluxes. f1(x,r,t) R4 (x,r,Y,t) represents the particle flux through physical space, and f1(x,r,t) X4 (x,r,Y,t) is the particle flux through internal coordinate space. Let ψ(x,r) be an extensive property (i.e., the property changes with the quantity of matter present) associated with a single particle of state (x, r). The amount of this extensive property contained in all of the particles in the material population above, denoted ψ(t), is given by

ψ(t) )

∫V(t)ψ(x,r) f1(x,r,t) dV′

(66)

where V(t) is considered a combined abstract material volume (i.e., containing Vr and Vx) of a hypothetical medium containing the embedded particles. Using the preceding nomenclature, a balance of the property (ψ) for the total population of entities simultaneously found in the volume, V, can be formulated in a suitable frame of reference. Randolph124 and Randolph and Larson,125 as well as Ramkrishna,126,127 adopted the Lagrangian point of view. In accordance with the standard approach for deriving the multifluid model,100 the Eulerian point of view is used in this contribution. In this framework, the control volume (CV) is fixed with no local translational velocity (vCV ) 0).

The balance principle applied to an Eulerian control volume (CV) is, by definition, expressed as a balance of accumulation (the term on the left-hand side of the equation), net inflow by convection (the first term on the right-hand side of the equation), and volumetric production (the second term on the right-hand side of the equation). The Eulerian transport equations can thus be cast in the generalized form

(rate of ψ accumulation in the CV) ) (net inflow rate of ψ by convection across the CV surface) + (net source of ψ within the CV) This approach thus consider the balance principle of the quantity ψ within a fixed abstract volume (V) bounded by a fixed abstract surface area A

d dt

∫V(ψf1) dV′ ) -∫A(vψf1)‚n dA′ + ∫Vφ dV′

(67)

In this notation, n is an abstract outward-directed unit vector normal to the surface of an abstract control volume, d/dt is an extended total time derivative operator, v is the phase-space velocity, ψ is the balanced quantity, and φ is the net source term. The term on the left-hand side of eq 67 can be transformed using an extended Leibnitz theorem into an integral over only volume because the contributions from the fixed abstract surfaces (A) vanish as the surface phase-space velocities become zero. In terms of the previously defined nomenclature, Randolph and Larson125 expressed a generalized theorem in the form

d dt

∫V(t)(ψf1) dV′ ) ∫V(t)

[

]

∂(ψf1) + ∇‚(vCVψf1) dV′ ∂t

(68)

where vCV denotes a surface phase-space velocity vector (≡ X4 + R4 ). Adopting an Eulerian frame instead, the theorem simplifies considerably

d dt

∂(ψf1) dV′ ∂t

∫V(ψf1) dV′ ) ∫V

(69)

The convective transport terms in eq 67 can be rewritten as a volume integral using an extended version of Gauss’ theorem

∫A(vψf1)‚n dA′ ) ∫V∇‚(vψf1) dV′

(70)

Equation 67 can then be rewritten as a volume integral

(

∫V

)

∂(ψf1) + ∇‚(vψf1) - φ dV′ ) 0 ∂t

(71)

Equation 71 must be satisfied for any V, so the expressions inside the volume integral must be equal to zero. The local instantaneous equation for a general balanced quantity ψ then yields (eq 72)

∂(ψf1) + ∇‚(vψf1) - φ ) 0 ∂t

(72)

This relation can be referred to as a general transport equation.


The corresponding source term φ contains the birth and death terms representing the net rate of production of particles of state (x, r) at time t in an environment of state Y, thus φ(x,r,Y,t). Considering the bubble coalescence and breakage phenomena only, the function φ(x,r,Y,t) can be expressed as φ(x,r,Y,t) ) BC(x,r,Y,t) DC(x,r,Y,t) + BB(x,r,Y,t) - DB(x,r,Y,t). With ψ(x,r) ) 1, the generalized population balance framework becomes

∂f(x,r,t) + ∇r‚(f(x,r,t)vr) + ∇x‚(f(x,r,t)X4 ) ∂t ) BB(x,r,Y,t) - DB(x,r,Y,t) + BC(x,r,Y,t) -

(73)

DC(x,r,Y,t) where f(x,r,t) is the bubble number probability density [1/m3 (x)], which varies with internal coordinates (x), spatial position (r), and time (t). The first term on the LHS of eq 73 is the change in the bubble number probability density with time. The second term is the change in the bubble number probability density due to convection. The third term is the change in the bubble number probability density due to convection in the internal coordinates denoting several particle growth phenomena. The RHS is the source terms. The source terms represent the change in the number probability density, f1, due to particle breakup and coalescence. BB and DB represent the birth and death, respectively, of bubbles due to breakage, and BC and DC, respectively, are the birth and death of bubbles due to coalescence. To close this model formulation, constitutive relations are needed for both the growth and source term functions. This process is of extreme importance as it represents the weakest part of population balance modeling. Possibly, one could expect that these terms depend on the number density function (kinetic theory) at the instant in question, but the detailed nature of these relationships are generally unknown and require further attention. To some degree, the complexity of the growth terms is also dependent on the choice of internal coordinates or particle properties used to characterize the dispersed phase, as mentioned earlier. In most reports, the particle volume (diameter) is used, so the gas expansion due to pressure, temperature, and composition changes as well as interfacial mass-transfer fluxes has to be incorporated through this term. In contrast, using the particle mass as the internal coordinate, only the interfacial mass-transfer fluxes have to be incorporated through this term. However, it is still not feasible to solve the microscopic population balance equations derived above by direct numerical simulations, so some kind of averaging is required to enable solution with reasonable time and physical space resolutions. In accordance with previous work,8 the microscopic transport equation is averaged over a physical averaging volume and thereafter over time. The result is

∂〈f(x,t)〉 + ∇r‚(〈f(x,t)〉〈vr〉) + ∇x‚(〈f(x,t)〉〈X4 〉) ∂t ) 〈BB(x,Y,t)〉 - 〈DB(x,Y,t)〉 + 〈BC(x,Y,t)〉 - (74) 〈DC(x,Y,t)〉

One can now introduce number-density-weighted average velocities, i.e.

〈f(x,t)〉〈vr〉 ) 〈f(x,t)〉 〈vr〉f

(75)

〈f(x,t)〉〈X4 〉 ) 〈f(x,t)〉 〈X4 〉f

(76)

and

This makes the velocities in the population balance different from the mass- or phase-weighted average velocities obtained solving the multifluid model. According to Dorao and Jakobsen,134 this discrepancy is an argument for the formulation of a mass density population balance instead of a number balance to achieve a consistently integrated population balance within the multifluid modeling framework. One can also perform standard Reynolds decomposition of the variables and then time average the equation. The typical turbulence modeling of the resulting covariance terms gives rise to diffusive terms in the balance equation.121 The physics involved in these terms is not understood, however, making this modeling issue a challenging task

〈f(x,t)〉〈vr〉 ≈ 〈f(x,t)〉 〈vr〉 + 〈f(x,t)〉′〈vr〉′

(77)

〈f(x,t)〉〈X4 〉 ≈ 〈f(x,t)〉 〈X4 〉 + 〈f(x,t)〉′〈X4 〉′

(78)

and

Provided that sufficient relations for the contact area are available, the bubble growth term due to interfacial mass transfer can be modeled in accordance with the well-known two-film theory and the ideal gas law.32 The modeling of the source and sink terms due to bubble breakage and coalescence are more difficult from a physical point of view, and the existing theory is rather complex and not easily available. In the following subsections, these source terms are thus elucidated in further detail. Microscopic Source Term Closures. In any twophase flow field, the initial bubble size is determined in terms of the mechanisms of fluid particle generation such as formation of bubbles at an orifice, perforated or sintered distributor plates, or other sparger devices. However, in most dispersions, the initial fluid particle size might be too large or too small to be stable. In these cases, the fluid particle size is further determined by either the breakage or coalescence mechanism, respectively. When a fluid particle exceeds a critical value, the particle interface becomes unstable, and breakage is likely to occur. Similarly, when fluid particles are smaller than a certain critical dimension, coalescence is likely to occur through a series of collision events. It has thus been common to relate particle breakage to a maximum attainable size of the particle and particle coalescence to a minimum size. For determining the maximum and minimum sizes of fluid particles, several criteria for breakage and coalescence have been reported in the literature. However, these criteria do not treat a distribution of fluid particle sizes, but describe the particle size limits of breakage and coalescence. In the following discussion, a brief description of the physical mechanisms that have been considered for formulating closure laws for the more elaborated source

5136


and sink terms within the population balance equation is provided. Emphasis is placed on a selected number of the novel closures suitable for bubble column modeling. A fairly general framework has also been formulated for the source terms considering bubble breakage and coalescence.67,121,122,127,169 However, no standard notation has yet emerged in the literature, so to avoid unnecessary confusion, the following subsections are presented in accordance with the textbook of Ramkrishna.127 Considering the breakage terms, the breakage behaviors of different particles are assumed to be independent of each other. The average loss rate of particles of state (x, r) per unit time by breakage yields

DB(x,r,Y,t) ) bB(x,r,Y,t) f1(x,r,t)

(79)

where bB(x,r,Y,t) has the dimension of reciprocal time and is often called the breakage frequency. It represents the fraction of particles of state (x, r) breaking per unit time. The breakage processes are often considered to be random in nature, so the modeling work usually adopts probabilistic theory, as will be outlined later. The average production rate for particles of state (x, r) originating from breakage of particles of all other particle states, considering both internal and external coordinates is given by

BB(x,r,Y,t) )

∫V ∫V ν(x′,r′,Y,t) bB(x′,r′,Y,t) r

x

PB(x,r|x′,r′,Y,t) f1(x′,r′,t) dV′r dV′x (80) where ν(x′,r′,Y,t) denotes the average number of particles formed from the breakage of a single particle of state (x′, r′) in an environment of Y at time t. PB(x,r|x′,r′,Y,t) denotes the probability density function for particles from the breakage of a particle of state (x′, r′) in an environment of state Y at time t that have state (x, r). In the engineering literature, PB is commonly referred to as the daughter size distribution function {i.e., a probability density with units of 1/[m3 (x)] and not a probability that is dimensionless} denoting the size distribution of daughter particles produced upon breakage of a parent particle. The integrand on the RHS represents the rate of formation of particles of state (x, r) formed by breakage of particles of state (x′, r′). That is, the number of particles of state (x′, r′) that breaks per unit time is bB(x′,r′,Y,t) f1(x′,r′,t) dV′r dV′x. Multiplying the above expression by ν(x′,r′,Y,t) yields the number of new particles resulting from the breakage processes, ν(x′,r′,Y,t) bB(x′,r′,Y,t) f1(x′,r′,t) dV′r dV′x, of which a fraction PB(x,r|x′,r′,Y,t) dV′r dV′x represents particles of state (x, r). Usually, binary breakage is assumed, for which ν(x′,r′,Y,t) ) ν ) 2. However, experimental determination of this variable is highly recommended. The function PB(x,r|x′,r′,Y,t) should also be determined from experimental observations or, if future understanding allows, by detailed modeling of the breakage process. PB(x,r|x′,r′,Y,t) is a probability density function and should satisfy the normalization condition (in the internal coordinates)

∫V PB(x,r|x′,r′,Y,t) dV′x ) 1 x

(81)

This function also needs to fulfill other obvious conservation law properties that constrain the breakage process.127 The coalescence functions are considered next. The population balance source term due to coalescence is usually defined as127

BC(x,r,Y,t) )

∫V ∫V δ1aC(x˜ ,r˜ ;x′,r′,Y) r

x

|

∂(x˜ ,r˜ ) f2(x˜ ,r˜ ;x′,r′,t) ∂(x,r)

(x′,r′))

dV′x dV′r (82)

where δ represents the number of times identical pairs have been considered in the interval of integration, so 1/δ corrects for the redundancy. aC(x˜ ,r˜ ;x′,r′,Y,t) denotes the coalescence frequency or the fraction of particle pairs of states (x˜ , r˜ ) and (x′, r′) that coalesce per unit time. The coalescence frequency is defined for an ordered pair of particles, although from a physical viewpoint, the ordering of particle pairs should not alter the value of the frequency. In other words, aC(x˜ ,r˜ ;x′,r′,Y,t) satisfies a symmetry property: aC(x˜ ,r˜ ;x′,r′,Y,t) ) aC(x′,r′;x˜ ,r˜ ,Y,t). It is thus essential to consider only one of the above orders for a given pair of particles. To proceed, it is necessary to assume that it is possible to solve for the particle state of one of the coalescing pair given those of the other coalescing particle and the new particle (as the three variables are not independent). Thus, given the state (x, r) of the new particle and the state (x′, r′) of one of the two coalescing particles, the states of the other coalescing particle is known and denoted by [x˜ (x,r|x′,r′),r˜ (x,r|x′,r′)]. As in the classical kinetic theory of gases, the pair density function f2 is impossible to determine analytically, and some closure approximation has to be made. In the population balance derivation, one follows a quite common procedure in kinetic theory and makes the coarse approximation f2(x˜ ,r˜ ;x′,r′,t) ≈ f1(x˜ ,r˜ ,t) f1(x′,r′,t). This assumption implies that there is no statistical correlation between particles of states (x′, r′) and (x˜ , r˜ ) at any instant t. The sink term due to coalescence is usually defined as

DC(x,r,Y,t) )

∫V ∫V aC(x′,r′;x,r,Y) r

x

f2(x′,r′;x,r,t) dV′x dV′r

[

]

1 (83) s m3 (x)

Although Carrica et al.13 and Hagesaether146 suggested that the most convenient choice of internal coordinate for the bubble column case is the particle mass, m, for pedagogic and practical reasons, we use the bubble diameter, d, as the only internal coordinate in this review. The expressions for the continuous source term functions thus yield

BB(d,r,Y,t) )

∫V ∫d∞ν(d′,r′,Y) bB(d′,r′,Y) r

PB(d,r|d′,r′, Y) f1(d′,r′,t) dd′ dV′r DB(d,r,Y,t) ) bB(d,r,Y) f1(d,r,t)


BC(d,r,Y,t) )

∫V ∫0daC(d˜ ,r˜ ;d′,r′,Y) f1(d˜ ,r˜ ,t)

1 2

1 Vr

r

∫V bB(d,r,Y) f1(d,r,t) dVr ≈ r

f1(d′,r′,t) dd′ dV′r DC(d,r,Y,t) ) f1(d,r,t)

∫V′r∫0

∞

1 Vr

aC(d,r,d′,r′,Y)

〈BB(d,Y,t)〉 ) )

1 Vr

r

f1(d′,r′,t) dV′r dd′

Vr

r

r

∫V BC(d,r,Y,t) dVr 1 Vr

d

a (d ˜ ,r˜ ;d′,r′,Y) f1(d ˜ ,r˜ ,t) V C

r

∫V ∫V aC(d,r;d′,r′,Y) dV′r dVr r

r

Vr

∫V DC(d,r,Y,t) dVr

1 Vr

∫0∞V1r∫V ∫V f1(d,r,t) aC(d,r,d′,r′,Y) r

f1(d′,r′,t) dV′r dVr dd′ To proceed Ramkrishna,127 (pp 59 and 74) introduced the usual assumptions in volume averaging,170 neglecting all correlation terms occurring within the source term expressions, meaning that the averages of products are approximated as products of averages. The following approximations are introduced

∫V ν(d′,r′,Y) bB(d′,r′,Y)

∫V PB(d,r|d′,r′,Y) dVr r

Vr

r

×

f1(d′,r′,t) dV′r

∫V ν(d′,r′,Y) dV′r ∫V bB(d′,r′,Y) dV′r r

1 Vr

∫V ν(d′,r′,Y) dV′r

1 Vr

∫V bB(d′,r′,Y) dV′r

(1s)

r

1 (86) (m)

( )

m3 (87) s

∫V ∫V aC(d,r;d′,r′,Y) dVr dV′r r

r

1 Vr

〈 f1(d′,t)〉 )

∫V f1(d′,r′,t) dV′r r

(85)

[ ]

〈〈aC(d;d′,Y)〉〉 ) 1 Vr

(84)

r

〈〈PB(d|d′,Y)〉〉 ) 1 P (d,r|d′,r′,Y) dVr dV′r Vr Vr Vr B

[

1 m3 (m)

]

(88)

The approximate volume-average source terms can thus be expressed as

∫d∞〈ν(d′,Y)〉〈bB(d′,Y)〉〈〈PB(d|d′,Y)〉〉〈f1(d′,t)〉 dd′ 〈DB(d,Y,t)〉 ) 〈bB(d,Y)〉〈f1(d,t)〉

∫0d〈〈aC(d˜ ;d′,Y)〉〉〈f1(d˜ ,t)〉〈f1(d′,t)〉 dd′

〈BC(d,Y,t)〉 )

1 2

〈DC(d,Y,t)〉 )

∫0∞〈f1(d,t)〉〈〈aC(d,d′,Y)〉〉〈f1(d′,t)〉 dd′

Performing turbulence modeling on these terms is a very difficult task. One usually assumes that the physical time scales involved are much longer than the time scales of the turbulent fluctuations so that no correlation functions need to be considered. The time-after-volume-averaged source terms are listed below, dropping the averaging symbols for simplicity in the proceeding model derivation

BB(d,Y,t) )

Vr

∫V ∫V PB(d,r|d′,r′,Y) dVr dV′r ∫V f1(d′,r′,t) dV′r r

r

where

r

Vr

∫V f1(d′,r′,t) dV′r.

〈BB(d,Y,t)〉 )

r

r

1 Vr

r

f1(d′,r′,t) dVr dV′r dd′

r

dVr

r

∫0 ∫V ∫

≈

r

∫∫

∫V DB(d,r,Y,t) dVr ) V1r∫V bB

1 Vr

1 ) 2

)

∫V f1(d,r,t)

〈bB(d′,Y)〉 )

(d,r,Y) f1(d,r,t) dVr

〈DC(d,Y,t)〉 )

1 Vr

r

bB(d′,r′,Y)

r

r

〈ν(d′,Y)〉 )

∫d∞∫V ν(d′,r′,Y) ∫V PB(d,r|d′,r′,Y) dVr

1 〈BC(d,Y,t)〉 ) Vr

≈

r

∫V BB(d,r,Y,t) dVr r

〈DB(d,Y,t)〉 )

r

∫V ∫V f1(d,r,t)aC(d,r,d′,r′,Y) f1(d′,r′,t) dV′r dVr

1 Vr

f1(d′,r′,t) dd′ dV′r where d ˜ ) (d3 - d′3)1/3. All of the source terms have the common units 1/[s m3 (m)]. The local source term formulation given above is very seldom used in practical simulations; some kind of averaging is required. In this work, the time-aftervolume averaging procedure is employed. Integrating the source terms over an averaging volume, Vr, yields

∫V bB(d,r,Y) dVrV1r∫V f1(d,r,t) dVr

∫d∞ν(d′,Y) bB(d′,Y) PB(d|d′,Y) f1(d′,t) dd′

(89)

r

Vr

Vr

DB(d,Y,t) ) bB(d,Y) f1(d,t)

(90)

5138


BC(d,Y,t) )

∫0daC(d˜ ;d′,Y) f1(d˜ ,t) f1(d′,t) dd′

1 2

DC(d,Y,t) )

∫0∞f1(d,t) aC(d,d′,Y) f1(d′,t) dd′

(91) (92)

All of the source terms still have the common units, 1/[s m3 (m)]. Comparing the source term expressions, eqs 89-92, with eqs 29-32, it is clearly seen that the two formulations give rise to identical expressions for the source terms only under certain conditions. The macroscopic formulation is explicitly expressed in terms of a discrete discretization scheme and is very difficult to convert to other schemes. Coalescence Frequency Closures, aC(d; d′, Y). The literature contains several models and experimental analyses on bubble coalescence.67,77,130,140,147,171,172 Intuitively, bubble coalescence is related to bubble collisions. The collisions are caused by the existence of spatial velocity differences among the particles themselves. However, not all collisions necessarily lead to coalescence. Thus, the modeling of bubble coalescence on these scales means the modeling of bubble collision and coalescence probability or efficiency mechanisms. The pioneering work on coalescence of particles to form successively larger particles was carried out by Smoluchowski.173,174 As for the collision density in the macroscopic model formulation, the average collision frequency of fluid particles is usually described assuming that the mechanisms of collision are analogous to those of collisions between molecules as in the kinetic theory of gases.140 Furthermore, as pointed out before, the units of the volume-average coalescence efficiency, aC(d;d′,Y), correspond to the units of the effective swept volume rate, hC(d;d′,Y), in the kinetic theory of gases.138 Therefore, rather than the actual collision frequency, the volume swept by a moving bubble of diameter d is calculated, from which the number of other bubbles with diameter d′ that are hit by this moving bubble of diameter d can be indirectly determined. The volume-average coalescence frequency, aC(d;d′,Y), can thus be defined as the product of an effective swept volume rate, hC(d;d′,Y), and the coalescence probability, pC(d;d′,Y) (Prince and Blanch;77 Kolev,130 p 169; Coulaloglou and Tavlarides;135 Venneker et al.;138 Tsouris and Tavlarides169)

aC(d;d′,Y) ) hC(d;d′,Y) pC(d;d′,Y)

( ) m3 s

(93)

which expresses the fact that not all collisions lead to coalescence. Modeling coalescence thus means finding adequate physical expressions for hC(d;d′,Y) and pC(d;d′,Y). Kamp et al.,140 among others, suggested that microscopic closures can be formulated in line with the macroscopic population balance approach; thus, we can define

( ) m3 s

π hC(d;d′,Y) ) (d + d′)2(vj t,d2 + vj t,d′2) 4

(

)

∆tcoal ∆tcol

BC(d,Y,t) ) 1 d h (d ˜ ;d′,Y) pC(d ˜ ;d′,Y) f1(d ˜ ,t) f1(d′,t) dd′ (96) 2 0 C

∫

DC(d,Y,t) ) f1(d,t)

∫0∞hC(d;d′,Y) pC(d;d′,Y) f1(d˜ ,t) dd′

(97)

All of the terms in the population balance equation thus have common units, 1/[m3 s (m)], noting that the collision densities for the average microscopic and the macroscopic frameworks have different units. However, by use of a discrete numerical scheme for the solution of the average microscopic model, the two formulations become very similar. Models for the Breakage Frequency, bB(d). Formulating average microscopic source term closure laws for the breakage processes usually means deriving relations for the time-after-volume-average breakage frequency, bB(d); the average number of daughter bubbles produced, ν; and the average daughter particle size distribution, PB(d|d′). During the past several decades, a few models have been developed for the bubble breakage frequency, bB(d), under turbulent conditions.137,158 Several different categories of breakage frequency models are distinguished in the literature; models based on reaction kinetic concepts,175 phenomenological models based on the turbulent nature of the system,169 and models based on purely kinematic ideas.165,176 In an interesting attempt to overcome the limitations found in the turbulent breakage models described above, Martıńez-Bazań et al.165 (see also Lasheras et al.158) proposed an alternative model in the kinetic theory (microscopic) framework based on purely kinematic ideas to avoid the use of the incomplete turbulent eddy concept and the macroscopic model formulation. The basic premise of the model of Martıńez-Bazań et al.165 is that, for a particle to break, its surface has to be deformed and, further, this deformation energy must be provided by the turbulent stresses produced by the surrounding fluid. The minimum energy needed to deform a particle of size d is its surface energy

Es(d) ) πσId2

(J)

(98)

The surface restoring pressure is Es per unit volume

(94)

and

pC(d;d′,Y) ) exp -

where suitable continuous closures are needed for the model variables such as the mean turbulent bubble approach velocities and the time scales. Multifluid models can provide improved estimates for the mean turbulent bubble approach velocities as individual velocity fields for each of the bubble phases are calculated. However, as discussed earlier, the existing closures for the time scales are very crude and need further consideration. The coalescence terms in the average microscopic population balance can then be expressed in terms of the local effective swept volume rate and the coalescence probability variables as

(95)

6Es σI σs(d) ) 3 ) 6 d πd

(Pa)

(99)

The size of these particles is assumed to be within the inertial subrange of turbulence, so the average deformation stress, which results from velocity fluctuations


existing in the liquid between two points separated by a distance d, was estimated as

1 σt(d) ) Fc∆v2(d) 2

(Pa)

(100)

where ∆v2(d) is the mean value of the velocity fluctuations between two points separated by a characteristic distance d and Fc is the density of the continuous phase. When the turbulent stresses are equal to the confining stresses, σt(d) ) σs(d), a critical diameter, dc, is defined such that particles with d < dc are stable and will not break. A particle with d > dc has a surface energy smaller than the deformation energy, σs(d) < σt(d), and thus, such a particle deforms and might eventually break up in a time tb. Martıńez-Bazań et al.165 applied Kolmogorov’s universal theory valid for homogeneous and isotropic turbulent flow to estimate the mean value of the velocity fluctuations, ∆v2(d) ) β(d)2/3. The critical diameter, dc ) (12σI/βF)3/5-2/5, is defined by the crossing point of the two curves determined by σs(d) and σt(d). Martıńez-Bazań et al.165 postulated, in accordance with Newton’s law, that the acceleration of the particle interface during deformation is proportional to the difference between the deformation and confinement forces acting on it. In other words, the probability of breaking a particle of size d in time tb increases as the difference between the pressure produced by the turbulent fluctuations on the surface of the particle, 1 /2Fc∆v2(d), and the restoring pressure caused by surface tension, 6σI/d, increase. On the other hand, the breakage frequency should decrease to a zero limiting value as this difference of pressures vanishes. Thus, the particle breakage time can be estimated as

tb ∝

d ) vbreakage

d

x

σI ∆v2(d) - 12 Fc d

(101)

where vbreakage is the characteristic velocity of the particle breakage process. The breakage frequency bB(d,) is given by

x∆v (d) - 12Fσd I

2

1 bB(d,) ) ) Kg tb

c

d

)

x

Kg

σI β(d)2/3 - 12 Fc d (102) d

where the values β ) 8.2 and Kg ) 0.25 were found experimentally for bubbly flows. The breakage frequency is zero for particles of size d e dc, and it increases rapidly for particles larger than the critical diameter, d > dc. After reaching a maximum at dbB,max ) (9/4)3/5dc ≈ 1.63dc, the breakage frequency decreases monotonically with the particle size. The maximum breakage frequency, achieved at dbB,max, is given by bB,max() ) 0.538Kgβ1/23/5{12[σI/(Fcβ)]}-2/5. Although this approach avoids the eddy concept, it is still restricted to homogeneous and isotropic turbulent flows, it contains several hypotheses that are not yet verified for bubble column flows, and it contains a few additional adjustable parameters that need to be fitted to many sets of experimental data. Furthermore, this

model concept has only been applied to systems (i.e., turbulent water jets165,176-178) where the turbulent dissipation rates are 2-3 orders of magnitude larger than what is observed in bubble columns. The bubble sizes considered were also very small, up to about 1 mm only. On the other hand, after careful investigations of the turbulent breakage processes in bubble columns, it might be possible to extend the application of this approach to bubble column systems. The experimental techniques developed for investigating the breakage of bubbles on the centerline within the fully developed region of turbulent water jets might initiate ideas on how to study these phenomena in bubble columns.177,178 The digital image-processing technique177,178 used in these studies gives satisfactory results in very dilute two-phase flows, but the method is still questionable when the bubble concentration increases considerably as for bubble columns operated in the heterogeneous flow regime. A severe drawback for integrated multifluid/population balance models is that all of the kernels suggested in the literature are very sensitive to the turbulent energy dissipation rate (). As mentioned earlier, the local variable is difficult to calculate from the k- turbulence model because the equation for the dissipation rate merely represents a fit of a turbulent length scale to single-phase pipe flow data. Therefore, further work is highly needed to elucidate the mechanisms of bubble breakage in turbulent flows. However, an alternative to the eddy concept has been reported that might make the work on model validation easier. Perhaps the greatest advantages lies in the fact that this closure is formulated within the average microscopic modeling framework, thereby avoiding the limitations of the macroscale formulation. At this point, care should be taken as the discrete macroscopic breakage rate models [ΩB(di)] are often erroneously assumed to be equal to the average microscopic breakage frequency (bB(b)). Number of Daughter Bubbles Produced, ν. In the average microscopic formulations, this parameter determines the average number of daughter particles produced by breakage of a parent particle of size d. The breakage of parent bubbles into two daughter bubbles is assumed in most investigations reported (ν ) 2). In a recent paper, Risso and Fabre145 observed that the number of fragments depended on the specific shape of the parent bubble during the deformation process and varied in a wide range. In their experiments, a videoprocessing technique was used for studies of two groups of bubbles, one group being of sizes in the range between 2 and 6 mm and the other of sizes in the range 12.421.4 mm. Two fragments were obtained in 48% of the cases, between 3 and 10 in 37% of the cases, and more than 10 fragments in 15% of the cases. The breakage process is certainly not purely binary. Implementing this effect instead of binary breakage is expected to significantly alter the number of smaller bubbles and the interfacial area concentration predicted by the population balance equation. The experimental data of Risso and Fabre145 also indicate that an equal-size daughter distribution is more common for bubble breakage than an unequal one. In contrast, Hesketh et al.179,180 investigated bubble breakage in turbulent flows in horizontal pipes and concluded that an unequal-size daughter distribution is more probable than an equal-size one. The daughter particle

5140


size distribution model of Luo and Svendsen71 relies on the assumption that unequal-size breakage is more probable than equal-size breakage in accordance with the observations of Hesketh et al., as will be discussed in the next subsection. The deviating experimental observations indicate that further experimental investigations are needed to determine reliable correlations for the physical systems and flow conditions for which the equal and unequal breakage outcomes occur. Daughter Particle Distribution, PB(d|d′). Among the papers adopting the average microscopic formulation, two of the earliest models for the daughter particle size distribution were proposed by Valentas et al.181,182 The daughter particle probability distribution functions suggested by Valentas et al. were purely statistical relations. The first was a discrete model, in which a parent particle of diameter d was assumed to split into equally sized daughter particles of diameter d/ν, where ν is the number of daughter particles formed. The second model proposed by Valentas et al. was a logical, continuous analogue of the discrete daughter particle particle distribution function (PDF). In this case, it was assumed that the daughter particle sizes were normally distributed about a mean value. Ross and Curl,175 Coulaloglou and Travlarides,135 and Prince and Blanch77 argued that the particle breakage frequency should be a function of the difference in surface energy between a parent particle and daughter particles produced and the kinetic energy of a colliding eddy. Although Ross and Curl,175 Coulaloglou and Travlarides,135 and Prince and Blanch77 developed more physical models for the particle breakage frequency, the daughter particle size distributions still relied upon statistical relations. Among the most widely used phenomenological models based on surface energy considerations is the one proposed by Tsouris and Tavlarides.169 Their expression for the daughter particle PDF yields

PB(d|d′) )

emin + [emax - e(d′)]

∫d

d min

(m1 ) (103)

emin + [emax - e(d′)] dd′

Considering binary breakage only, Tsouris and Tavlarides169 postulated that the probability of formation of a daughter particle of size d′ is inversely proportional to the energy required to split a parent particle of size d into a particle of size d′ and its complementary particle of size d ˜ ) (d3 - d′3)1/3. This energy requirement is proportional to the excess surface area generated by splitting the parent particle

[(d′d ) + [1 - (d′d ) ] - 1] (104) 2

˜ 2 - πσId2 ) πσId2 e(d′) ) πσId′2 + πσId

3 2/3

To avoid the singularity present in their daughter size distribution model for d′ ) 0 (i.e., the parent particle does not break), a minimum particle size was defined, dmin. The excess surface area relation reaches a minimum value when a particle of minimum diameter, dmin, and a complementary one of maximum size, dmax ) (d3 - dmin3)1/3, are formed. The minimum energy, emin, is given by

[( ) [ ( ) ]

e(dmin) ) πσId2

dmin d

2

+ 1-

dmin d

3 2/3

-1

]

(105)

This expression gives a minimum probability for the formation of two particles of the same size and a maximum probability for the formation of a pair made up of a very large particle and a complementary very small one. The maximum energy corresponding to the formation of two particles of equal volumes or of diameters d′ ) (d3 - d′3)1/3 ) d /21/3 is emax ) πσId2[21/3 - 1]. Luo and Svendsen71 derived a discrete expression for the breakage density of a particle of diameter di into two daughter particles of size dj and (di3 - dj3)1/3 using energy arguments similar to those employed by Tsouris and Tavlarides.169 The significant difference between the Luo and Svendsen model and its predecessors is that the former gives both a partial breakage rate, that is, the breakage rate for a particle of size di splitting into a particle of size dj and its complementary bubble, and an overall breakage rate. The previous surface energy models provided only an overall breakage rate. An expression for the daughter particle size distribution function can thus be calculated by normalizing the partial breakage rate by the overall breakage rate. As mentioned earlier, this variable was not required for the population balance closure but is rather a spin-off from the primary closures. The Luo and Svendsen daughter particle size distribution is thus determined from the expression

PB(di|dj) )

ΩB(di:dj) ΩB(di)

[-]

(106)

Note that the units of the discrete daughter bubble size distribution variable are different from the units obtained by deriving the continuous daughter size distribution function from the microscopic formulations (1/m). It is thus not trivial how the model of Luo70 should be adopted within a more fundamental modeling approach. However, similarly to the model of Tsouris and Tavlarides,169 the model of Luo and Svendsen predicted that the probability of breaking a parent particle into a very small particle and a complementary large particle is larger than the probability of equal-size breakage (Figure 15). The distribution has a U-shape, with a minimum probability for the formation of two equally sized daughter particles and a maximum probability for the formation of a very large daughter particle and its complement dmin. Although Hinze159 and Risso and Fabre,145 among others, showed and discussed the diversity of shapes of the bubbles that can be found in turbulent flows, at present, sufficient experimental data do not exist in the literature to assist in developing adequate daughter bubble probability density functions valid for bubble columns. In an attempt to overcome the limitations found in the models for the daughter particle size distribution function described above based on the eddy collision concept, Martıńez-Bazań et al.176 proposed an alternative statistical model based on energy balance principles. The model was originally intended for the prediction of air bubble breakage at the centerline of a high-Reynolds-number, turbulent water jet. It was


PB(d/d′) ∝

[21F β(d′)

2/3

- 6σI/d

c

]{21F β[(d 3

c

}

d′3)1/3]2/3 - 6σI/d (107) Relating d′ and (d3 - d′3)1/3 through the mass balance given above yields

PB(d/d′) ∝

[21F β(d) ] [(d′/d) 2/3

2

c

2/3

- (dc/d)5/3]{[1 -

(d′/d )3]2/9 - (dc/d)5/3} (108) or

PB(d*) ∝

Figure 15. Sketch of the breakage kernel of Luo and Svendsen.71 An unequal daughter size distribution is predicted by this model. In this case, the parent diameter size is, di ) 0.006 m, and the turbulent energy dissipation rate is set at ) 1 m2/s3.

assumed that, when an air bubble is injected into the turbulent water jet, the velocity fluctuations of the underlying turbulence result in pressure deformation forces acting on the bubble’s surface that, when greater than the confinement forces due to surface tension, will cause breakage. The model assumes that, when a parent particle breaks, two daughter particles are formed (ν ) 2) with diameters d′ and (d3 - d′3)1/3. The validity of this assumption is supported by the high-speed video images, presented by Martıńez-Bazań et al.176 The diameters are related through the conservation of mass, (d3 - d′3)1/3 ) d[1 - (d′/d)3]1/3. The particle-splitting process was not considered purely random, as the pressure fluctuations and thus the deformation stress (σt) in homogeneous and isotropic turbulence is not uniformly distributed over all scales. It was further assumed that there is a minimum distance, dmin, over which the turbulent stresses acting between two points separated by this distance, 1/2Fcβ(dmin)2/3, are just equal to the confinement pressure due to surface tension: σt(dmin) ) σs(d). In other words, at this distance, the turbulent pressure fluctuations are exactly equal to the confinement forces for a parent particle of size d. The probability of breaking off a daughter particle with size d′ < d′min ) (12σI/βFcd)3/2-1 should therefore be zero. The fundamental postulate of the Martıńez-Bazań et al.176 model is that the probability of splitting off a daughter particle of any size such that d′min < d′ < d is proportional to the difference between the turbulent stresses over a length d′ and the confinement forces holding the parent particle of size d together. For the formation of a daughter particle of size d′, the difference in stresses is given by ∆σt,d′ ) 1/2Fcβ(d′)2/3 - 6σI/d. For each daughter particle of size d′, a complementary daughter particle of size (d3 - d′3)1/3 is formed with a difference of stresses given by ∆σt,(d3-d′3)1/3 ) 1/2Fcβ[(d3 - d′3)1/3]2/3 - 6σI/d. The model states that the probability of forming a pair of complementary daughter particles of sizes d′ and (d3 - d′3)1/3 from the splitting of a parent particle of size d is related to the product of the excess stresses associated with the length scales corresponding to each daughter particle size. That is

[21F β(d) ] (d 2/3

c

2

2/3

- Λ5/3)[(1 - d3)2/9 - Λ5/3] (109)

where Λ ) dc/d and dc is the critical diameter defined as dc ) (12σI/βFc)3/5-2/5. The critical parent particle diameter defines the minimum particle size for a given dissipation rate of turbulent kinetic energy for which breakage can occur. The minimum daughter diameter defines the distance over which the turbulent normal stresses just balance the confinement forces of a parent particle of size d. The minimum diameter, therefore, gives the minimum length over which the underlying turbulence can pinch off a piece of the parent particle. This length is not arbitrarily selected; rather, its determination is based on kinematics. This model further assumes that the size of the parent particle is in the inertial subrange of turbulence. Therefore, it implies that dmin e d′ e dmax provided that dmin > λd, where λd is the Kolmogorov length scale of the underlying turbulence. Otherwise, dmin is taken to be equal to λd. However, no assumption needs to be made about the minimum and maximum eddy size that can cause particle breakage. All eddies with sizes between the Kolmogorov scale and the integral scale are taken into account. The daughter particle probability density function can be obtained from the probability expression given above * max P(d*) by utilizing the normalization condition: ∫ddmin * d(d*) ) 1. The PDF of the ratio of diameters d* ) d′/d, P/B(d*), can then be written as

P/B(d*) )

(d*2/3 - Λ5/3)[(1 - d*3)2/9 - Λ*5/3]

∫dd

/ max / min

(d*2/3 - Λ5/3)[(1 - d*3)2/9 - Λ*5/3]d(d*) (110)

Note that PB(d′|d) ) P/B(d*)/d, and Λ ) dc/d ) (dmin/ d)2/5. In contrast to the collision-based phenomenological models for the daughter particle size distribution, this model predicts a symmetric distribution with the highest probability for equal-size particle breakage in accordance with the pure statistical model of Konno et al.183 (Figure 16). To summarize, purely statistical models for the daughter particle size distribution lack physical support. As for the breakage frequency models, the existing daughter particle size models based on eddy collision arguments rely on the assumption that turbulence consists of a collection of eddies that can be treated using relationships from the kinetic theory of gases. Martıńez-Bazań et al.,176 on the other hand, made a first

5142


Figure 16. Sketch of the breakage kernel of Konno et al.183 An equal daughter size distribution is predicted by this model. The kernel by Martıńez-Bazań176 is quite similar to this one, but the shape of the profile is slightly different close to the minimum and maximum values of d′3/d3.

approach toward developing an average microscopic kinematic/statistical model, where the physics were included through a daughter size distribution function expressed on the basis of kinematic arguments. Alternatively, one can develop kinematic/statistical models in which an effective daughter size distribution function represents an average of a large number of simulations of the physical breakage phenomena for a series of operating conditions. The deviating experimental observations reported in the literature indicate that further experimental investigations are needed to determine reliable correlations for which physical systems and under which flow conditions the equal and unequal breakage outcomes occur.

For most engineering calculations in the past, fluxes and integral results were considered sufficient interpreting growth, agglomeration, nucleation, and breakage data and for fitting the mechanistic kernels for solid particle suspensions. Therefore, only the first moments of the size density distribution function of the population (i.e., the mean, the standard deviation, etc.) are required. For the more complex fluid particle applications investigated in recent studies, local or semilocal size distributions are required, enabling better understanding of the physical phenomena involved and proper model validation. However, many of the present CFD codes still resort to integral formulations to minimize the computational time and memory requirements. It thus seems convenient to divide the numerical solution methods into two groups; methods solving the models for the moments of the size distribution and those solving the models for the size distribution itself. Methods for the Moments of the Distribution. Basically, the method of moments converts the set of PIDEs into a set of PDEs in which each PDE represents a given moment of the population size density function. Defining the jth moment of the population density function as

µj )

(111)

the PBE can be reduced to a set of PDEs by integrating the basic transport equation over the whole particle size distribution

∫0∞[∂t∂f + ∇r‚(R˙ f ) + ∇x‚(X˙ f ) - DC + BC - DB +

]

BB dj dd

for j ) 0, 1, ... (112)

or in terms of the moments

∂µj + ∂t

6. Numerical Discretization of the Population Balance Equation For the average microscopic formulations, the closures and the numerical discretizations are split, so that an optimal numerical solution method has to be found after the closure laws are derived. The numerical solution methods are, in general, problem-dependent and optimized for particular applications. In mathematical terms, the population balance equation (PBE) is classified as a nonlinear partial integrodifferential equation (PIDE). Because analytical solutions of this equation are not available for most cases of practical interest, several numerical solution methods have been proposed during the past two decades, as discussed by Williams and Loyalka128 and Ramkrishna.127 In general, the numerical solution of PIDEs consists of a three-step procedure. First, the basic set of PIDEs is discretized in the internal coordinates and expressed as a set of partial differential equations (PDEs). The PDEs are then discretized in time and in the physical space coordinates using standard PDE discretization techniques in a second step. Finally, the resulting set of algebraic equations is solved by use of a suitable solver.

∫0∞djf(r,d,t) dd

∫0∞[∇r‚(R˙ f)]dj dd + ∫0∞[∇x‚(X˙ f )]dj dd ) ∫0∞(DC + BC - DB + BB)dj dd

for j ) 0, 1, ... (113)

Assuming that the convective velocities are independent of the property coordinate (i.e., the particle diameter), one can integrate by parts to obtain

∂µj + djfR˙ |∞0 - ∇r‚R˙ ∂t

)

∫0∞jdj-1f dd + djX˙ f|∞0 ∞ ∇x‚X˙ ∫0 jdj-1 f dd

∫0∞(DC + BC - DB + BB)dj dd

(114)

for j ) 0, 1, ...

Because the f(r,d,t) function goes to zero at the maximum and minimum boundaries, the equation can be simplified and rewritten as

∂µj ∞ - j∇r‚(R˙ µj-1) - j∇x‚(X˙ µj-1) ) 0 (DC + BC - DB + ∂t for j ) 0, 1, ... (115) BB)dj dd

∫

with initial conditions µj(r,d,t ) 0) ) µˆ j(r,d) for j ) 0, 1, ...


The mathematical model characteristics of the source terms, i.e., DC, BC, DB, and BB, determine whether the system of equations is closed or not. To illustrate this problem, a pure breakage process is considered, i.e., DC ) BC ) 0, but DB * 0 and BB * 0. In addition, as an example only, a generalized formulation of the breakage death term (DB) is analyzed in the following discussion. Assuming that the breakage frequency, bB(d), can be represented by a polynomial N relation such as bB(d) ) ∑i)0 bBiPBi(d), where PBi(d) ) j d , integration of the source term gives

∫0

∞

N

bBi f(r,d,t)dj dd )

bB ∫0 dj f(r,d,t)dj dd ) ∑ i)0 ∞

i

N

bB µi+j ∑ i)0 i

(116)

This implies that higher-order moments are introduced, so the system of PDEs cannot be closed analytically. It is possible to show that similar effects will occur for the other source terms as well. This problem limits the application of the exact method of moments to the particular case where one has constant kernels only. In other cases, one has to introduce approximate closures in order to eliminate the higher-order moments, thereby ensuring that the transport equations for the moments of the PSD can be expressed in terms of the lower-order moments only (i.e., a modeling process very similar to turbulence modeling). For bubbly flows, most of the early papers adopted either a macroscopic population balance approach with an inherent discrete discretization scheme as described earlier or rather semiempirical transport equations for the contact area and/or the particle diameter. Actually, very few consistent source term closures exist for the microscopic population balance formulation. The existing models are usually solved using discrete semiintegral techniques, as will be outlined in the next subsection. Therefore, the approximate integral method is not widely used to solve the population balance model for bubbly flows, as the kernels for these systems are rather complex, so it is very difficult to eliminate the higherorder moments developing closures with sufficient accuracy. Hounslow et al.184 suggested that the difficulty incurred by the inclusion of complex closures for the source terms investigating nucleation, growth, and aggregation of particulate suspensions can be relaxed by discretizing the population balance equation. The moments of the size distribution were calculated from the discrete data resulting from the simulations. However, the direct solution of the equations containing physical source term parametrizations by use of conventional finite-difference techniques was complicated by the presence of the integral terms and resulted in computational loads far greater than desirable. Hounslow et al. thus placed substantial emphasis on correctly determining both the moments and the particle size distribution. The novel numerical technique developed, later referred to as the class method, reinforces the conservative properties of a few lower-order moments by transforming the equation through the introduction of a constant-volume correction factor that was found valid for several cases. Further details will be given later in discussing discrete semi-integral and differential

methods in general. The class method containing volume correction factors has not been widely used for bubbly flows mainly because of the complexity of the bubble breakage closures. When the population balance is written in terms of one internal coordinate (e.g., particle diameter or particle volume), the closure problem mentioned above for the moment equation can be successfully relaxed for solid particle systems by the use of a quadrature approximation. In the quadrature method of moments (QMOM) developed by McGraw185 for the description of sulfuric acid-water aerosol dynamics (growth), a certain type of quadrature function approximation is introduced to approximate the evolution of the integrals determining the moments. Marchisio et al.186,187 extended the QMOM for the application to aggregation-breakage processes. For the solution of crystallization and precipitation kernels, the size distribution function is expressed using an expansion in delta functions186,187 P

f(d,t) ≈

∑ ωRδ(d - dR)

(117)

R)0

Defining the ith moment of the population density function as

µi )

∫0∞di f(d,t) dd

(118)

and using a quadrature rule to approximate the integral yields

µi )

∫

∞ i d f(d,t) 0

P

dd ≈

ωR(t)diR ∑ R)1

(119)

Inserting the same quadrature approximation into the source term integrals provides an approximate numerical type of closure avoiding the closure problem of higher-order moments at the cost of model accuracy.185 This numerical approximation actually neglects the physical effects of the higher-order moments. No reports applying this procedure to bubbly flows have been found so far. Methods for the Distribution. This group of methods contains a large variety of solution strategies, but only a few popular techniques such as the finite-volume (FVM) methods and the direct quadrature method of moments (DQMOM) are discussed here. The main premise of these methods is that, in practice, one is not necessarily interested in the number density probability f1, but rather in the number density Ni (i.e., the number of bubbles of a particular size or size interval per unit volume). The methods within this category either divide the size domain into finite regions on which the number density is integrated to provide balances between the number of particles in each “bin” or solve the balance equations accurately for a limited number of discretization points on the size domain. A group of discrete techniques is based on the conventional finite-volume and finite-difference concepts and can be classified as finite-volume methods (FVMs). The wide use of the FVMs to solve PBEs is due to their simple construction and conservative characteristics, in a similar manner as the FVM is often preferred for formulating the CFD discretization algorithms. That is, the flux of a given property leaving through one face of

5144


the grid volume must equal the flux entering the neighboring ones. Adopting this method thus enforces the conservation of the conserved quantities. The growth terms are simply treated as flux processes, whereas the coalescence and breakage processes are more difficult to handle because the particles can disappear and appear in different regions of the domain. To fulfill the conservation properties considering these processes, it is necessary to employ some kind of ad hoc numerical “tricks” analogous to the linear source term manipulations in the CFD algorithms. The different tricks employed in the literature on PBEs have been characterized as the so-called discrete, class, and multigroup methods. In principle, the discrete and multigroup methods are very similar. A slight difference is observed in the way the number density is approximated. The discrete method approximates the functional values at discrete pivotal points within the size interval using delta functions, whereas the multigroup method divides the size spectrum into a number of continuous subintegrals or groups. In addition, the mean value theorem is sometimes applied in different ways. The multigroup method substitutes the mean value of the population density in each interval in each integrand by fi ) Ni/ (di+1 - di) and withdraws this factor from the integral, whereas the discrete method can be based on the multigroup approach or alternatively derived by estimating the kernels at the pivots and assigning the size interval definition to the remaining factors of the integrals. The two names used to refer to these methods originate from the parallel historical developments within two different fields of science. The discrete method was developed in chemical engineering applications investigating biopopulations, whereas the multigroup approach was developed in neutron and nuclear reactor physics. A well-established approach in nuclear engineering to solving the neutron density transport equations resorts to the multiple energy group equations using transfer cross sections. Similar ideas were later used to determine a detailed description of the dynamic interactions among groups of particles. Because some versions of the discrete or multigroup FVMs are adopted in nearly all multifluid CFD simulations, the main steps in the PDE discrete discretization procedures are outlined here.127 First, the dispersed phase is divided into size intervals according to some criterion; in this case, we have chosen the bubble diameter di. The interval [di, di+1] is denoted Ii, and the total number of particles in interval Ii is given by

Ni(t) ≡

∫dd

i+1

fi(d,t) dd

i

(120)

The population balance is integrated in particle size over the subinterval

∫dd

∂ ∂t )

i

)

∫d

i+1

∂f(d,t) dd + ∂t

∫dd i

i+1

(f(d,t)v) dd

i

(122)

di+1

[BB(d,t) - DB(d,t) + BC(d,t) - DC(d,t)] dd

i

∫dd

∫dd

[f(d,t)v] dd )

i+1

i

f(d,t) dd 〈v〉Ni ) Ni〈v〉Ni (123)

i+1

i

Inserting eqs 120 and 123 into eq 122 yields

∂Ni + ∇‚(Ni〈v〉Ni) ∂t )

∫dd

[BB(d,t) - DB(d,t) + BC(d,t) - DC(d,t)] dd (124)

i+1

i

The integrals (with respect to d′) in the source terms (see eqs 89-92) are expressed as the sums of integrals over subintervals. Replacing the integral in each source terms with a sum yields M

BB(d,t) )

∫d ∑ j)i

dj+1

v(d′) bB(d′)PB(d|d′) f1(d′,t) dd′ (125)

j

DB(d,t) ) bB(d) f1(d,t) BC(d,t) )

1i-1

∑∫ 2j)0 d

j

aC[(d3 - d′3)1/3,d′] f1[(d3 d′3)1/3,t] f1(d′,t) dd′ (127) M

DC(d,t) ) f1(d,t)

∫d ∑ j)0

dj+1 j

aC(d,d′) f1(d′,t) dd′

The class boundaries are independent of time and spatial position, so ∂/∂t and ∇ can be moved outside the integral

(128)

The above expressions contain the continuous bubble number probability density, f(d,t). The source terms must be expressed entirely in terms of the dependent variable Ni. This can be achieved by using the mean value theorem.127 Note, as mentioned before, at this point, the discrete method of Ramkrishna127 might deviate slightly from the multigroup method. The mean value theorem is used to cast the eqs 125129 entirely in terms of Ni and Nj. In each interval, the variables are replaced by a mean value. For example, using the discrete method, the mean value theorem can be used to express

∫dd

DC(d,t) dd )

i+1

∫dd i

M

i+1

f1(d,t)

∫d ∑ j)0

dj+1 j

aC(d,d′) f1(d′,t) dd′ dd

∇‚(f(d,t)vp) dd

(BB(d,t) - DB(d,t) + BC(d,t) - DC(d,t)) dd (121)

(126)

dj+1

i+1

di+1 i

∫dd

f(d,t) dd + ∇‚

The result of the integration is a balance equation for the number densities, Ni, in terms of the unknown number density probability, f1(d,t), and is, hence, still unresolvable. A number-average particle velocity can be introduced to simplify the expression for the convective terms

i

∫dd

∫d

i+1

i

∫d

≈ aC(xi,xj)

di+1 i

M

f1(d,t) dd

∫d ∑ j)0

dj+1 j

f1(d′,t) dd′ ) M

Ni

aC(xi,xj)Nj ∑ j)0

(129)

where xi and xj are pivot points in Ii and Ij, respectively.


The pivot concentrates the particles in the interval at a single representative point. Thus, one can write the number probability density f1(d,t) as being given by

[ ]

M

f1(d,t) )

∑ i)0

1

Niδ(d - xi)

m3 (m)

∫d

M

BB(d,t) dd )

i

∫d

Njv(xj) bB(xj)∫d ∑ j)i

i

PB(d|xj) dd

∫dd

BC(d,t) dd )

i+1

i

DB(d,t) dd ) bB(xi)Ni 1i-1

∑Nj ∑ 2j)0

NkaC(xk,xj)

∫d

(133)

i

NjaC(xi,xj) ∑ j)0

(134)

+

Njv(xj) bB(xj)∫d ∑ j)i

di+1 i

PB(d|xj) dd - bB(xi)Ni

i-1

∑Nj ∑

2j)0

M

+

NkaC(xk,xj) - Ni

NjaC(xi,xj) ∑ j)0 i ) 0, 1, ..., M (135)

Equation 135 is not necessarily conservative because of the finite (i.e., in practice, rather coarse) size grid resolution, and some sort of numerical trick must be used to enforce the conservative properties. It is mainly at this point in the formulation of the numerical algorithm that the class method of Hounslow et al.,184 the discrete method of Ramkrishna,127 and the multigroup approach used by Carrica et al.,13 among others, differ to some extent, as discussed earlier. The problem in question is related to the birth terms only. The formation of a bubble of size d′ in size range (xi, xi+1) due to breakage or coalescence is represented by assigning fractions γ1(d,xi) and γ2(d,xi+1) to bubble population at pivots xi and xi+1, respectively. This is necessary because not all coalescence and breakage result in a bubble that has a legitimate size. For the nonvalid daughter particles, one has to distribute the new bubble in fractions γ1(d,xi) and γ2(d,xi+1) of the two neighboring size pivots. To determine these two fractions, one needs two balance equations. The first balance equation relates to the total volume or mass of the bubbles, to ensure mass conservation. The second balance equation is the number balance for the bubbles involved in the breakage and coalescence processes

γ1(d,xi) Vb,i(d) + γ2(d,xi+1) Vb,i+1(d) ) Vb(d) (136) γ1(d,xi) + γ2(d,xi+1) ) 1

xi+1 i

+

1i-1

∑Nj ∑

+

PB(d|xj) dd - bB(xi)Ni

Nk γ1(d,xi)a(xk,xj)

(xj+xk)∈Ii

1i-1

M

∑Nj ∑ 2j)0

Nkγ2(d,xi)a(xk,xj) - Ni

(137)

NjaC(xi,xj) ∑ j)0

(139)

The breakage function, PB(d|d′) dd, is the number of particles formed between size d′ and d′ + dd′ divided by the total number of size d particles broken. At the pivots, the corresponding breakage function is defined as PB(xi|xk)

PB(xi|xk) )

M

(xj+xk)∈Ii

γ2(d,xi) Njv(xj) bB(xj)∫x ∑ j)i

(xj+xk)∈Ii-1

dNi + ∇‚(Ni〈v〉Ni) dt M

xi i

M

DC(d,t) dd ) Ni

γ1(d,xi) Njv(xj) bB(xj)∫x -1 PB(d|xj) dd ∑ j)i

2j)0

and inserting into eq 124 gives

1

)

(xj+xk)∈Ii

di+1

)

dNi + ∇‚(Ni〈v〉Ni) dt

(131) (132)

(138)

The fractions γ1(d,xi) and γ2(d,xi+1) need to be embedded in the source terms as part of the discretization procedure. Thus, the source terms are modified to

M

di+1

di+1 i

γ1(d,xi) ab,i + γ2(d,xi+1) ab,i+1 ) ab(d)

(130)

The RHS of eq 131 is zero for d * xi and equals Ni when d ) xi. Applying the mean value theorem to the source terms yields di+1

Alternatively, the second balance equation is sometimes substituted by a bubble surface balance as one wants to ensure good estimates of the contact area

∫xx

i

PB(d|xk) dd +

i-1

∫xx i

i+1

PB(d|xk) dd (140)

In this method, the smallest bubbles do not break up, and the largest bubbles are not involved in the coalescence process. To cover a broad range in bubble volume (diameter), a large number of classes is needed, making this algorithm rather time-consuming. However, the method provides information on the bubble size distribution that is needed in the multifluid framework and can also be used for proper validation of the source term closures. The simpler method of moments containing transport equations for only a few moments such as the mean bubble size, the variance, etc., cannot be used in a multifluid framework because of the lack of any bubble size resolution and can thus be validated only in an average sense even if experimental data on the size distribution are provided. For very simple distribution functions, it might be possible to reconstruct the physical size distributions with sufficient accuracy utilizing the information provided by a few moments only (i.e., as provided by the moment models). However, as the existing integral models provide two to five moments only, it is believed that the complex bubble size distributions cannot be reproduced with sufficient accuracy. It is then an open question as to whether the computational costs of solving transport equations for about 10-15 moments, ensuring sufficient accuracy, are less than the corresponding costs of solving for the whole size distributions using spectral methods. Further work in our group continues to elucidate this issue. Another method representing an extension of the QMOM method has obtained increasing attention for particulate systems during the past several years. According to Fan et al.,101 one of the main limitations

5146


of the QMOM is that the solid phase is represented through the moments of the distribution, so the phaseaverage velocity of the different solid phases must be used to solve the transport equations for the moments. Thus, to use this method in the context of multiphase flows, it is necessary to extend QMOM to handle cases where each particle size is convected by its own velocity. To address these issues, a direct quadrature method of moments (DQMOM) has been formulated by Fan et al.101 DQMOM is based on the direct solution of the transport equations for weights and abscissas of the quadrature approximation (i.e., quadrature approximations for the moments). Moreover, each node of the quadrature approximation is treated as a distinct solid phase. DQMOM is thus used as a module in multifluid models describing polydisperse solid phases undergoing segregation, growth, aggregation, and breakage processes in the context of CFD simulations. However, in the literature, there seems to be some inconsistency regarding whether the QMOM is really well posed or not.188 To date, no reports have been published evaluating the behavior of this procedure for bubbly flows. 7. Concluding Remarks The past two decades have seen substantial developments in both experimental description and modeling of bubbly flows. However, the physical understanding of the local spatial and temporal scales is still very limited, and thus, the modeling tools are restricted to very low gas fractions. Average multifluid models determining pressure and phase velocities combined with an integrated population balance module predicting the phase and bubble size distributions have been found to represent a tradeoff between accuracy and computational efforts for practical applications. During the past decade, it has become clear that, to resolve the physics of bubble column flows, dynamic and three-dimensional fluid dynamic models are required. Considering the fluid dynamic part of the model, the general picture from the literature is that time-averaged velocity fields are reasonably well-predicted with models of this type. The prediction of phase distribution phenomena, on the other hand, is still a problem, particularly at high gas flow rates. The limiting steps in the model development are the formulation of closure relations or closure laws determining turbulence effects, interfacial transfer fluxes, and bubble coalescence and breakage processes. Improved formulations of the boundary conditions are also needed for flows with higher gas volume fractions, especially at the inlet because of the unresolved interaction between the fluid dynamics and the bubble size and phase distributions and at the outlet because of the highly recirculating flow phenomena observed at the top of the column. In a similar manner, the wall boundary needs further development because of the complex interaction between multiphase fluid dynamics, colloid chemistry, bubble coalescence and breakage, and the turbulence boundary layer structure. Finally, the tremendous computational requirements in terms of both CPU time and memory have to be reduced by developing more efficient parallel algorithms. The numerical accuracy and stability of these schemes are also of critical importance, as the numerical solution algorithm deter-

mines an integrated part of the bubble column model development. Considering the population balance modeling of bubbly flows, the general picture from the literature is that these models are at a very early stage of development and that the physicochemical fluid dynamics involved are not yet sufficiently understood. For bubbly flows, there is even some confusion regarding the formulation of the governing population balance equation. Most of the reports on bubbly flows resort to a rather empirical macroscopic formulation, whereas only a few papers adopt microscopic and more generalized concepts. Most of the existing coalescence and breakage closures are thus preferably formulated directly within the macroscopic framework. The breakage closures are often based on an incomplete discrete (eddy) interpretation of the turbulence structure. These formulations are very difficult or impossible to validate directly on the bubblescale level and are considered semiempirical in nature. In addition, the commonly used exponential relations for both the coalescence and breakage probabilities have not been validated in a sufficient manner for bubble column flows. The coalescence closures are also limited because of the lack of understanding in formulating accurate coalescence criteria. Proper validation of the existing breakage criteria for bubble column flows is also recommended. Both the turbulent breakage and coalescence closures are very sensitive to the turbulence energy dissipation rate, which is difficult to estimate with sufficient accuracy using the standard turbulence models. Further, by the use of number probability density functions within the population balance framework, an unfortunate inconsistency occurs regarding the interpretation of the velocity variable involved in the integrated multifluid/population balance model. In the multifluid model formulation, the velocity variable is usually a mass-average quantity, whereas the corresponding variable within the population balance is a number-density-average variable. It would probably be advantageous to formulate a population balance for the particle mass rather than the number density to ensure better consistency between the two model parts having somewhat different origins. However, it is expected that the next generation of population balances will be formulated adopting the more fundamental statistical concepts on the microscopic level and that the coalescence and breakage closures will be reformulated in a consistent manner. The formulation of optimal numerical solution methods solving the integrated multifluid/population balance model is not a trivial task. In most cases, the population balance part of the model is solved by adopting the simplest methods originally developed for solid particle analysis. Because the mathematical properties of the physical processes and closures involved can deviate considerably, the choice of solution strategy needs further consideration. In parallel to the modeling work, intensive efforts have to be focused on well-planned experimental programs to enable better understanding of the phenomena to be described by the microscopic modeling process and for model validation. In this respect, extended integrations and knowledge transfer between the modeling and experimental groups will be beneficial.


Acknowledgment This work received support from the Research Council of Norway (Program of Supercomputing) through a grant of computer-time. The Ph.D. fellowships (H.L. and C.A.D.) financed by the Research Council of Norway through a Strategic University Program (CARPET) are gratefully appreciated. Nomenclature Latin Letters aC(x,r;x′,r′,Y,t) ) coalescence frequency (general), i.e., fraction of particles of sizes d and d′ coalescing per unit time (l/s) aC(d,d′,Y) ) average coalescence frequency (m3/s) A ) parameter in the drag coefficient relation bB(d) ) breakage frequency (l/s) BB(x,r,Y,t), BC(x,r,Y,t) ) birth rates due to breakage and coalescence (general), respectively {1/[s m3 (x)]} BB(d), BC(d) ) average birth rates due to breakage and coalescence, respectively, with d as the internal coordinate {1/[s m3 (m)]} c ) concentration of surfactant species (kmol/m3) CD ) drag coefficient Cf ) wall friction coefficient Ck ) Kolmogorov energy spectrum parameter ()1.5) CL ) lift coefficient CV ) added-mass force coefficient Cwb ) wall boundary coefficient Cw1 ) wall lift coefficient Cw2 ) wall lift coefficient d ) bubble diameter (m) dc ) critical diameter (m) di ) diameter of particle in interval i (m) d(i) ) diameter of particle in class i (m) d′ ) diameter of daughter bubble (m) d′′ ) diameter of the smallest daughter bubble (m) d ˜ ) diameter of complementary daughter particle, d ˜ ) (d3 - d′3)1/3 (m) dmax, dmin ) boundary diameters of class to be split (m) ds ) Sauter mean diameter (m) dVx ) infinitesimal volume in property space dVr ) infinitesimal volume in physical space (m3) DB(x,r,Y,t), DC(x,r,Y,t) ) death rates due to breakage and coalescence (general), respectively {1/[s m3 (x)]} DB(d), DC(d) ) average death rates due to breakup and coalescence, respectively, with d as the internal coordinate {1/[s m3 (m)]} ej(di,λ) ) mean kinetic energy of an eddy (J) es(di,dj) ) increase in surface energy due to binary breakage (J) E ) constant of integration Eeddies(λ), Espectra(λ) ) kinetic energy contained within eddies of size between λ and λ + dλ {J/[m3 (m)]} E(k) ) kinetic energy contained within eddies of wavenumbers between k and k + dk [m2 (m)/s2] Es ) minimum energy needed to deform a particle (J) E0 ) Eo¨tvo¨s number f(d,t) ) particle number density probability {1/[m3 (m)]} f(x,r,t) ) particle number density (general) {1/[m3 (m)]} fvij ) breakage volume fraction ()dj3/di3) fλ ) number of eddies of size between λ and λ + dλ {1[m3 s (m)] h ) film thickness (m) hC ) effective swept volume rate (m3/s) hf ) final film thickness (m) h0 ) initial film thickness (m) k ) turbulent kinetic energy (m2/s2), or wavenumber (1/ m) ni ) number density in the interval i (1/m3)

N ) total number of particles per unit volume of physical space (1/m3) Ni ) number density of particles in size interval i (1/m3) p ) pressure (N/m2) pB(di:dj) ) breakage probability (or breakage efficiency) pB(di:dj,λ) ) breakage probability (or breakage efficiency), dependent on eddy size pC ) coalescence probability (or coalescence efficiency) PB(d|d′,Y) ) average probability density function, or daughter size distribution function [1/(m)] PB(x,r|x′,r′,Y,t) ) probability density function (general), or daughter size distribution function (general) {1/[m3 (X)]} rb ) bubble radius (m) R ) molar gas constant [J/(mol K)] Rd ) radius of liquid disk between coalescing bubbles (m) ReB ) particle Reynolds number Sc ) Schmidt number t ) time coordinate (s) tb ) particle breakage time (s) T ) temperature (K) vj λ ) mean eddy velocity expressed by the Maxwell distribution function (m/s) v2(λ) ) second-order structure function (m2/s2) V ) volume, physical or abstract Vi ) volume of a bubble in class i We ) Weber number xi, xj ) pivotal points in Ii and Ij y ) distance from the wall (m) y+ ) inner boundary layer length scale (m) z ) axial coordinate (m) Vectors Fi,k ) interfacial force (i ) D, V, L, TD, W) [kg/(m2 s2)] g ) acceleration of gravity, or gravity force per unit mass (m/s2) Mk,l ) interfacial momentum transfer to phase k from phase l [kg/(m2 s2)] nw ) unit vector normal to the wall n ) outward-directed normal unit vector r ) position vector in physical space (m) R4 ) external coordinates (m) v ) velocity (m/s) vd ) velocity of dispersed phase (m/s) vk ) velocity of phase k (m/s) vp ) velocity of particles (m/s) vrel ) relative velocity between the phases (m/s) x ) (internal) property vector X4 ) internal coordinates (m) Y ) vector representing the continuous phase variables Greek letters Rk ) volume fraction of phase k β ) empirical model parameter β˜ ) empirical model parameter in turbulence theory γ ) bubble fraction Γ ) incomplete gamma function Γk,l ) net mass-transfer flux to phase k from all other l phases δ ) correction of redundancy factor, or delta function ) turbulent energy dissipation rate λ ) eddy size, or turbulent length scale (m) Λ ) bubble diameter ratio µj ) jth moment of the population density function µk ) viscosity of phase k [kg/(m s)] ν(d,t) ) average number of particles formed from the breakup ξ ) eddy/bubble size ratio ()λ/di) ξij ) bubble/bubble size ratio () dj/di) Fk ) density of phase k (kg/m3) σA ) effective collision cross-sectional area (m2)

5148


σI ) surface tension (N/m) σk, σ ) Prandtl numbers σ j k ) viscous stress tensor (Pa) σs ) surface restoring pressure (Pa) σt ) average deformation stress (Pa) φ ) net source term χ ) energy ratio [es(di,dj)/ej(di,λ)] χc ) critical energy ratio ψ ) extensive property ωB(di:dj) ) eddy-bubble collision density [1/(s m3)] ωIB(di:λ) ) eddy-bubble collision probability density {1/[s m3 (m)]} ωC ) bubble-bubble collision density {1/[m3 s (m) (m)] ωBC ) buoyancy collision density [1/(s m3)] 3 ωLS C ) collision density due to laminar shear [1/(s m )] ωf1 ) single-particle collision density [1/s (m)] ωTC ) collision density resulting from turbulent motion (s m3) ΩB ) breakage rate [1/(s m3)] ΩC ) coalescence rate [1/(s m3)]

Literature Cited (1) Deckwer, W.-D.; Schumpe, A. Bubble columnsThe state of the art and current trends. Int. Chem. Eng. 1987, 27, 405. (2) Deckwer, W.-D. Bubble Column Reactors; John Wiley & Sons: Chichester, U.K., 1992. (3) Deckwer, W.-D.; Schumpe, A. Improved Tools for Bubble Column Reactor Design and Scale-up. Chem. Eng. Sci. 1993, 48, 889. (4) Dudukovic, M. P. Trends in Catalytic Reaction Engineering. Catal. Today 1999, 48, 5. (5) Dudukovic, M. P.; Larachi, F.; Mills, P. L. Multiphase reactorssRevisited. Chem. Eng. Sci. 1999, 54, 1975. (6) Dudukovic, M. P.; Larachi, F.; Mills, P. L. Multiphase Catalytic Reactors: A Perspective on Current Knowledge and Future Trends. Catal. Rev. 2002, 44, 123. (7) Shah, Y. T.; Sharma, M. M. Gas-Liquid-Solid Reactors. In Chemical Reaction and Reaction Engineering; A Series of Reference Books and Textbooks; Carberry, J. J., Varma, A., Eds.; Marcel Dekker: New York, 1987; p 667. (8) Jakobsen, H. A.; Sannæs, B. H.; Grevskott, S.; Svendsen, H. F. Modeling of Vertical Bubble-Driven Flows. Ind. Eng. Chem. Res. 1997, 36, 4052. (9) Joshi, J. B. Computational flow modeling and design of bubble column reactors. Chem. Eng. Sci. 2001, 56, 5893. (10) Rafique, M.; Chen, P.; Dudukovic, M. P. Computational modeling of gas-liquid flow in bubble columns. Rev. Chem. Eng. 2004, 20, 225. (11) Jakobsen, H. A. Phase distribution phenomena in twophase bubble column reactors. Chem. Eng. Sci. 2001, 56, 1049. (12) van den Akker, H. E. A. Computational Fluid Dynamics: More Than a Promise to Chemical Reaction Engineering, Paper Presented at the CHISA 2000 Conference, Praha, Czech Republic, Aug 27-31, 2000. (13) Carrica, P. M.; Drew, D.; Bonetto, F.; Lahey, R. T., Jr. A Polydisperse Model for Bubbly Two-Phase Flow Around a Surface Ship. Int. J. Multiphase Flow 1999, 25, 257. (14) Oey, R. S.; Mudde, R. F.; van den Akker, H. E. A. Sensitivity Study on Interfacial Closure Laws in Two-Fluid Bubbly Flow Simulations. AIChE 2003, 49, 1621. (15) Sokolichin, A.; Eigenberger, G.; Lapin, A. Simulation of Buoyancy Driven Bubbly Flow: Established Simplifications and Open Questions. AIChE J. 2004, 50, 24. (16) Behzadi, A.; Issa, R. I.; Rusche, H. Modelling of dispersed bubble and droplet flow at high phase fractions. Chem. Eng. Sci. 2004, 59, 759. (17) Svendsen, H. F.; Jakobsen, H. A.; Torvik, R. Local Flow Structures in Internal Loop and Bubble Column Reactors. Chem. Eng. Sci. 1992, 47, 3297. (18) Sanyal, J.; Vasquez, S.; Roy, S.; Dudukovic, M. P. Numerical simulation of gas-liquid dynamics in cylindrical bubble column reactors. Chem. Eng. Sci. 1999, 54, 5071. (19) Krishna, R.; van Baten, J. M. Rise Characteristics of Gas Bubbles in a 2D Rectangular Column: VOF Simulations vs Experiments. Int. Commun. Heat Mass Transfer 1999, 26, 965.

(20) Jakobsen, H. A. On the Modeling and Simulation of Bubble Column Reactors Using a Two-Fluid Model. Dr. Ing. Thesis, Norwegian Institute of Technology, Trondheim, Norway, 1993. (21) Tzeng, J.-W.; Chen, R. C.; Fan, L.-S. Visualization of Flow Characteristics in a 2D Bubble Column and Three-Phase Fluidized Bed. AIChE J. 1993, 39, 733. (22) Lin, T.-J.; Reese, J.; Hong, T.; Fan, L.-S. Quantitative Analysis and Computation of Two-Dimensional Bubble Columns. AIChE J. 1996, 42, 301. (23) Chen, R. C.; Reese, J.; Fan, L.-S. Flow Structure in ThreeDimensional Bubble Column and Three-Phase Fluidized Bed. AIChE J. 1994, 40, 1093. (24) Sokolichin, A.; Eigenberger, G. Gas-liquid flow in bubble columns and loop reactors: Part I. Detailed modeling and numerical simulation. Chem. Eng. Sci. 1994, 49, 5735. (25) Sokolichin, A.; Eigenberger, G.; Lapin, A.; Lu¨bbert, A. Dynamical numerical simulation of gas-liquid two-phase flows. Euler/Euler versus Euler/Lagrange. Chem. Eng. Sci. 1997, 52, 611. (26) Sokolichin, A.; Eigenberger, G. Applicability of the k- turbulence model to the dynamic simulation of bubble columns: Part I. Detailed numerical simulations. Chem. Eng. Sci. 1999, 54, 2273. (27) Becker, S.; Sokolichin, A.; Eigenberger, G. Gas-Liquid Flow in Bubble Columns and Loop Reactors: Part II. Comparison of Detailed Experiments and Flow Simulations. Chem. Eng. Sci. 1994, 49, 5747. (28) Becker, S.; De Bie, H.; Sweeney, J. Dynamics flow behaviour in bubble columns. Chem. Eng. Sci. 1999, 54, 4929. (29) Krishna, R.; van Baten, J. M. Scaling up Bubble Column Reactors with Highly Viscous Liquid Phase. Chem. Eng. Technol. 2002, 25, 1015. (30) Pfleger, D.; Becker, S. Modeling and simulation of the dynamic flow behaviour in a bubble column. Chem. Eng. Sci. 2001, 56, 1737. (31) Olmos, E.; Gentric, C.; Vial, Ch.; Wild, G.; Midoux, N. Numerical simulation of multiphase flow in bubble column reactors. Influence of bubble coalescence and breakup. Chem. Eng. Sci. 2001, 56, 6359. (32) Lehr, F.; Millies, M.; Mewes, D. Bubble-Size Distributions and Flow Fields in Bubble Columns. AIChE J. 2002, 48, 2426. (33) Olmos, E.; Gentric, C.; Midoux, N. Numerical description of flow regime transitions in bubble column reactors by multiple gas-phase model. Chem. Eng. Sci. 2003, 58, 2113. (34) Bertola, F.; Vanni, M.; Baldi, G. Application of Computational Fluid Dynamics to Multiphase Flow in Bubble Columns. Int. J. Chem. React. Eng. 2003, 1, 1. (35) Krishna, R.; van Baten, J. M. Scaling up Bubble column reactors with the aid of CFD. Trans. Inst. Chem. Eng. 2001, 79, 283. (36) Lehr, F.; Mewes, D. A transport equation for the interfacial area density applied to bubble columns. Chem. Eng. Sci. 2001, 56, 1159. (37) Mudde, R. F.; Lee, D. J.; Reese, J.; Fan, L.-S. Role of Coherent Structures on Reynolds Stresses in a 2-D Bubble Column. AIChE J. 1997, 43, 913. (38) Mudde, R. F.; Simonin, O. Two- and three-dimensional simulations of bubble plume using a two-fluid model. Chem. Eng. Sci. 1999, 54, 5061-5069. (39) Pfleger, D.; Gomes, S.; Gilbert, N.; Wagner, H.-G. Hydrodynamic simulations of laboratory scale bubble columns fundamental studies of the Eulerian-Eulerian modelling approach. Chem. Eng. Sci. 1999, 54, 5091. (40) Deen, N. G.; Solberg, T.; Hjertager, B. H. Large Eddy Simulation of the Gas-Liquid Flow in a Square Cross-Sectioned Bubble Column. Chem. Eng. Sci. 2001, 56, 6341. (41) Buwa, V. V.; Ranade, V. V. Dynamics of gas-liquid flow in a rectangular bubble column: Experimental and single/multigroup CFD simulations. Chem. Eng. Sci 2002, 57, 4715. (42) Buwa, V. V.; Ranade, V. V. Mixing in Bubble Column Reactors: Role of Unsteady Flow Structures. Can. J. Chem. Eng. 2003, 81, 402. (43) Bove, S.; Solberg, T.; Hjertager, B. H., Numerical Aspects of Bubble Column Simulations. Int. J. Chem. React. Eng. 2004, 2. (44) Lopez de Bertodano, M. Turbulent Bubbly Two-Phase Flow in a Triangular Duct. Ph.D. Thesis, Rensselaer Polytechnic Insitute, Troy, New York, 1992. (45) Anglart, H.; Andersson, S.; Podowski, M. Z.; Kurul, N. An Analysis of Multidimensional Void Distribution in Two-Phase Flows. In Proceedings of the Sixth International Topic Meeting on

Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5149 Nuclear Reactor Thermal HydraulicssNURETH 6; Courtaud, M., Delhaye, J. M., Eds.; 1993; Vol. 1, pp 139-153. (46) Antal, S. P.; Lahey, R. T., Jr.; Al-Dahhan, M. H. Simulating Churn-Turbulent Flows in a Bubble Column using a Three Field, Two-Fluid Model. Paper Presented at the 5th International Conference on Multiphase Flow, ICMF′04, Yokohama, Japan, May 30-Jun 4, 2004; Paper No. 182. (47) Dhotre, M. T.; Joshi, J. B. Two-Dimensional CFD Model for the Prediction of Flow Pattern, Pressure Drop and Heat Transfer Coefficient in Bubble Column Reactors. Trans. Inst. Chem.Eng., Chem. Eng. Res. Des. 2004, 82, 689. (48) Tchen, C. M. Mean value and correlation problems connected with the motion of small particles suspended in a turbulent fluid. Ph.D. Thesis, TU Delft, Delft, The Netherlands, 1947. (49) Simonin, O. Eulerian formulation for particle dispersion in turbulent two-phase flows. In Fifth Workshop on Two-Phase Flow Predictions; Sommerfeld, M., Wennerberg, P., Eds.; FRG: Erlangen, Germany, 1990; pp 156-166. (50) Simonin, O.; Viollet, P. L. Prediction of an oxygen droplet pulverization in a compressible subsonic coflowing hydrogen flow. Numer. Methods Multiphase Flows 1990, 91, 73. (51) Viollet, P. L.; Simonin, O. Modelling dispersed two-phase flows: Closure, validation and software development. Appl. Mech. Rev. 1994, 47, S80. (52) Tomiyama, A.; Matsouka, T.; Fukuda, T.; Sakaguchi, T. A Simple Numerical Method for Solving an Incompressible TwoFluid Model in a General Curvelinear Coordinate System. Adv. Multiphase Flow 1995, Elsevier: 241. (53) Tomiyama, A. Struggle with computational bubble dynamics. Presented at the Third International Conference on Multiphase Flow, ICMF′98; Lyon, France, Jun 8-12, 1998. (54) Tomiyama, A.; Shimada, N. A Numerical Method for Bubbly Flow Simulation Based on a Multi-Fluid Model. J. Pressure Vessel Technol. 2001, 23, 510. (55) Sato, Y.; Sekoguchi, K. Liquid Velocity Distribution in TwoPhase Bubbly Flow. Int. J. Multiphase Flow 1975, 2, 79. (56) Sato, Y.; Sadotomi, M.; Sekoguchi, K. Momentum and Heat Transfer in Two-Phase Bubble Flow. Int. J. Multiphase Flow 1981, 7, 167. (57) Jakobsen, H. A. Numerical Convection Algorithms and Their Role in Eulerian CFD Reactor Simulations. Int. J. Chem. React. Eng. 2003, 1, 1. (58) Bech, K. Modeling of bubble-driven flow applied to electrolysis cells. Paper Presented at the 5th International Conference on Multiphase Flow, ICMF′04, Yokohama, Japan, May 30-Jun 4, 2004; Paper 404. (59) Cartland Glover, G. M.; Generalis, S. C. The modeling of buoyancy driven flow in bubble columns. Chem. Eng. Process. 2004, 43, 101. (60) van den Akker, H. E. A. Coherent structures in multiphase flows. Powder Technol. 1998, 100, 123. (61) Harteveld, W. K.; Mudde, R. F.; van den Akker, H. E. A. Dynamics of a Bubble Column: Influence of Gas Distribution of Coherent Structures. Can. J. Chem. Eng. 2003, 81, 389. (62) Politano, M. S.; Carrica, P. M.; Larreteguy, A. E. Gas/ Liquid Polydisperse Two-Phase Flow Modeling for Bubble Column Reactors. Paper Presented at the 3rd Mercosur Congress on Process Systems Engineering 1st Mercosur Congress on Chemical Engineering, Sante Fe, Argentina, Sep 16-20, 2001. (63) Politano, M. S.; Carrica, P. M.; Converti, J. A model for turbulent polydisperse two-phase flow in vertical channels. Int. J. Multiphase Flow 2003, 29, 1153. (64) Lo, S. Application of population balance to CFD modelling of bubbly flow via the MUSIG model. AEA Technol. 1996, AEAT1096. (65) Lo, S. Some recent developments and applications of CFD to multiphase flows in stirred reactors. Presented at the AMIFESF Workshop: Computing Methods for Two-Phase Flow, Aussois, France, Jan 12-14, 2000. (66) Bertola, F.; Grundseth, J.; Hagesaether, L.; Dorao, C.; Luo, H.; Hjarbo, K. W.; Svendsen, H. F.; Vanni, M.; Baldi, G.; Jakobsen, H. A. Numerical Analysis and Experimental Validation of Bubble Size Distribution in Two-Phase Bubble Column Reactors. Multiphase Sci. Technol. 2005, 17 (1-4), 1-23. (67) Kocamustafaogullari, G.; Ishii, M. Foundation of the interfacial area transport equation and its closure relations. Int. J. Heat Mass Transfer 1995, 38, 481.

(68) Millies, M.; Mewes, D. Interfacial area density in bubbly flow. Chem. Eng. Sci. 1999, 38, 307. (69) White, A. J.; Hounslow, M. J. Modelling droplet size distributions in polydispersed wet-stream flows. Int. J. Heat Mass Transfer 2000, 43, 1873. (70) Luo, H. Coalescence, break-up and liquid circulation in bubble column reactors. Dr. Ing. Thesis, Department of Chemical Engineering, The Norwegian Institute of Technology, Trondheim, Norway, 1993. (71) Luo, H.; Svendsen, H. F. Theoretical Model for Drop and Bubble Breakup in Turbulent Dispersions. AIChE J. 1996, 42, 1225. (72) Hagesaether, L.; Jakobsen, H. A.; Svendsen, H. F. Theoretical Analysis of Fluid Particle Collisions in Turbulent Flow. Chem. Eng. Sci. 1999, 54, 4749. (73) Hagesaether, L.; Jakobsen, H. A.; Hjarbo, K.; Svendsen, H. F. A Coalescence and Breakup Module for Implementation in CFD-codes. Comput.-Aided Chem. Eng. 2000, 8, 367. (74) Hagesaether, L.; Jakobsen, H. A.; Svendsen, H. F. A Model for Turbulent Binary Breakup of Dispersed Fluid Particles. Chem. Eng. Sci. 2002, 57, 3251. (75) Hagesaether, L.; Jakobsen, H. A.; Svendsen, H. F. Modeling of the Dispersed-Phase Size Distribution in Bubble Columns. Ind. Eng. Chem. Res. 2002, 41, 2560. (76) Jakobsen, H. A.; Hjarbo, K. W.; Svendsen, H. F. Bubble size distributions, bubble velocities and local gas fractions from conductivity measurements; Report STF21 F93103; SINTEF Applied Chemistry: Trondheim, Norway, 1994. (77) Prince, M. J.; Blanch, H. W. Bubble Coalescence and Break-Up in Air-Sparged Bubble Columns. AIChE J. 1990, 36, 1485. (78) Ranade, V. V. Computational Flow Modeling for Chemical Reactor Engineering; Academic Press: San Diego, 2002. (79) Ramakrishnan, S.; Kumar, R.; Kuloor, N. R. Studies in bubble formation-I: Bubble formation under constant flow conditions. Chem. Eng. Sci. 1969, 24, 731. (80) Geary, N. W.; Rice, R. G. Bubble Size Prediction for Rigid and Flexible Spargers. AIChE J. 1991, 37, 161. (81) Gresho, P. M.; Sani, R. L. On pressure boundary conditions for the incompressible Navier-Stokes equations. Int. J. Numer. Methods Fluids 1987, 7, 1111. (82) Padial, N. T.; VanderHeyden, W. B.; Rauenzahn, R. M.; Yarbro, S. L. Three-dimensional simulation of a three-phase drafttube bubble column. Chem. Eng. Sci. 2000, 55, 3261. (83) Troshko, A. A.; Hassan, Y. A. Law of the wall for two-phase turbulent boundary layers. Int. J. Heat Mass Transfer 2001, 44, 871. (84) Troshko, A. A.; Hassan, Y. A. A two-equation turbulence model of turbulent bubbly flows. Int. J. Multiphase Flow 2001, 27, 1965. (85) Oliveira, P. J.; Issa, R. I. On the numerical treatment of interfaphase forces in two-phase flow. In Numerical Methods in Multiphase Flows; Crowe, C. T., Ed.; ASME Press: New York, 1994; Vol. 185, p 131. (86) Spalding, D. B. The calculation of free-convection phenomena in gas-liquid mixtures. In Turbulent Buoyant Convection; Hemisphere: Washington, DC, 1977; p 569. (87) Spalding, D. B. Numerical computation of multiphase fluid flow and heat transfer. In Recent Advances in Numerical Methods in Fluids; Taylor, C., Morgan, K., Eds.; Pineridge Press: Swansea, U.K., 1980; p 139. (88) Kelly, J. M.; Stewart, C. W.; and Cuta, J. M. VIPRE-02s A two-fluid thermal-hydraulics code for reactor core and vessel analysis: Mathematical modeling and solution methods. Nucear Technol. 1992, 100, 246. (89) Kunz, R. F.; Siebert, B. W.; Cope, W. K.; Foster, N. F.; Antal, S. P.; Ettorre, S. M. A coupled phasic exchange algorithm for three-dimensional multi-field analysis of heated flows with mass transfer. Comput. Fluids 1998, 27, 741. (90) Lathouwers, D.; van den Akker, H. E. A. A numerical method for the solution of two-fluid model equations. In Numer. Methods Multiphase Flows 1996, 236, 121. (91) Lathouwers, D. Modeling and simulation of turbulent bubbly flow. Ph.D. Thesis, The Technical University of Delft, Delft, The Netherlands, 1999. (92) Antal, S. P.; Ettorre, S. M.; Kunz, R. F.; Podowski, M. Z. Development of a Next Generation Computer Code for the Prediction of Multicomponent Multiphase Flows. Presented at the International Meeting on Trends in Numerical and Physical

5150


Modeling for Industrial Multiphase Flow 2000, Cargese, France, Sep 27-29, 2000. (93) Kunz, R. F.; Yu, W.-S.; Antal, S. P.; Ettorre, S. M. An Unstructured Two-Fluid Method Based on the Coupled Phasic Exchange Algorithm; Report AIAA-2001-2672; American Institute of Aeronautics and Astronautics: Herndon, VA, 2001. (94) Vasquez, S. A.; Ivanov, V. A. A phase coupled method for solving multiphase problems on unstructured meshes. In Proceedings of ASME FEDSM′00; ASME Press: New York, 2000; p 1. (95) van Santen, H.; Lathouwers, D.; Kleijn, C. R.; van den Akker, H. E. A. Influence of segregation on the efficiency of finite volume methods for the incompressible Navier-Stokes equations. In Proceedings of the ASME Fluid Engineering Division Conference; ASME Press: New York, 1996; Vol. 238, N. 3, p 151. (96) Lindborg, H.; Eide, V.; Unger, S.; Henriksen, S. T.; Jakobsen, H. A. Parallelization and performance optimization of a dynamic PDE ractor model for practical applications. Comput. Chem. Eng. 2004, 28, 1585. (97) Caia, C.; Minev, P. A finite element method for an averaged multiphase flow model. Int. J. Comput. Fluid Dyn. 2004, 18, 111. (98) Jakobsen, H. A.; Bourg, I.; Hjarbo, K. W.; Svendsen, H. F. Interaction Between Reaction Kinetics and Flow Structure in Bubble Column Reactors. In Parallel Computational Fluid DynamicssTrends and Applications; Periaux, J., Jenssen, C. B., Pettersen, B., Ecer, A., Eds.; Elsevier Science B. V.: Amsterdam, 2001; pp 543-550. (99) Jakobsen, H. A.; Svendsen, H. F.; Hjarbo, K. W. On the prediction of local flow structures in internal loop and bubble column reactors using a two-fluid model. Comput. Chem. Eng. 1993, 17S, S531. (100) Jakobsen, H. A. General REActor Technology FUNdamentals. In Lecture Notes in Subject KP8103 Advanced Reactor Modeling; NTNU: Trondheim, Norway, 2004; p 890. (101) Fan, R.; Marchisio, D. L.; Fox, R. O. Application of the direct quadrature method of moments to polydisperse gas-solid fluidized beds. Powder Technol. 2004, 139, 7. (102) Tomiyama, A.; Miyoshi, K.; Tamai, H.; Zun, I.; Sakaguchi, T. A Bubble Tracking Method for the Prediction of SpatialEvolution of Bubble Flow in a Vertical Pipe. Presented at the Third International Conference on Multiphase Flow, ICMF′98, Lyon, France, Jun 8-12, 1998. (103) Ishii, M.; Mashima, K. Two Fluid Model and Hydrodynamic Constitutive Relations. Nucl. Eng. Des. 1984, 82, 107. (104) Drew, D. A.; Lahey, R. T. Application of General Constitutive Principles to the Derivation of Multidimensional Two-Phase Flow Equations. Int. J. Multiphase Flow 1979, 7, 243. (105) Antal, S. P.; Lahey, R. T.; Flaherty, J. E. Analysis of Phase Distribution in Fully Developed Laminar Bubbly Two-Phase Flow. Int. J. Multiphase Flow 1991, 17, 635. (106) Marie´, J. L.; Moussali, E.; Lance, M. A First Investigation of a Bubbly Stress Measurements. Presented at the International Symposium on Turbulence Modification in Multiphase Flows, ASME/JSME Annual Meeting, Portland, OR, Jun 23-27, 1991. (107) Jakobsen, H. A.; Lindborg, H.; Handeland, V. A numerical study of the interactions between viscous flow, transport and kinetics in fixed bed reactors. Comput. Chem. Eng. 2002, 26, 333. (108) van Leer, B. Towards the Ultimate Conservation Difference Scheme 2. Monotonicity and Conservation Combined in a Second-Order Scheme. J. Comput. Phys. 1974, 14, 361. (109) van Leer, B. Towards the Ultimate Conservation Difference Scheme 4. A New Approach to Numerical Convection. J. Comput. Phys. 1977, 23, 276. (110) Grienberger, J. Untersuchung und Modellierung von Blasensaülen. Dr. Ing. Thesis, Der Technischen Fakulta¨t der Universita¨t Erlangen-Nu¨rnberg, Erlangen, Germany, 1992. (111) Kurul, N.; Podowski, M. Z. On the Modeling of Multidimensional Effects in Bioling Channels. In Proceedings of the 1991 National Heat Transfer Conference; ASME Press: New York, 1991; Vol. 5, pp 30-40. (112) Morsi, S. A.; Alexander, A. J. An Investigation of Particle Trajectories in Two-Phase Flow Systems. J. Fluid. Mech. 1972, 55, 193. (113) Schiller, L.; Naumann, Z. U ¨ ber die grundlegenden Berechnungen bei der Schwerkraftaufbereitung. Z. Ver. Deutsch. Inc. 1935, 77, 318. (114) Schumann, U. Simulations and parametrizations of large eddies in convective atmospheric boundary layers. In Proceedings of a workshop held at the European Centre for medium-range

weather forecasts (ECMWF) on fine scale modelling and the development of parametrization schemes; ECMWF: London, 1992; pp 21-51. (115) Menzel, T. Die Reynolds-Schubspannung als wesentlicher Parameter zur Modellierung der Stro¨ mungssstruktur in Blasensa¨ ulen and Airlift-Schlaufenreaktoren; VDI Verlag: Du¨sseldorf, Germany, 1990; VDI Reihe 3, Nr. 198. (116) Schlichting, H.; Gersten, K. Boundary-Layer Theory; Springer-Verlag: Berlin, 2000. (117) Telionis, D. P. Unsteady Viscous Flows; SpringerVerlag: New York, 1981. (118) Johanson, S. H.; Davidson, L.; Olsson, E. Numerical simulation of vortex shedding past triangular cylinders at high Reynolds number using a k- turbulence model. Int. J. Numer. Methods Fluids 1993, 16, 859. (119) Williams, F. A. Combustion Theory, 2nd ed.; AddisonWesley Publishing Company: Redwood City, CA, 1985. (120) Gidaspow, D. Multiphase Flow and Fluidization-Continuum and Kinetic Theory Descriptions; Academic Press, Harcourt Brace & Company, Publishers: Boston, MA, 1994. (121) Friedlander, S. K. Smoke, Dust, and Haze: Fundamentals of Aerosol Dynamics, 2nd ed.; Oxford University Press: New York, 2000. (122) Hulburt, H. M.; Katz, S. Some problems in particle technology: A statistical mechanical formulation. Chem. Eng. Sci. 1964, 19, 555. (123) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; John Wiley & Sons: New York, 1960. (124) Randolph, A. D. A Population Balance for Countable Entities. Can. J. Chem. Eng. 1964, 42, 280. (125) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes: Analysis and Techniques of Continuous Crystallization, 2nd ed.; Academic Press, Harcourt Brace Jovanovich, Publishers: San Diego, CA, 1988. (126) Ramkrishna, D. The Status of Population Balances. Rev. Chem. Eng. 1985, 3, 49. (127) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering; Academic Press: San Diego, CA, 2000. (128) Williams, M. M. R.; Loyalka, S. K. Aerosol Science Theory and Practice: With Special Applications to the Nuclear Industry; Pergamon Press: Oxford, U.K., 1991. (129) Kolev, N. I. Fragmentation and Coalescence Dynamics in Multiphase Flows. Exp. Therm. Fluid Sci. 1993, 6, 211. (130) Kolev, N. I. Multiphase Flow Dynamics 2: Mechanical and Thermal Interactions; Springer: Berlin, Germany, 2002. (131) Guido, Lavalle G.; Carrica, P. M.; Clausse, A.; Qazi, M. K. A bubble number density constitutive equation. Nucl. Eng. Des. 1994, 152, 213. (132) Pilon, L.; Fedorov, A. G.; Ramkrishna, D.; Viskanta, R. Bubble transport in three-dimensional laminar gravity-driven flowsmathematical formulation. J. Non-Cryst. Solids 2004, 336, 71. (133) Prasher, C. L. Crushing and Grinding Process Handbook; John Wiley & Sons Limited; Chichester, U.K., 1987. (134) Dorao, C.; Jakobsen, H. A. The Kernel Effect in the Population Balance Equation. Poster Presented at CHISA 2004, Prague, Czech Republic, Aug 22-26, 2004. (135) Coulaloglou, C. A.; Tavlarides, L. L. Description of Interaction Processes in Agitated Liquid-Liquid Dispersions. Chem. Eng. Sci. 1977, 32, 1289. (136) Lee, C.-K.; Erickson, L. E.; Glasgow, L. A. Bubble Breakup and Coalescence in Turbulent Gas-Liquid Dispersions. Chem. Eng. Commun. 1987, 59, 65. (137) Wang, T.; Wang, J.; Jin, Y. A novel theoretical breakup kernel function for bubbles/droplets in a turbulent flow. Chem. Eng. Sci. 2003, 58, 4629. (138) Venneker, B. C. H.; Derksen, J. J.; van den Akker, H. E. A. Population Balance Modeling of Aerated Stirred Vessels Based on CFD. AIChE J. 2002, 48, 673. (139) Fleischer, C.; Bierdel, M.; Eigenberger, G. Prediction of Bubble Size Distributions in G/L-Contactors with Population Balances. Poster Presented at the Third German/Japanese Symposium on Bubble Columns, Schwerte, Germany, Jun 13-15, 1994. (140) Kamp, A. M.; Chesters, A. K.; Colin, C.; Fabre, J. Bubble coalescence in turbulent flows: A mechanistic model for turbulenceinduced coalescence applied to microgravity bubbly pipe flow. Int. J. Multiphase Flow 2001, 27, 1363.

Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5151 (141) Smoluchowski, M. Versuch einer mathematischen theorie der koagulationskinetik kolloider lo¨sungen. Z. Physik Chem. 1918, XCII, 129. (142) Kolmogorov, A. N. Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds Number. Dokl. Akad. Nauk. SSSR 1941, 30, 229. (143) Pope, S. B. Turbulent Flows; Cambridge University Press: Cambridge, U.K., 2001. (144) Batchelor, G. K. The Theory of Homogeneous Turbulence; Cambridge University Press: Cambridge, U.K., 1982. (145) Risso, F.; Fabre, J. Oscillations and breakup of a bubble immersed in a turbulent field. J. Fluid. Mech. 1998, 372, 323. (146) Hagesaether, L. Coalescence and Break-up of Drops and Bubbles; Dr. Ing. Thesis, Department of Chemical Engineering, The Norwegian University of Science and Technology, Trondheim, Norway, 2002. (147) Oolman, T. O.; Blanch, H. W. Bubble Coalescence in Stagnant Liquid. Chem. Eng. Commun. 1986, 43, 237. (148) Chester, A. K. The modelling of coalescence processes in fluid-fluid dispersions: A review of current understanding. Trans. Inst. Chem. Eng. 1991, 69, 259. (149) Kuboi, R.; Komasawa, I.; Otake, T. Behavior of dispersed particles in turbulent liquid flow. J. Chem. Eng. Jpn. 1972, 5, 349. (150) Kuboi, R.; Komasawa, I.; Otake, T. Collision and coalescence of dispersed drops in turbulent liquid flow. J. Chem. Eng. Jpn. 1972, 5, 423. (151) Brenn, G.; Braeske, H.; Durst, F. Investigation of the unsteady two-phase flow with small bubbles in a model bubble column using phase-Doppler anemometry. Chem. Eng. Sci. 2002, 57, 5143. (152) Klaseboer, E.; Chevallier, J. Ph.; Masbernat, O.; Gourdon, C. Drainage of the liquid film between drops colliding at constant approach velocity. Paper Presented at the Third International Conference on Multiphase Flow, ICMF98, Lyon, France, Jun 8-12, 1998. (153) Klaseboer, E.; Chevallier, J. Ph.; Gourdon, C.; Masbernat, O. Film drainage between colliding drops at constant approach velocity: Experimental and modeling. J. Colloid Interface Sci. 2000, 229, 274. (154) Orme, M. Experiments on Droplet Collisions, Bounce, Coalescence and Disruption. Prog. Energy Combust. Sci. 1997, 23, 65. (155) Havelka, P.; Gotaas, C.; Jakobsen, H. A.; Svendsen, H. F. Droplet Formation and Interactions under Normal and High pressures. Paper Presented at the 5th International Conference on Multiphase Flow, ICMF′04, Yokohama, Japan, May 30-Jun 4, 2004. (156) Saboni, A.; Gourdon, C.; Chesters, A. K. The influence of inter-phase mass transfer on the drainage of partially mobile liquid films between drops undergoing a constant interaction force. Chem. Eng. Sci. 1999, 54, 461. (157) Saboni, A.; Alexandrova, S.; Gourdon, C.; Chesters, A. K. Inter-drop coalescence with mass transfer: Comparison of the approximate drainage models with numerical results. Chem. Eng. J. 2002, 88, 127. (158) Lasheras, J. C.; Estwood, C.; Martıńez-Bazań, C.; Montan˜e´s, J. L. A review of statistical models for the break-up of an immiscible fluid immersed into a fully developed turbulent flow. Int. J. Multiphase Flow 2002, 28, 247. (159) Hinze, J. O. Fundamentals of the Hydrodynamic Mechanism of Splitting in Dispersion Processes. AIChE J. 1955, 1, 289. (160) Drazin, P. G.; Reid, W. H. Hydrodynamic Stability; Cambridge Univerity Press: Cambridge, U.K., 1981. (161) Kolmogorov, A. N. On the Breakage of Drops in a Turbulent Flow. Dokl. Akad. Navk. SSSR 1949, 66, 825. (162) Politano, M. S.; Carrica, P. M.; Balinõ, J. L. About bubble breakup models to predict bubble size distributions in homogeneous flows. Chem. Eng. Commun. 2003, 190, 299. (163) Azbel, D. Two-Phase Flows in Chemical Engineering; Cambridge Univerity Press: Cambridge, U.K., 1981. (164) Azbel, D.; Athanasios, I. L. A Mechanism of Liquid Entrainment. In Handbook of Fluids in Motion; Cheremisinoff, N., Ed.; Ann Arbor Science Publishers: Ann Arbor, MI, 1983; p 473. (165) Martıńez-Bazań, C.; Montan˜e´s, J. L.; Lasheras, J. C. On the breakup of an air bubble injected into a fully developed turbulent flow. Part 1. Breakup frequency. J. Fluid. Mech. 1999, 401, 157.

(166) Batchelor, G. K. The Theory of Homogeneous Turbulence; Cambridge University Press: Cambridge, U.K., 1956. (167) Batchelor, G. K. The Theory of Homogeneous Turbulence; Cambridge University Press: Cambridge, U.K., 1953. (168) Angelidou, C.; Psimopoulos, M.; Jameson, G. J. Size Distribution Functions of Dispersions. Chem. Eng. Sci. 1979, 34, 671. (169) Tsouris, C.; Tavlarides, L. L. Breakage and Coalescence Models for Drops in Turbulent Dispersions. AIChE J. 1994, 40, 395. (170) Soo, S. L. Particulates and Continuum; Hemisphere Publishing Corporation: New York, 1989. (171) Marrucci, G. A Theory of Coalescence. Chem. Eng. Sci. 1969, 24, 975. (172) Thomas, R. M. Bubble Coalescence in Turbulent Flows. Int. J. Multiphase Flow 1981, 7, 709. (173) Smoluchowski, M. Drei vortrage uber diffusion, Brownsche molekularbewegung und koagulation von kolloidteilchen. Phys. Z. 1916, 17, 557. (174) Smoluchowski, M. Versuch einer mathematischen theorie der koagulationskinetik kolloider lo¨sungen. Z. Phys. Chem. 1917, 92, 129. (175) Ross, S. L.; Curl, R. L. Measurement and models of the dispersed phase mixing process. In Joint Chemical Engineering Conference; Symposium Series 139; AIChE Press: New York, 1973; Paper 29b. (176) Martıńez-Bazań, C.; Montan˜e´s, J. L.; Lasheras, J. C. On the breakup of an air bubble injected into a fully developed turbulent flow. Part 2. Size PDF of the resulting daughter bubbles. J. Fluid. Mech. 1999, 401, 183. (177) Rodrı´guez-Rodrı´guez, J.; Martıńez-Bazań, C.; Montan˜es, J. L. A novel particle tracking and break-up detection algorithm: application to the turbulent break-up of bubbles. Meas. Sci. Technol. 2003, 14, 1328. (178) Eastwood, C. D.; Armi, L.; Lasheras, J. C. The breakup of immersible fluids in turbulent flows. J. Fluid Mech. 2004, 502, 309. (179) Hesketh, R. P.; Etchells, A. W.; Russell, T. W. F. Experimental observations of bubble breakage in turbulent flows. Ind. Eng. Chem. Res. 1991, 30, 845. (180) Hesketh, R. P.; Etchells, A. W.; Russell, T. W. F. Bubble breakage in pipeline flows. Chem. Eng. Sci. 1991, 46, 1. (181) Valentas, K. J.; Amundson, N. R. Breakage and Coalescence in Dispersed Phase Systems. Ind. Eng. Chem. Fundam. 1966, 5, 271. (182) Valentas, K. J.; Amundson, N. R. Breakage and Coalescence in Dispersed Phase Systems. Ind. Eng. Chem. Fundam. 1966, 5, 533. (183) Konno, M.; Aoki, M.; Saito, S. Simulations model for break-up process in an agitated tank. J. Chem. Eng. Jpn. 1980, 13, 67. (184) Hounslow, M. J.; Ryall, R. L.; Marshall, V. R. A Discretized Population Balance for Nucleation, Growth, and Aggregation. AIChE J. 1988, 34, 1821. (185) McGraw, R. Description of Aerosol Dynamics by the Quadrature Method of Moments. Aerosol Sci. Technol. 1997, 27, 255. (186) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. Quadrature method of moments for aggregation-breakage processes. J. Colloid Interface Sci. 2003, 258, 322. (187) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. Implementation of the quadrature method of moments in CFD codes for aggregation-breakage problems. Chem. Eng. Sci. 2003, 58, 3337. (188) Press: W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran 77. The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1992; Vol. 1.

Received for review June 24, 2004 Revised manuscript received September 16, 2004 Accepted September 30, 2004 IE049447X

Modeling of Bubble Column Reactors: Progress and Limitations

Recommend Documents