Numerical Simulation of Bubble Column Flows in Churn-Turbulent

Article pubs.acs.org/IECR

Numerical Simulation of Bubble Column Flows in Churn-Turbulent Regime: Comparison of Bubble Size Models Lijia Xu, Boruo Yuan, Haoyin Ni, and Caixia Chen* Key Laboratory of Coal Gasification and Energy Chemical Engineering of Ministry of Education, East China University of Science and Technology, Shanghai, China ABSTRACT: Numerical simulations of cylindrical bubble column operating in the churn-turbulent regime have been simulated using Euler−Euler approach incorporated with the RNG k−ε model for liquid turbulence. Single-size bubble model, double-size bubble model, and the multiple size group model (MUSIG), including the homogeneous and inhomogeneous discrete methods, are employed in the simulations. Mass conserved formulations of breakup and coalescence rates were used in the computation of bubble size distributions. The Schiller−Naumann drag force was used in the single-size model, and the Ishii−Zuber drag force was used for the MUSIG simulations. For the double-size bubble model, an empirical drag formulation was adapted. The simulation results of time-averaged axial velocity and gas holdup obtained with the three models were compared with reported experimental data in the literature. The results showed that only MUSIG models with lift force can reproduce the measured radial distribution of gas holdup in the fully developed flow regime and that the inhomogeneous MUSIG model performs a little better than other models in the prediction of axial liquid velocity. The RNG k−ε model was used in all simulations, and the results confirmed that this version of k−ε model did yield relatively high turbulence dissipation rates and high bubble breakup rates and, thus, resulted in a rational bubble size distribution. The ad hoc manipulation of the breakup rates was avoided. The simulation results indicated profound mutual effects of drag force, mean bubble sizes, and turbulence characteristics. An increase in drag force yielded a decrease in the relative velocity between phases, the later could result in decreases in k and ε. A large Sauter diameter results from a low bubble breakup rate which was directly connected to the dissipation rates of turbulence. The change of Sauter diameter, in turn, influenced the drag force.

1. INTRODUCTION Gas−liquid bubble columns are widely used in the petrochemical and biochemical processes due to their numerous advantages of simple structures, easy operation, and good heat and mass transfer characteristics.1 Depending on the superficial gas velocity used in the reactors, the flow in the buoyancy-driven gas−liquid bubble column is represented by two regimes: the well-known bubbly flow regime and the churn-turbulent flow regime. In the turbulent gas liquid flows, the flow structures can become very complex due to substantial coalescence and breakup of bubbles. Therefore, a thorough understanding of the two-phase flow characteristics of the bubble column is very important for a quantitative description of the heat and mass transfer in reactors operating in the churnturbulent regime. Computational fluid dynamics (CFD) has provided a useful technique for the description of fluid dynamics of bubbling flows (e.g., Chen and Fan,2 Buwa et al.,3 Hu and Celik4). However, CFD simulation of multiphase physics of fluids is not yet as satisfactory as single-phase flow. The Eulerian multifluid model has been identified as a practical approach to resolving the flow structures and used in a variety of numerical simulations of bubble column reactors (e.g., the work of Jakobsen et al.,5 Joshi,6 and Rampure et al.7). Nevertheless, a computational method for the gas−liquid turbulent flows is not well-established, mainly due to the uncertainties related to closure models of gas−liquid interactions. These closure models are (1) description of bubble induced turbulence of liquid phase, (2) definition of interaction forces between gas and liquid phases, (3) determination of local bubble size © 2013 American Chemical Society

distribution, which is closely connected to the modeling of both the turbulence and the interaction forces. The interfacial forces are normally modeled using the relative velocity between the liquid and gas phases. Although many effects, such as lift force, added mass force, and turbulent dispersion force, could contribute to the interfacial momentum transfer, it is generally agreed that drag force is more important than other contributors (Oey et al.,8 Olmos et al,9 Chen et al.10). There are various empirical models proposed for the computation of drag force coefficient in the literature (e.g., the work of Schiller and Naumann,11 Grace et al.,12 Ishii and Zuber,13 Ma and Ahmadi,14 and Zhang and VanderHeyden15). Each model is suitable for limited cases of flow conditions. For example, the Schiller and Naumaan model is typically used for small, spherical bubbles, and the Ishii and Zuber model is formulated for the cases of flow systems with a wide range of bubble size distribution. It is argued that the lift force is insignificant compared to the drag force; hence, the lift force effect is frequently neglected in the computations. The classical lift force formulation for gas−liquid flows uses a constant positive coefficient in the range of 0−0.5 and acts in the direction of decreasing liquid velocity. Tabib et al.17 suggested the lift coefficient was not a constant but bubble size dependent. In a inhomogeneous bubble column, Lucas et al.16 observed that the small bubbles tend to move toward the Received: September 24, 2012 Accepted: April 22, 2013 Published: April 22, 2013 6794

dx.doi.org/10.1021/ie4005964 | Ind. Eng. Chem. Res. 2013, 52, 6794−6802

Industrial & Engineering Chemistry Research

Article

In the present paper, we report simulation results of bubble columns using the RNG k−ε model for the liquid turbulence. The experimental setups of Chen.30 are used as a benchmark case. The homogeneous and inhomogeneous MUSIG models implemented in the ANSYS Fluent software are used to model the gas bubble phase. At the same time, simulations using single-size and double-size bubble models are performed, respectively. The computed results of different models are compared with each other and the experimental data as well. Various drag force formulations were evaluated computationally in connection to the inhomogeneous MUSIG models, and the interactions between bubble sizes and liquid turbulence characteristics are discussed.

wall and the large bubbles move toward the center of the bubble column in a experimental study. These findings support the lift force model proposed by Tomiyama et al.,18 where the lift force coefficient of small bubbles is positive and that of large bubbles is negative. Determining the local mean bubble sizes and their distribution is also crucial for the modeling of drag force. Traditionally, a single-size bubble model has been used in simulations (e.g., the work of Sanyal et al.,19 Zhang et al.,20 and Laborde-Boutet et al.21). Nevertheless, bubbles with different dimensions have very different rise velocities, leading to substantially different drag forces. A uniform bubble size could make significant error for the simulation of some systems like the churn-turbulent bubble columns where the bubbles are represented in a wide range of sizes. Krishna proposed a double-size bubble model, which treated the gas phase as a mixture of two classes of bubbles, i.e.: “large” and “small” bubbles. Specifically, the small bubbles represent bubbles of 1− 6 mm in diameter, and the large bubbles are in the range 20−80 mm.22 Using this double-size model, Krishna and his group performed a series of simulations and obtained reasonable results.22−24 In the past decade, substantial efforts have been made in the simulation of the multiple size group (MUSIG) model. Taking the bubble coalescence and breakup into account, various versions of bubble coalescence and breakup models have been implemented into computational fluid dynamics (CFD).25−27 Chen et al.28 investigated the sensitivity of different bubble coalescence and breakup models and concluded that the choice of the closure model did not have significant impact on their simulation results. They arguably introduced a 10-fold increase of bubble break-up rate in order to obtain a reasonable bubble size distribution. Later, Laborde-Boutet et al.21 pointed out that the underestimation of ε by the standard k−ε model might be a major reason that causes the imbalance of bubble coalescence and breakup. They validated different k−ε formulations by comparing the RANS simulation with experimental data of Chen.30 The results clearly indicate that only an RNG k−ε model was able to reasonably capture the flow characteristics in bubble columns operating in the churn-turbulent regime. Laborde-Boutet et al.21 further presented results showing that a use of the RNG k−ε model provided the highest turbulent dissipation rate than other k−ε versions like the standard and realizable models. Because in most models, the bubble break-up rates are highly depending on the local turbulent dissipation rate, ε, this encouraging finding hence implied that the RNG k−ε model can potentially improve the bubble breakup modeling, especially the “try-and-error” 10-fold manipulation of bubble breakup rates can be avoided when a bubble coalescence and breakup model was implemented into CFD simulation.21 In the homogeneous MUSIG model, bubbles with different sizes are assigned to the same gas phase and therefore move with the same velocity. In the inhomogeneous MUSIG model, the gas phase is divided into two or more gas phases, which allows different groups of bubbles to be adverted by different phase velocities. Krepper et al.29 used the inhomogeneous MUSIG model to simulate the upward gas−liquid flow in vertical pipes. Good agreement was found for gas holdup and velocity profiles against experimental data. In their work, the difference of the homogeneous and inhomogeneous MUSIG models was not mentioned.

2. DESCRIPTION OF COMPUTATIONAL METHODS 2.1. Governing Equations. For the gas−liquid flow system, liquid is treated as a continuous phase, and the gas bubble is a dispersed secondary phase. The continuity and momentum equations of two fluids are represented in the Eulerian−Eulerian formulation as follows: ∂ρi αi + ∇·(αi u i) = 0 (1) ∂t ∂(αiρi u i) ∂t

+ ∇·(αiρi u iu i)

= −αi·∇p − ∇(αiτi) + Fi , j + αiρi gαi

(2)

Where, αi is the volume fraction of the ith phase. The sum of the volume fractions of all phases is equal to 1. The terms ρi, τi, ui, and Fi,j represent density, viscous stress tensor, velocity, and interface momentum exchange term, respectively. The RNG turbulence model is recommended by Laborde-Boutet et al.,21 and the k and ε of each phase are computed by the following equations: ∂(αiρi ki) ∂t

+ ∇·(αiρi u iki)

⎛ μ ⎞ = ∇·⎜αi i ∇ki⎟ + αiGk , i − αiρi εi + αiρi Πk , i ⎝ σk ⎠ ∂(αiρi εi) ∂t

(3)

+ ∇·(αiρi u iεi)

⎛ μ ⎞ ε = ∇·⎜αi i ∇εi⎟ + αi i (C1Gk , i − C2ρi εi) + αiρi Πε , i σ k ⎝ ε ⎠ i (4)

Where, Πk,i and Πε,i are defined as follows: 3 αjCd|ui , j|(ki , j − 2ki + ui , j ·udr) 4 ε = C 3 i Πk , i ki

Πk , i =

Πε , i

(5)

(6)

Where, σk, C1, C2, and C3 are constants. The term ki,j is the covariance of the velocities of phase i and phase j. Also, ui,j is the relative velocity between the liquid and dispersed phases. Here, Gk,i is the production rate of turbulent energy. 2.2. Bubble Size Models and Drag Law. 2.2.1. SingleSize Bubble Model. The single-size bubble model is represented by a uniform bubble size. The diameter is calculated using Wilkinson’s correlation:31 6795


Industrial & Engineering Chemistry Research dg = 3g −0.44 σ 0.34μ l 0.22 ρl −0.45 ρg −0.11U −0.02

Article

Where, BC, DC, BB, DB represent birth rate and death rate due to coalescence and breakup, respectively. These rates are calculated following the Hageseather method:33

(7)

For our air−water system, the diameter of single bubble is 4.5 mm. The drag force coefficient is determined by Schiller− Naumann’s drag model: ⎧ 24 ⎪ (1 + 0.15Re 0.687) Re ≤ 1000 Cd = ⎨ Re ⎪ ⎩ 0.44 Re > 1000

N

BC, i =

gdb,small 2σ + ρl db,small 2

ub,large = 0.71 gdb,large (SF)(AF)(DF)

∂αlarge ∂t

∑ Si i=1

∑

xi + 1, j Ω B(vi + 1 , vj)

j = 1, i ≠ N

i

∑

+

xi , j Ω B(vi , vj) (18)

j = 1, i ≠ N i−1

DB, i =

(9)

∑ ΩB(vi , vj) (19)

j=1 1+j−i

Here xi,j = 2 , v is the bubble volume, n is the bubble number density function, and the correction factors x and ξ are used to conserve the mass. In the present work, the breakup rate ΩB(vi,vj) proposed by Luo and Svendsen,27 and the coalescence rate cij of Luo34 are used, respectively. The formulations are briefly described as follows. The correction factors xkj and ξkj are defined as follows:

(10) (11)

⎧1 Vi − 1 < Vag < Vi + 1 ξkj = ⎨ ⎩ 0 otherwise

(20)

⎧ (Vag − Vi − 1) ⎪ Vi − 1 < Vag ≤ Vi ⎪ (Vi − Vi − 1) ⎪ ⎪ (Vag − Vi + 1) xkj = ⎨ Vi < Vag < Vi + 1 ⎪ (Vi − Vi + 1) ⎪ ⎪ Vag VN < Vag ⎪V ⎩ N

(21)

Vag is the particle volume resulting from the aggregation of particles k and j. The coalescence rate cij is defined as

The breakup rate ΩB(vi,vj) can be written as ⎛ ε ⎞1/3 Ω B(vi , vj) = 0.9238αlni⎜ 2 ⎟ ⎝ di ⎠

(13)

∫ξ

1

min

(1 + ξ)2 11/3

ξ

exp( −bξ −11/3) dξ

b = 12(f 2/3 + (1 − f )2/3 − 1)σρ−1ε−2/3di−5/3

(14)

Where, Si is the source term of ith bubble group generated by bubble breakage and coalescence, which is expressed as Si = BC, i − DC, i + BB, i − DB, i Vi

Ω B(vj , vi) +

j = i + 1, i ≠ N

Nlarge

+ ∇·(αlargeularge) =

i

∑

BB, i =

Nsmall i=1

(17) N

In our simulation, a constant drag coefficient of 0.839 is used for small bubbles of 4 mm in diameter and 0.289 is used for the large bubbles. 2.2.3. MUSIG Model. In the homogeneous MUSIG model, bubbles with different size move with the same velocity. In the inhomogeneous MUSIG model, the gas phase is divided into two or more gas phases with different momentum fields. It is suggested that two or three gas phases are sufficient to capture gas−liquid flow characteristics.29 For the inhomogeneous MUSIG model, gas phase is divided into small bubble and large bubble groups. The continuity equations for the two bubble groups are represented as follows:

∑ Si

∑ cijninj j=1

Where, SF, AF, and DF represent scale correction factor, acceleration factor, and density correction factor, respectively. The drag force coefficients of bubbles are computed by 4 ρl − ρg 1 CD = gdb 2 3 ρl ub (12)

∂αsmall + ∇·(αsmallusmall) = ∂t

(16)

N

DC, i =

(8)

Where, Utrans is superficial gas velocity at regime transition. For an air−water system, Utrans takes a value of 0.034 m/s. When U = 0.1 m/s, dlarge = 25 mm. Krishna22 introduced different methods for the computation of rise velocities for small and large bubbles, respectively: ub,small =

∑ ∑ ckjnknjxkjξkj k=1 j=1

2.2.2. Double-Size Bubble Model. The double-size model is a self-contained bubble size model that divided the bubbles into two classes: small bubbles and large bubbles. According to Harmathy,32 the rise velocity of small bubbles is independent of bubble sizes, thus it is not necessary to provide a specific diameter for the small bubbles. Small bubbles are in the range 1−6 mm, and in this case, the diameter of small bubbles is set to 4 mm. The diameter of large bubbles is calculated by dlarge = 0.069(U − Utrans)0.376

N

(23) (24)

In the breakup rate of eq 24, f is the volume fraction of one daughter bubble, and ξ is the ratio of eddy size to parent bubble. In eq 22, tdrainage and tcontact represent film drainage time and bubble contact time, respectively. C0 is an adjustable parameter and took very different values in the existing collision

(15) 6796



Article

frequency models. For example, in the models of Lee et al.,35 Prince and Blanch,25 and Luo,34 the values of C0 are 1.4142π/4, 0.28−1.11, and 1.43π/4, respectively. In our computation, the bubble coalescence rate has been reduced by setting C0 = 0.5. When the MUSIG model is used, the Sauter mean diameters of bubbles change in a wide range. The Ishii−Zuber model (nonspherical bubble drag model) is used for the computation of the drag force coefficient: Cd,spherical =

24 (1 + 0.1Re 0.75) Re

(25)

2 g (ρl − ρg ) ⎛ 1 + 17.67f 6/7 ⎞ 2 α ⎜ ⎟ = d ⎜ ⎟ , 18.67fα σ 3 ⎝ ⎠

Cd,distorted

fα = αl1.5

Cd,cap =

is 0.1 m/s. No-slip wall boundary condition is used for the side wall and a pressure outlet is used for the top of the column. For bubbly flow, 3D time dependent simulation is recommended. The RNG k−ε formulation is used for the modeling of bubble induced turbulence as recommended by Laborde-Boutet et al.21 The QUICK discretization scheme is used for volume fraction equation, and the second order upwind scheme is employed to discretize other equations. In order to capture the gradients of flow field near the wall, the near wall grid is refined with four rows of thinner cells. Hexahedral, structured computational cells are employed in computations. The total grid number is 67 392 for all cases except the case simulated by two-size model. For case 4, a special treatment of the inlet condition brings the grid number to 98 496. A time step of 0.01 s recommended by Chen et al.10 is used. The simulation results are averaged over 80 s after the flow pattern reaches a quasi-steady state.

(26)

8 2 al 3

3. RESULTS AND DISCUSSION 3.1. Single-Size Bubble Model. Simulations were performed using the single-size bubble model. A uniform bubble size of 4.5 mm was used in the computations. Two different sizes, 2.5 and 6.5 mm, were also used for comparison. Figure 1 compares the instaneous distributions of gas holdup

(27)

If Cd,spherical > Cd,distorted , Cd = Cd,spherical else Cd = min(Cd,distorted , Cd,cap)

(28)

For comparison, the Schiller−Naumann model of eq 8 is also used in the simulations. The lift force is taken into account and the Tamiyama lift force18 can be calculated as Flift = −C liftρl αg(ul − ug) × (∇ × ul)

(29)

Where

C lift

⎧ min[0.288tanh(0.121Re) Eo < 4 ⎪ ⎪ , f (Eo)] =⎨ ⎪ f (Eo) 4 < Eo < 10 ⎪ ⎩−0.27 10 < Eo

(30)

f (Eo) = 0.00105Eo3 − 0.0159Eo2 − 0.0204Eo3 + 0.474 (31)

Here Eo is the Eötvös number. The Sauter mean bubble diameter, which is used to calculate the drag force and the lift force, is calculated by N

db =

∑i =gas1 nidi 3 N

∑i =gas1 nidi 2

(32)

For the inhomogeneous MUSIG model, the overall Sauter mean bubble diameter is calculated in the postprocessing: agas asmall + alarge db = π N = a alarge 2 small + large ∑i =small nidi + d 1 6 d small

large

Figure 1. Distribution of instant gas holdup simulated with single-size bubble model using different bubble sizes Db, (a) Db = 2.5 mm; (b) Db = 4.5 mm; (c) Db = 6.5 mm.

(33)

2.3. Numerical Scheme and Benchmark Cases. Simulations are performed based on the experimental setup of Chen.30 The superficial gas velocity in Chen’s experimental study changed in range 0.02−0.1 m/s, covering the homogeneous and inhomogeneous bubble regimes. In the present work only the turbulent condition at superficial gas velocity 0.1 m/s is used as a benchmark case. In the computational cases, the diameter and height of air−water cylindrical bubble column are 0.44 and 2.43 m, respectively. Initially, the static liquid height of the column is 1.76 m. The gas holdup is 1 above this level, and the superficial gas velocity

for three bubble sizes. For the three simulation cases, the height of the bubble column is 2.43 m, and an initial liquid bed height of 1.76 m is used. Figure 1 clearly shows that the bed expansion increases with decreasing bubble size, indicating a lower averaged gas holdup is obtained when a larger bubble size is presented in the system. In order to compare the simulated results with experimental data of Chen,30 the time averaged simulation results are further circum-averaged at the plane z = 1.32 m. Specifically, the radial 6797



Article

distributions of the axial liquid velocity and the gas holdup were obtained by averaging data sampled from eight radial lines uniformly distributed in the circumferential direction. Figure 2a

Figure 3. Inlet boundary condition setup for double-size model simulation.

Figure 2. Time-averaged radial distribution of axial liquid velocities and gas holdup simulated with single-size model using different bubble sizes Db= (a) 2.5; (b) 4.5; (c) 6.5 mm.

and b show the radial axial liquid velocity and gas holdup. It is found that the axial liquid velocities are sensitive to bubble size. The results using a bubble diameter of 4.5 mm agree with the experiments better than a larger size and a smaller size of bubbles. As shown in Figure 2b, the gas holdup increases drastically with decreasing bubble size. Comparing with the experimental data, only a bubble size of 4.5 mm can predict a reasonable level of gas holdup, especially in the column center. For the three cases, the gas holdup decreases to zero on the wall due the use of the no slip boundary condition. Except the near wall region, the gas holdup of three bubble sizes are uniformly distributed along the radial direction, which is not inconsistent with the experiment, where a gradual decrease is presented from about half of a column diameter. 3.2. Double-Size Bubble Model. The simulations were also performed by using the double-size bubble model. In the simulation, the inlet boundary conditions used in Krishna’s work were adapted. That is, the small and large bubbles are introduced into the system through different inlets as shown in Figure 3. The small bubbles are sprayed from inlet-A and inletB, while the large bubbles are only sprayed from inlet-A. By doing so, the large bubbles are restricted in the inner circle area and move rapidly along the centerline. The small bubbles are introduced in the outer circle uniformly. The instaneous gas holdup of small and large bubbles is shown in Figure 4a and b, respectively. The sum up of the two gas holdups is shown in Figure 4c. Comparing with the results of single-size bubble

Figure 4. Distribution of instant gas holdup simulated with double-size bubble model: (a) small bubbles; (b) large bubbles; (c) total gas holdup.

model as shown in Figure 1, the flow structure computed from the double-size model is simpler, featuring in less fine structures of eddies. The flow pattern is dominantly controlled by the large bubbles. The time averaged liquid axial velocity and the gas holdup are shown in Figure 5a and b. The same experiment data shown in Figure 2 are compared. In Figure 5a, the gas velocities of the large and small bubbles are also presented. Note that the gas velocity of the small bubbles show a slightly negative value in the near wall region, indicating that the small gas bubbles are entrained by the liquid phase when a down flow is formed in response to the strong up flow of liquid in the column center. From Figure 5a, the axial liquid velocity agrees with the experiments relatively well; the radial distribution of gas holdup, however, shows significant differences between the simulation and the experimental data as indicated in Figure 5b. The simulated total gas holdup decreases abruptly below the measured points just away from the centerline. Ignoring the mass and momentum transfer between the two bubble classes may contribute to the discrepancy. Artificial treatments of the large and small bubble sizes without considering the median size bubbles may also pose a problem. 6798



Article

Figure 5. Time-averaged radial distribution of axial velocities of gas and liquid and gas holdup simulated with the double-size bubble model.

Figure 6. Comparison of time-averaged radial distribution of axial liquid velocities and gas holdup simulated with single-size model, double-size model, and MUSIG models: (black squares) exp; (green dotted line) homogeneous; (green solid line) homogeneous + lift; (red dotted line) inhomogeneous; (red solid line) inhomogenous + lift; (blue dotted line) single-size model Db = 4.5 mm; (black dotted line) double-size model.

3.3. MUSIG Model. The computed distributions of timeaveraged liquid velocities and the gas holdups obtained with homogeneous and inhomogeneous MUSIG models are shown in Figure 6a and b, respectively. For all MUSIG runs, the Ishii− Zuber drag force formulation was used. Parallel runs were performed including and excluding the lift force term as indicated in the figures. The simulation results obtained with the single-size and the double-size bubble model, and the experimental data as shown in Figures 2a and b and 5a and b are redrawn in the same picture. It is noted that the predicted liquid velocities based on homogeneous and inhomogeneous MUSIG models are relatively close. When the lift force was excluded in the simulations, the predicted profiles of gas holdup show large discrepancies compared to the experiment data. The homogeneous and inhomogeneous MUSIG models without lift force follow the single-size model very closely, while the homogeneous and inhomogeneous MUSIG models with lift force show a significant improvement. The gas holdup distribution calculated by the MUSIG models with lift force decreases more steeply than that of the MUSIG models without lift force. Both the homogeneous and inhomogeneous MUSIG models agree with the measured gas holdups in the fully developed flow region. Comparing the three bubble size models, i.e., the single-size, the double-size, and the MUSIG models, the single-size model predicts lower liquid velocities in the center area and very high velocities near the wall. On the other side, the double-size model predicts much higher liquid velocities in the center, and much lower liquid velocities in the wall region. Observing the simulated gas holdups of the three models as shown in Figure 6b, the results obtained with MUSIG models show the same pattern as experimental data, as long as the lift force is included

in the computation. While the effect of lift force is considered, the inhomogeneous MUSIG model perform a little better than the homogeneous MUSIG model. The impacts of different drag force models are investigated in combination with inhomogeneous MUSIG model. Here two drag force models, i.e., Schiller−Naumann and Ishii−Zuber models are employed in the computation of the drag forces while keep other setups the same. The simulated radial distributions of the Sauter diameters and the axial gas velocities are shown in Figure 7a and b. It is interesting to find that the simulated bubble size of Ishii−Zuber model is significantly larger than the Schiller−Naumann model. The gas phase velocity of the Ishii−Zuber model, on the contrary, is slightly lower than that of the Schiller−Naumann model. As is wellknown, the rise velocities of large bubbles are usually higher than that of small bubbles. The simulated lower mass-averaged gas velocity is due to a significantly larger drag force computed with the Ishii−Zuber model. It is therefore concluded that a suitable drag force model is critically necessary for a reasonable simulation of the dynamic flow field when the inhomogeneous MUSIG model is incorporated to the CFD simulation. Figure 8a and b compare the turbulent kinetic energy k and turbulent dissipation rate ε computed from the simulations with different drag models, respectively. When the Schiller− Naumann model was used in the simulation, both k and ε are higher than those of the Ishii−Zuber model. Recalling that the production terms in the conservation equations of k and ε, i.e.: 6799



Article

Figure 8. Time-averaged radial profiles of turbulent kinetic energy and turbulent dissipation rate simulated with inhomogeneous MUSIG model using different drag models.

Figure 7. Time-averaged radial distributions of the Sauter diameters and axial gas velocities simulated with inhomogeneous MUSIG model using different drag models.

When the simulations were conducted using the RNG k−ε model in combination with the population balance model as implemented in ANSYS Fluent 14.0, a rational bubble size distribution can be obtained without an ad hoc empiric manipulation of the bubble breakup rate. The simulation results indicated profound mutual effects of drag force, mean bubble sizes, and turbulence characteristics. An increase in drag force yielded a decrease in the relative velocity between phases; the later could result in decreases in k and ε. A small Sauter diameter results from a high bubble breakup rate which was directly connected to the dissipation rate of turbulence. The change of Sauter diameter, in turn, influenced the drag force.

Πk,i and Πε,i in eqs 3 and 4, since Πk,i and Πε,i are proportional to the relative velocity between the liquid and dispersed phase (eqs 5 and 6), the difference in the relative velocity between the phases may have important impacts on the production rates of k and ε. Specifically, for large bubbles the Ishii−Zuber model produces a drag force coefficient of 2.667, which is much higher than the value of 0.44 produced by the Schiller−Naumann model. Therefore, the Ishii−Zuber model produces a smaller relative velocity and leads to the decreases in k and ε. A lower value of ε is directly connected to lower breakup rates of bubbles and a slightly larger Sauter diameter. These results reinforce the importance of accurate simulations of k and ε.

■

AUTHOR INFORMATION

Corresponding Author

4. CONCLUSION We have compared the performance of three bubble size models, namely, a uniform single-size model, a double-size model, and the MUSIG models with experiments for a bubble column operating in the churn-turbulent regime. We found that all of these models obtained a satisfactory prediction of the axial liquid velocity distribution. However, both the single-size model and the double-size model could not reproduce the measured radial distribution of gas holdup in the fully developed flow region. Only by taking the lift force into account, the homogeneous and the inhomogeneous MUSIG models could produce comparable gas holdup as the experimental results. Furthermore, the inhomogeneous MUSIG model was found to perform better than other models in the prediction of gas holdup distribution in general, and the axial liquid velocity near wall region, specifically.

*Tel.: 86-021-64252057. E-mail address: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by China’s National Science Foundation (NSFC) under grants 21276085 and 20876049. C.C. gratefully acknowledges the financial support of the Pujiang Excellent Talents Plan of Shanghai. We wish to express our thanks to Dr. Jasper Van Baten for his help advising us on the simulation of the double-size bubble model.

■ 6800

NOTATION BB = production of bubble by breakup, m−3·kg−1·s−1 DB = destruction of bubble by breakup, m−3·kg−1·s−1 dx.doi.org/10.1021/ie4005964 | Ind. Eng. Chem. Res. 2013, 52, 6794−6802


Article

BC = production of bubble by coalescence, m−3·kg−1·s−1 DC = destruction of bubble by coalescence, m−3·kg−1·s−1 c = coalescence rate, m−3·s−1 Cd = drag coefficient, dimensionless Clift = lift coefficient, dimensionless d = bubble diameter, m FD = drag force, N·m−3 Flift = lift force, N·m−3 Flg = interfacial momentum exchange term, N·m−3 FA = acceleration factor, dimensionless FD = density correction factor, dimensionless FS = scale factor, dimensionless g = gravity, m2·s−1 g(V′) = breakup frequency, s−1 G = production of turbulent energy, W·m−3 k = turbulent kinetic energy, m2·s−2 p = pressure, Pa n = number of particles per unit volume, m−3 t = time, s tdrainage = film drainage time, s tcontact = bubble contact time, s u = velocity vector, m·s−1 udr = drift velocity, m·s−1 ub = bubble rise velocity, m·s−1 U = superficial gas velocity, m·s−1 Utrans = superficial gas velocity at regime transition, m·s−1 Re = Reynolds number, dimensionless v = bubble volume, m3

(8) Oey, R.; Mudde, R.; Van den Akker, H. Sensitivity study on interfacial closure laws in two-fluid bubbly flow simulations. AlChE J. 2003, 49, 1621. (9) Olmos, E.; Gentric, C.; Midoux, N. Numerical description of flow regime transitions in bubble column reactors by a multiple gas phase model. Chem. Eng. Sci. 2003, 58, 2113. (10) Chen, P.; Sanyal, J.; Dudukovic, M. CFD modeling of bubble columns flows: implementation of population balance. Chem. Eng. Sci. 2004, 59, 5201. (11) Schiller, L.; Naumann, Z. Ü ber die grundlegenden Berechnungen bei der Schwerkraft-aufbereitung. Z. Ver. Deutsch. Ing. 1935, 77, 318. (12) Grace, J.; Wairegi, T.; Nguyen, T. Shapes and velocities of single drops and bubbles moving freely through immiscible liquids. Trans. Inst. Chem. Eng. 1976, 54, 167. (13) Ishii, M.; Zuber, N. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AlChE J. 1979, 25, 843. (14) Ma, D.; Ahmadi, G. A thermodynamical formulation for dispersed multiphase turbulent flows-II: Simple shear flows for dense mixtures. Int. J. Multiphase Flow. 1990, 16, 341. (15) Zhang, D. Z.; VanderHeyden, W. B. The effects of mesoscale structures on the macroscopic momentum equations for two-phase flows. Int. J. Multiphase Flow. 2002, 28, 805. (16) Lucas, D.; Krepper, E.; Prasser, H. Evolution of flow patterns, gas fraction profiles and bubble size distributions in gas-liquid flows in vertical tubes. Trans. Inst. Fluid-Flow Machinery. 2003, 112, 37. (17) Tabib, M. V.; Roy, S. A.; Joshi, J. B. CFD simulation of bubble column-an analysis of interphase forces and turbulence models. Chem. Eng. J. 2008, 3, 589. (18) Tomiyama, A.; Sou, I.; Zun, I.; Kanami, N.; Sakaguchi, T. Effects of Eötvös number and dimensionless liquid volumetric flux on lateral motion of a bubble in a laminar duct flow. Adv. Multiphase Flow 1995, 3. (19) Sanyal, J.; Roy, S.; Dudukovic, M. Numerical simulation of gasliquid dynamics in cylindrical bubble column reactors. Chem. Eng. Sci. 1999, 54, 5071. (20) Zhang, D.; Deen, N.; Kuipers, J. Numerical simulation of the dynamic flow behavior in a bubble column: A study of closures for turbulence and interface forces. Chem. Eng. Sci. 2006, 61, 7593. (21) Laborde-Boutet, C.; Larachi, F.; Dromard, N.; Delsart, O.; Schweich, D. CFD simulation of bubble column flows: Investigations on turbulence models in RANS approach. Chem. Eng. Sci. 2009, 64, 4399. (22) Krishna, R.; Van Baten, J.; Urseanu, M. Three-phase Eulerian simulations of bubble column reactors operating in the churnturbulent regime: a scale up strategy. Chem. Eng. Sci. 2000, 55, 3275. (23) Van Baten, J.; Ellenberger, J.; Krishna, R. Scale-up strategy for bubble column slurry reactors using CFD simulations. Catal. Today 2003, 79, 259. (24) Van Baten, J.; Krishna, R. Eulerian simulation strategy for scaling up a bubble column slurry reactor for Fischer−Tropsch synthesis. Ind. Eng. Chem. Res. 2004, 43, 4483. (25) Prince, M. J.; Blanch, H. W. Bubble coalescence and break-up in air-sparged bubble columns. AlChE J. 1990, 36, 1485. (26) Chesters, A. The modelling of coalescence processes in fluidliquid dispersions: a review of current understanding. Chem. Eng. Res. Des. 1991, 69, 259. (27) Luo, H.; Svendsen, H. F. Theoretical model for drop and bubble breakup in turbulent dispersions. AlChE J. 1996, 42, 1225. (28) Chen, P.; Sanyal, J.; Duduković, M. Numerical simulation of bubble columns flows: effect of different breakup and coalescence closures. Chem. Eng. Sci. 2005, 60, 1085. (29) Krepper, E.; Lucas, D.; Frank, T.; Prasser, H. M.; Zwart, P. J. The inhomogeneous MUSIG model for the simulation of polydispersed flows. Nucl. Eng. Des. 2008, 238, 1702. (30) Chen, J.; Li, F.; Degaleesan, S M.; et al. Fluid dynamic parameters in bubble columns with internals. Chem. Eng. Sci. 1999, 54, 2187.

Greek Letters

α = void fraction, dimensionless β(V,V′) = probability density function (pdf), m−3·kg−1 ε = turbulent dissipation rate, m2·s−3 μ = viscosity, kg·m−1·s−1 Πk,i = bubble-induced turbulence term, W·m−3 Πε,i = bubble-induced turbulence term, W·s−1·m−3 ρ = density, kg·m−3 σ = surface tension, N·m−1 ΩB = breakup rate, m−3·s−1

Subscripts

b = bubble index g = gas index l = liquid index i,j = phase index small = small bubble index large = large bubble index

■

REFERENCES

(1) Dudukovic, M. Trends in catalytic reaction engineering. Catal. Today. 1999, 48, 5. (2) Chen, C.; Fan, L. S. Discrete simulation of gas-liquid bubble columns and gas-liquid-solid fluidized beds. AIChE J. 2004, 50, 288. (3) Buwa, V. V.; Deo, D. S.; Ranade, V. V. Eulerian-Lagrangian simulations of unsteady gas-liquid flows in bubble columns. Int. J. Multiphase Flow. 2006, 32, 864. (4) Hu, G.; Celik, I. Eulerian-Lagrangian based large-eddy simulation of a partially aerated flat bubble column. Chem. Eng. Sci. 2008, 63, 253. (5) 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. (6) Joshi, J. Computational ow modelling and design of bubble column reactors. Chem. Eng. Sci. 2001, 56, 5893. (7) Rampure, M.; Mahajani, S.; Ranade, V. CFD Simulation of Bubble Columns: Modeling of Nonuniform Gas Distribution at Sparger. Ind. Eng. Chem. Res. 2009, 48, 8186. 6801



Article

(31) Wilkinson, P. M. Physical aspects and scale-up of high pressure bubble columns. Ph.D. Thesis, University of Groningen, Nederlands, 1991. (32) Harmathy, T. Z. Velocity of large drops and bubbles in media of infinite or restricted extent. AlChE J. 1960, 6, 281. (33) Hagesaether, L.; Jakobson, H. A.; Svendsen, H. F. A model for turbulent binary breakup of dispersed fluid particles. Chem. Eng. Sci. 2002, 57, 3251. (34) Luo, H. Coalescence, breakup and liquid circulation in bubble column reactors. D.Sc. Thesis, Norwegian Institute of Technology, Norway, 1993. (35) Lee, C. H.; Erickson, L.; Glasgow, L. Bubble breakup and coalescence in turbulent gas-liquid dispersions. Chem. Eng. Commun. 1987, 59, 65.

6802


Numerical Simulation of Bubble Column Flows in Churn-Turbulent

Recommend Documents