Simulations of bubble column reactors using a volume of fluid approach: effect of air distributor.

Akhtar, M. Abid ; Tade, Moses ; Pareek, Vishnu 等

INTRODUCTION

Bubble column reactors have a variety of applications in the chemical, biochemical and petrochemical industries due to their relatively simple construction; favourable heat and mass transfer properties and low operating cost (Degaleesan et al., 2001; Deckwer, 1992). In such reactors, the gas rises in the form of bubbles through a continuous liquid phase. These rising bubbles have an important role in the liquid mixing and mass transfer. Based on bubble size distribution, the flow in these systems is characterized by the homogeneous and heterogeneous regimes. The homogeneous regime (also called bubbly flow regime) is obtained at low superficial gas velocities (<5 cm/s). In the homogeneous regime, bubbles have more or less equal size (generally in mm), with a uniform gas holdup in the reactor. Because of their relatively low rise-velocity, the gas throughput with small bubbles is limited. As superficial gas velocity increases, the heterogeneous regime (also known as churnturbulent regime) is obtained, in which coalescence and breakup occurs and large bubbles with sizes (1-10 cm or more), depending on the column diameter, are formed. These large bubbles have much higher rise velocities (up to 1.5 m/s) and are mainly responsible for high throughputs of gas at high velocities (Ranade, 2002; Krishna and van Baten, 2002). The characteristics and hydrodynamics of bubbles, bubble size distribution, superficial gas velocity and gas distributor configuration are a few factors which govern the performance of a bubbling system. Among these, the bubble size is the most important parameter because it not only affects the bubble rise velocity but also has a direct influence on the gas holdup and interfacial area; consequently, this is an important criterion for evaluating the efficiency of a bubble column reactor. Therefore, a study on the effect of bubble size is needed for the better understanding of the gas dispersion mechanisms (Lehr et al., 2002; Krishna and van Baten, 1999). Finally, in view of its importance in the heat and mass transfer rates, it is necessary to thoroughly understand the effect of different parameters (such as superficial gas velocity and distributor configuration) that affect the bubble size (Liu et al., 2005a; Chen et al., 2005; Karamanev, 1994; Clift et al., 1978).

Among the different modelling techniques for bubble column reactor, the Eulerian-Lagrangian model can be used to track the individual bubble using the equation of motion, with external forces accounted for using a set of constitutive equations (Van Wachem and Almstedt, 2003; Lain et al., 1999). The advantage of this approach is that bubble dynamic characteristics can be accessed, including the bubble trajectory and bubble-bubble interaction. However, in this approach, the gas bubbles are assumed to be spherical, which limits its application to the bubble columns with small bubbles in the homogeneous regime. In the heterogeneous regime, the consideration of large, deformable bubble becomes important in order to simulate different aspects of the large bubbles, such as bubble breakage and coalescence (Li et al., 2000). This drawback may be overcome by using the interface tracking methods.

In the interface tracking methods, two different approaches can be used: the front-tracking approach based on the Lagrangian grid system and the volume-tracking approach based on the Eulerian framework. The front-tracking approach constructs a mesh system on the interface to track the movement of bubble surfaces (Van Sint Annaland et al., 2005b; Hirt and Nichols, 1981) and the entire Lagrangian grid mesh system moves with the interface. This approach can achieve a high degree of accuracy but requires substantial CPU time.

The volume of fluid (VOF) (Hirt and Nichols, 1981) is one of the most well known methods for the volume tracking. It relies on the fact that two or more phases are not interpenetrating. In this approach, the motion of all phases is modelled by solving a single set of transport equations using appropriate jump boundary conditions at the interface (Krishna and van Baten, 1999; Tomiyama et al., 1998; Delnoij et al., 1997). The fluid location is recorded by employing a volume of fluid function, or colour function, which is defined as unity within the field regions and zero elsewhere. In practical numerical simulations employing a VOF algorithm, this function is unity in computational cells occupied completely by fluid of phase 1, zero in region occupied completely by phase 2, and a value between these limits in cells which contain a free surface. In this approach, some interface related forces, such as surface or adhesion forces can also be modelled accurately (Van Wachem and Almstedt, 2003). The only drawback of VOF method is the so-called artificial (or numerical) coalescence of gas bubbles which occurs when their mutual distances is less than the size of the computational cell. Therefore, very fine meshes should be used, which makes this approach memory intensive for the simulation of the dispersed multiphase flows in large equipment (Ranade, 2002).

Hydrodynamic modelling of bubble column reactor using the Eulerian-Eulerian (Akhtar et al., 2004; Sokolichin et al., 2004) and Eulerian-Lagrangian approaches (Zhang and Ahmadi, 2005; Chen and Fan, 2004; Li et al., 2000) has been a popular subject. In comparison, VOF has received little attention and most of the available studies deal with a single bubble. This is because the motion of single bubbles is relatively well understood and extensive experimental data on shape and terminal velocity are available in the literature (Clift et al., 1978). Using these experimental data, simulations have been performed for single bubbles rising in stagnant fluid by many researchers mostly in twodimensions (Liu et al., 2005a; Essemiani et al., 2001; Krishna and van Baten, 1999) and some in three-dimensions (Van Sint Annaland et al., 2005a; Olmos et al., 2001). The rise trajectories of bubbles, their size and shape, the rise velocity, the effect of fluid properties on the bubble dynamics, and the gas holdup have been analyzed. The calculation of air bubble terminal velocities and shapes in stagnant water have also been investigated (Krishna et al., 1999; Rudman, 1997). The effect of column diameter on the bubble size was studied by Krishna (2000) but it was limited to the rise of a single bubble. For a precise prediction of ellipsoidal bubble properties, three-dimensional system with sufficient small grid scales was considered by Olmos et al. (2001) and small spherical bubbles were simulated using a two-dimensional axi-symmetric model. Sankaranarayanan and Sundaresan (2002) simulated a single bubble rising in water in a periodical domain. The effect of drag and lift forces, and virtual mass on gas holdup was also studied but their study was restricted to low values of gas fraction. Chen and Fan (2004) applied the VOF approach to simulate bubble motion for two- and three-phase system and studied the rise behaviour of a single gas bubble in liquids but they ignored the effect of superficial gas velocity on bubble size distribution. In a recent study, Bertola et al. (2004) have investigated the influence of bubble diameter and gas holdup on the hydrodynamics of bubble column reactor using VOF approach. The effect of turbulence models on the buoyant motion of fluids has been studied in detail by Cartland Glover and Generalis (2004) and gas-liquid simulation on an airlift bubble column reactor were successfully compared with the experimental data by Blazej et al. (2004). Although all of these studies are useful in predicting the dynamic behaviour of single bubbles, they cannot accurately model the behaviour of a bubble swarm. Also, because of the bubblebubble interaction in an actual reactor, a model for the interphase forces based on single bubble simulations does not provide good results when adapted to simulate a bubble swarm. Finally, a bubble stream causes an up-welling flow that results in the bubbles to rise faster than if they were alone therefore, the information obtained on the rising of one isolated bubble cannot be considered valid. Despite all of these significant efforts in predicting the characteristics of bubbly flows with single bubbles, the hindering trajectory of bubbles in an otherwise stagnant fluid is still one of least understood aspects of bubble hydrodynamics. Consequently, there is need to develop both qualitative as well as quantitative understanding about this aspect of bubbles by performing multi-bubble simulations.

In this study, we have simulated hydrodynamics of a continuous chain of bubbles for bubble column reactor. The VOF approach available in the commercial CFD code FLUENT has been used in two- and three-dimensional simulations. The effect of gas distributor (hole-size) and superficial gas velocity on the bubble size distribution, bubble rise velocity and its trajectory has been investigated by performing simulations on a 20 cm diameter and 1 m high cylindrical bubble column reactor. The predictions have been compared with the experimental studies reported in the literature.

GOVERNING EQUATIONS AND NUMERICAL SCHEME

FLUENT 6.1 has been used to simulate the motion of a continuous chain of bubbles rising in the stagnant liquid. In the FLUENT's VOF model, the movement of the gas-liquid interface is tracked based on the distribution of [[alpha].sub.G], the volume fraction of gas in a computational cell, where [[alpha].sub.G] = 0 in the liquid phase and [[alpha].sub.G] = 1 in the gas phase. Therefore, the gas-liquid interface exists in the cell where [[alpha].sub.G] lies between 0 and 1. The geometric reconstruction scheme that is based on the piece linear interface calculation (PLIC) method was applied to reconstruct the bubble-free surface. The surface tension was included by assigning it a constant value (0.072 N/m), while the turbulence was modelled with the k-[epsilon]?turbulence model.

Modelling Equations

The VOF approach used to simulate bubble motion is based on the Navier-Stokes equations for a mixture phase and are described as bellow:

The Continuity Equation

[delta]/[delta]t ([rho]) + [nabla].([rho] v) = 0 (1)

The Momentum Equation

A single momentum equation, which is solved throughout the domain and shared by all the phases, is given by:

[delta]/[delta]t ([rho]v) + [nabla].([rho] v v) = -[nabla] p + [nabla].[[micro] ([nabla]v + [[nabla]v.sup.-T]] + [rho]g + F (2A)

where F is the surface tension force and according to Brackbill et al. (1992), it may be expressed as:

F = [sigma]K[nabla][rho]/<[rho]> (2B)

where [sigma] is the surface tension, K is the curvature of the surface, and <[rho]> is the average density:

<[rho]> = [[rho].sub.L] + [[rho].sub.G]/2 (2C)

The Volume Fraction Equation

The tracking of the interface between the gas and liquid is accomplished by the solution of a continuity equation for the volume fraction of gas, which is:

[delta]/[delta]t ([[alpha].sub.G]) + v.[nabla] [[alpha].sub.G] = 0 (3)

The volume fraction equation is not solved for the liquid; the liquid volume fraction is computed based on the following constraint:

[[alpha].sub.G] + [[alpha].sub.L] = 1 (4)

where [[alpha].sub.G] and [[alpha].sub.L] is the volume fraction of gas and liquid phase, respectively.

Surface Tension

The surface tension model in FLUENT is the continuum surface force (CSF) model proposed by Brackbill et al. (1992) (FLUENT 6.1 Manual). In this approach, the surface tension in the VOF calculation is modelled as a source term in the momentum equation at fluid interfaces having finite thickness. In this study, simulations are performed with a constant value of the source term (0.072 N/m).

In the interfacial multiphase flows, the importance of surface tension effects is determined on the basis of the Weber number:

We = [[rho].sub.G][U.sub.b]D/[sigma] (5)

where [[rho].sub.G] is the density of gas, D is the diameter of the bubble, [U.sub.b] is the bubble rise velocity and [sigma] is the surface tension. Surface tension effects can be neglected if We >> 1 and in other cases the Weber number determines the shapes and motion of bubbles (Clift et al., 1978). However, because the motion of individual bubbles was not tracked hence the effect We on the bubble shapes and motion was not evaluated.

Turbulence Modelling

The turbulence in the continuous phase has been modelled using a modified k-[epsilon]?turbulence model available in FLUENT 6.1, which is the most widely used turbulence model to simulate turbulence eddies. This model accounts for the transport not only of the turbulence velocity scale but also of the length scale. It employs a transport equation for the length scale that allows the length scale distribution to be determined even in complex flow situations like those in bubble column reactor.

Differencing Schemes

A first-order up-wind differencing scheme was applied for the solution of momentum equation. The pressure-implicit with splitting of operators (PISO) pressure-velocity-coupling scheme, a member of the SIMPLE family of algorithms, was used for the pressure-velocity-coupling scheme, which is recommended for usual transient calculations. Use of PISO allows for a rapid rate of convergence without any significant loss of accuracy. Pressure was discretized with a PRESTO scheme. Other schemes (linear or second-order schemes) took longer to converge (a comparison is shown in Figure 14). As large body forces (namely, gravity and surface tension forces) exist in multiphase flows, the body force and pressure gradient terms in the momentum equation were almost in equilibrium, with the contributions of convective and viscous terms small in comparison. Segregated algorithms converge poorly unless partial equilibrium of pressure gradient and body forces is taken into account. FLUENT provides an optional "implicit body force" treatment that can account for this effect, making the solution more robust. The volume fraction equation for gas (Equation (3)) was solved using an explicit time-marching scheme and the maximum allowed Courant number was set to 0.25. Under relaxation factor used for pressure and momentum were 0.6 and 0.4, respectively. For turbulence parameters, intensity and hydraulic diameter specifications were used. A time step value of 5 x [10.sup.-4] s was used throughout the simulations. At start, for few time steps, the solution was converged in about 50 iteration per time step. Once the solution was converged then it took less than 30 iterations per time step. Simulations were run for 5 and 2.5 s in two- and three-dimensional simulations, respectively.

[FIGURE 1 OMITTED]

Physical Properties

The properties of air and water were used in the transport equations when the computational cell was in pure liquid or pure gas phase, respectively. At interface between the gas and liquid phases, the mixture properties of gas and liquid phases were based on the volume fraction weighted average. The density and viscosity in each cell at interface were computed by the application of following equations:

[rho] = [[alpha].sub.G] [[rho].sub.G] + (1 - [[alpha].sub.G]) [[rho].sub.L] (6A)

[micro] = [[alpha].sub.G] [[micro].sub.G] + (1 - [[alpha].sub.G]) [[micro].sub.L] (6A)

where [[rho].sub.G], [[rho].sub.L], [[micro].sub.G] and [[micro].sub.L] is density and viscosity of gas and liquid phase, respectively, while [[alpha].sub.G] is the volume fraction of gas.

Interface Tracking

To overcome the problem of numerical diffusion which most standard differencing schemes suffer, the geometric reconstruction was performed with piecewise-linear (PLIC) scheme, which is the most accurate scheme in FLUENT and is also applicable for general unstructured meshes (Fluent 6.1 Manual). It assumes that the interface between two fluids has a linear slope within each cell, and uses this linear shape for calculating the advection of fluid through the cell faces. The first step in this reconstruction scheme is to calculate the position of the linear interface relative to the centre of each partially filled cell, which is based on the information about the volume fraction and its derivatives in the cell. The second step is to calculate the adverting amount of fluid through each face using computed linear interface representations and the information about the normal and tangential velocity distribution on the face. The third step is to calculate the volume fraction in each cell using the balance of fluxes calculated during the previous step (Rider and Kothe, 1998).

Domain Description

Simulations were performed both in two- and three-dimensional manners with different size distributors. Table 1 summarizes different geometrical configurations used in the present work. Results with gas distributor

size ranging from 2-10 cm are reported in this paper. Simulations were run for different mesh sizes and reported numbers of cells are for a grid independent solution. Detail about mesh size used as well as number of grids and solution time taken by different cases has also been presented in Table 1. Velocity inlet and pressure outlet boundary conditions were applied at inlet and outlet, respectively. At the walls, no-slip boundary condition was imposed. The column was simulated as an open system, so the pressure in the gas space above the initial liquid column was equal to the ambient pressure (101.325 kPa).

In order to reduce the computational time for these memory intensive calculations, all simulations were performed on high-speed XeonTM dual processor computers (having four 2.66GHz processors). Simulations were performed till a fully developed flow field was ensured by examining the overall mass balance and time history of the relevant flow variables. Depending on the distributor size, each simulation required between approximately 100 000--250 000 iterations, which corresponds to 1-5 d of CPU time. Overall continuity equation and the volume fractions of the phases were used to monitor the convergence to fully developed flow.

RESULTS AND DISCUSSIONS

Simulations were performed with a VOF approach to study the effect of superficial gas velocity and gas distributor on the bubble size distribution and hydrodynamics (bubble rise velocity and bubble trajectory). The effect of superficial gas velocity on the bubble hydrodynamics and gas holdup was studied as a function of the gas distributor and superficial gas velocity. All simulations were carried out both in two- and three-dimensions. A quantitative comparison of the bubble rise velocity and volume average gas holdup was made with correlations from the literature (Deckwer, 1992; Hikita et al., 1980; Joshi and Sharma, 1979). The qualitative comparison of the bubble trajectory and bubble size distribution was also in good agreement with the experimental observations (Clift et al., 1978).

Effect of Superficial Gas Velocity on Bubble Size Distribution

In order to study the effect of superficial gas velocity on bubble size distribution, simulations were performed with a single-hole distributor (hole size = 10 cm) and different superficial gas velocities ranging from 1 - 10 cm/s. Figures 2 and 3 show snapshots of volume fraction of air for two- and three-dimensional simulations, respectively. It is clear that bubble size was a function of superficial gas velocity, which is in accordance with the experimental observations of Clift et al. (1978) and other researchers (Chen and Fan, 2004). Small bubbles produced at low superficial gas velocity (for example with 1 cm/s) had low bubble rise velocity, which was measured using the rise position of the nose of the initial bubble produced from the orifice (Figure 2). It was also noted that the coalescence times with the air-water interface (i.e time taken by the leading bubble to emerge from liquid) with 1.0, 5.0 and 10 cm/s superficial gas velocity were 4.5, 3.5 and 2.75 s, respectively (precise calculation of the time was made by saving data files with small intervals). The leading bubble produced from orifice had bigger diameter when compared to the trailing bubbles which might be due to less effects of wall and other surface forces on the trailing bubbles. The trailing bubbles were observed to move in a rectilinear manner at low superficial gas velocities (1.0 cm/s) while at higher superficial gas velocities these bubbles exhibit a slightly zigzag or oscillatory behaviour. This can be explained by considering the behaviour of a single rising bubble. A rising bubble pushes the liquid in front of it while the liquid behind it is sucked by the bubble wake which lies directly behind the bubble; therefore, it is due to the greater degree of turbulence and the formation wake behind the bubbles which is responsible for the oscillatory behaviour of the trailing bubbles behind large rising bubbles, which is in accordance with experimental observations of Olmos et al. (2001) and Chen and Fan (2004).

In the case of three-dimensional simulations, the effect of superficial gas velocity similar to that observed in two-dimensional simulations (Figure 3). For example, a similar qualitative increase in bubble size and zigzag trajectory of trailing bubbles was observed with the increase in gas velocity. However, under the same conditions, quantitative comparisons of bubble position have shown some interesting trends in bubble rise velocity in case of 3-D simulations when compared with those shown under 2-D simulations. For example, at 5 cm/s superficial velocity, the first bubble ejected in 2-D simulation had taken 3.5 s to reach to the top of the column (coalescence) while in case of 3-D simulations, it took only 2.5 s. This difference in 2-D and 3-D simulations may attributed to the way 2-D simulation are formulated in FLUENT (as well as in most CFD calculations) as the depth of the geometry in a 2-D simulation is considered as unity (in FLUENT 1 m) and therefore, a circle in 2-D is essentially a cylinder of unit length

In order to validate these qualitative results in Figures 2 and 3, bubble rise velocities were calculated from a linear regression of the z-coordinates of the nose of the bubble during its upward motion (a similar approach to that used by Krishna et al. (2000)). Figure 4 shows the position of z-coordinate of the nose of bubble at different times. These curves were used to determine the rise velocities for both two- and three-dimensional simulations. A fifth-order polynomial was fitted to get more accurate values of bubble rise velocity. Calculated values of bubble rise velocities are plotted against superficial gas velocity and are shown in Figure 5. This quantitative increase in bubble rise velocity with increase in superficial gas velocity is consistent with the qualitative behaviour. The bubble rise velocities were varied from 34-45 cm/s in 3-D and from 24.5-27 cm/s in 2-D simulations. Calculated rise velocities in 2-D simulations were ~ 30% lower than those in 3-D simulations which is in accordance with Krishna et al. (2000) results obtained with 2-D circular cap bubble and 3-D spherical cap bubble. This discrepancy with 2-D simulations can be explained on the basis of the 3-D wake exhibited by the bubbles (Olmos et al., 2001), which cannot be accurately modelled with 2-D simulations. Therefore, for a precise prediction of bubble rise velocity and other hydrodynamic characteristics of bubble, 3-D simulation is better choice. In order to provide further verification of CFD simulations, a comparison of bubble rise velocity was made with experimental correlation from Deckwer (1992). Dotted line plotted with Deckwer equation in Figure 5 was very close to the calculated rise velocities with 3-D simulations--particularly at lower superficial gas velocities. A slight deviation at higher superficial gas velocities could be attributed to the use of single 10 cm distributor instead of actual gas distributors. The results with 2-D simulations were significantly different from Deckwer's correlation which indicates the significance of 3-D simulations for bubble rise velocity prediction.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

Effect of Gas Distributor (Hole size) on Bubble Size Distribution

Effect of distributor hole size on bubble size distribution and hydrodynamics was studied at constant superficial gas velocity (1.0 cm/s). Three distributors with hole size ranging from 2-10 cm (Table 1) were selected for two- and three-dimensional simulations. Figure 6 shows two-dimensional effect of gas distributor on the bubble size distribution. Although these results appear similar if just looking at bubble sizes, but one can observe decrease in rise velocity with a decrease in the hole-size which for example, becomes quite apparent by comparing bubble rise trajectories for 5 and 2 cm size distributor (with the average plume rise velocity of 20.7 cm/s for 5 cm hole and 16.7 cm/s for 2 cm hole). However, there was relatively smaller difference in the rise velocities between the simulations for 5 and 10 cm holes (20.7 cm/s and 21.2 cm/s, respectively). This may be attributed to the predominance of the wall effect for the relatively bigger leading bubbles coming from these larger holes. Therefore, it can be inferred that the bubble size has similar dependency on the hole-size as superficial gas velocity. Small bubbles usually formed from small size distributor and have low rise velocities. A steady-state rise of these bubbles without much deflection from centre position can also be observed. Three-dimensional simulations shown in Figure 7 show a similar behaviour. Relatively bigger size bubbles were formed with 10 cm hole when compared with 5 and 2 cm and their rise time to the same height of column was less. A comparison between Figures 6 and 7 shows that at relatively low superficial gas velocity, bubbles exhibited rectilinear trajectory. These results are in good agreement with the earlier work reported by Cartland Glover and Generalis (2004) for an identical problem by using an algebraic slip mixture model. In that study, inlet velocity conditions were applied to 80% of the base of the column in contrast to the 50% in present work, still the rectilinear trajectory of bubbles is comparable for the range of superficial gas velocity considered in our simulations.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

Figure 8 plots calculated bubble rise velocities as a function of hole-size at constant superficial gas velocity (1.0 cm/s). For the investigated range of hole-sizes (2-10 cm), bubble rise velocities varied from 13-23 cm/s and 15-34 cm/s in two- and three-dimensional simulations, respectively. Interestingly, as in bubble rise velocity curves (Figure 5), there was an offset between 2-D and 3-D results (10 - 30%), however, this difference was not prominent with small size distributor at low superficial gas velocity. When compared with Deckwer's equation, both 2-D and 3-D simulations failed to agree at smaller distributor diameters. However, this apparent disagreement can be attributed to the use of single hole-size in these simulations that resulted in very low perforations compared to those used in Deckwer's work. Even so, 2-D simulations, in general gave a very poor comparison with Deckwer's equation throughout the entire domain. In contrast, 3-D simulations gave reasonable quantitative approximation of bubble rise velocities, particularly those with bigger distributor hole-sizes.

[FIGURE 7 OMITTED]

Effect on Bubble Shape and Rise Trajectory

Qualitative effects of superficial gas velocity and gas distributor on bubble shape and bubble rise trajectories are shown in Figure 9. Effect of gas velocity was studied for the distributor 5 cm diameter at two different velocities (0.25 and 1.0 cm/s). Comparison between Figures 9 (a) and (b) shows a formation of small spherical bubbles at low velocities and large ellipsoidal bubbles at higher velocities. The obtained bubble shapes are in agreement with previous experimental observations (Clift et al., 1978). In Figure 9 (a) an oscillatory or zigzagging behaviour was exhibited by small spherical bubbles which are produced at low superficial gas velocities. It can also be seen that the onset of path oscillation, which was dominated by the vortex behaviour, was much slower than the natural shape relaxation. The observed behaviour is consistent with experimental results reported by Liu et al. (2005b). This zigzagging behaviour could be attributed due to the ability of small bubbles to travel toward the wall as indicated by Tomiyama et al. (1993). It is also evident from Figure 9 (a) that nearly same size bubbles were produced at comparatively low superficial gas velocities.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

In Figure 9 (b) ellipsoidal cap bubbles can be observed at high superficial gas velocity as compared to Figure 9 (a). The bigger size bubbles exhibited rectilinear trajectory, which is also in accordance with Liu et al. (2005a) experimental observations.

Effect of distributor hole-size on the bubble shape and rise trajectory is shown in Figure 9 (c). For this simulation, gas was injected through a 20 cm distributor with a superficial velocity of 10 cm/s. Consequently, a very big ellipsoidal bubble was produced which moved in a plug flow manner. As discussed in the previous section, the large leading bubble has a rectilinear movement and a wake is expected to be formed behind it. Due to the wake formation, the trailing bubbles exhibit an oscillatory behaviour. One can also observe a shedding of small bubbles from large bubble, which is a limitation of the VOF approach and causes inaccurate results for interfaces in a high shear flow as also mentioned by Delnoij et al. (1999). Variation in the shape and size of trailing bubbles from the leading bubble is also obvious in Figure 9 (c). It can be explained on the basis of the fact that the trailing bubbles are rising in the wake region of the leading bubble; hence one can expect the effect of high degree of turbulence on the shape and size of the trailing bubbles. This is the main reason for the formation of bubbles with different size and shape with a wobbling motion in that region.

[FIGURE 10 OMITTED]

Effect of Grid Size and Discretization

The effect of mesh size on the accuracy of numerical solutions was investigated for the 2-D column by performing simulations with three different grids (size ranging from 3 to 6 mm) for one of the gas distributors (10 cm diameter). Comparative results for the calculated shape of the bubbles at t = 2.5 s are depicted in Figures 10 and 11 where contours of volume fraction of gas have been overlayed on meshed planes. Comparison of these snapshots shows only minor effect of the grid size on the qualitative behaviour of bubbles. To demonstrate the effect of mesh size quantitatively, its influence on the calculated bubble rise velocities was also studied. It is clear from Figure 12, that the difference between bubble rise velocity patterns in three plots was insignificant. Therefore, it may be concluded that the intermediate grid size (5 mm) yielded solutions with sufficient accuracy (with respect to the trade-off between computational cost and the features being studied), and was chosen for the results presented in this paper.

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

In order to study the effect of discretization scheme on the numerical solution, simulations were performed with both first-order and third-order Quick schemes available in FLUENT. Figure 13 shows a qualitative comparison between the discretization schemes and it is clear that the two plumes are very similar. Therefore, on a qualitative basis use of the third-order scheme did not provide any significant advantage although it took about 50% more time to converge than the first-order scheme. Furthermore, use of the PRESTO scheme for pressure discretization (in Figure 14) gave identical results to those with the first- or third-order schemes. However, with the body force weighted scheme for pressure discretization the bubble rise velocity was excessively low.

[FIGURE 13 OMITTED]

MODEL VALIDATION

In order to achieve quantitative validation of simulation results, volume-averaged gas holdups for different distributors were plotted as a function of superficial gas velocity. As shown in Figure 15, the gas holdup increased with superficial gas velocity and this was in accordance with the data available in the literature. Predicted values of volume-averaged gas hold ups for 3-D simulations have been compared with the experimental work of Joshi and Sharma (1979), and Hikita et al. (1980). It is clear that simulation results showed a near-complete agreement with the experimental work of Joshi and Sharma (1979) throughout the superficial velocity range studied. However, the correlation of Hikita et al. (1980) showed significant deviations particularly those at higher superficial gas velocities. It may be observed from these plots that there was slight increase in gas holdup value with a decrease in hole-size, which can be attributed due to the formation of smaller bubbles from smaller holes. Interestingly, this dependency of gas holdup on superficial gas velocity was less at relatively low and high superficial gas velocities, which is also in agreement with the Andou et al. (1996) experiments.

[FIGURE 14 OMITTED]

An experimental validation of the present VOF model was achieved by performing a series of experiments with 10 cm (ID) and 80 cm high column at Chemical Reaction Engineering Laboratory (CREL), Washington University in St. Louis, USA. A single hole distributor (5 mm) was utilized for this work. Figure 16 compared overall gas holdup results obtained from experiments and simulation. A close agreement between experiments and 3-D simulations is evident which has further validated the significance of VOF simulations with for bubble column reactors.

CONCLUSIONS

In this paper, a volume of fluid (VOF) approach was used to study the effect of air distributor on the hydrodynamics of a bubble column reactor. The effect of superficial gas velocity and distributor hole-size on the bubble hydrodynamics were studied both in two- and three-dimensional geometries. Both bubble rise velocity and bubble size increased with increase in gas superficial velocity. The bubble rise velocity, bubble shape, typical rise trajectories and gas holdup were compared both qualitatively and quantitatively with the experimental data available in the literature. In general, 3-D VOF model computed the shape and the motion of continuous bubble chain more accurately than those using 2-D simulations.

[FIGURE 15 OMITTED]

[FIGURE 16 OMITTED]

ACKNOWLEDGEMENT

The authors gratefully acknowledge the financial support provided by the Australian Research Council (ARC) to conduct this research (Grant Application ID: DP0451314). MAA also thanks the Western Australian Energy Research Alliance (WAERA) for providing a top-up scholarship. Authors also express their thanks to Prof Dudukovic and Prof Al-Dahan for granting permission to MAA to perform experiments in their labs at Washington University at St Louis.

NOMENCLATURE

Ug gas superficial velocity [[ms.sup.-1]]

v velocity vector [[ms.sup.-1]]

g acceleration due to gravity [[ms.sup.-2]]

F external body force [N]

t time [s]

p static pressure [Pa]

D column diameter [m]

l characteristic length [m]

d hole diameter [m]

Greek Symbols

[[alpha].sub.G] volume fraction of the gas phase in the computational cel[--]

[[alpha].sub.L] volume fraction of the liquid phase in the computational cell [--]

[[rho].sub.L] liquid density [kg[m.sup.-3]]

[[alpha].sub.G] gas density [kg[m.sup.-3]]

[[micro].sub.L] liquid viscosity [kg[m.sup.-1][s.sup.-1]]

[[micro].sub.G] gas viscosity [kg[m.sup.-1][s.sup.-1]]

[sigma] surface tension [N/m]

Manuscript received August 10, 2006; revised manuscript received January 11, 2007; accepted for publication January 15, 2007.

REFERENCES

Akhtar, M. A., V. K. Pareek and M. O. Tade, "Two-Fluid Eulerian Simulation of Bubble Column Reactors with Distributors," Journal of Chemical Engineering of Japan submitted (2004).

Andou, S., K. Yamagiwa and A. Ohkawa, "Effect of Gas Sparger Type on Operational Characteristics of a Bubble Column under Mechanical Foam Control," Journal of Chemical Technology and Biotechnology 66, 65-71 (1996).

Blazej, M., G. M. Cartland Glover, S. C., Generalis and J. Markos, "Gas-Liquid Simulation of an Airlift Bubble Column Reactor," Chemical Engineering and Processing 43, 137-144 (2004).

Bertola, F., G. Baldi, D. Marchisio and M. Vanni, "Momentum Transfer in a Swarm of Bubbles: Estimates from Fluid-Dynamic Simulations," Chemical Engineering Science 59, 5209-5215 (2004).

Brackbill, J. U., D. B. Kothe and C. Zemach, "A Continuum Method for Modelling Surface Tension," Journal of Computational Physics, 100, 335-354 (1992).

Cartland Glover, G. M. and S. C. Generalis, "The Modelling of Buoyancy Driven Flow in Bubble Columns," Chemical Engineering and Processing 43, 101-115 (2004).

Chen, C. and L.-S. Fan, "Discrete Simulation of Gas-Liquid Bubble Columns and Gas-Liquid-Solid Fluidized Beds," AIChE Journal 50, 288-300 (2004).

Chen, P., J. Sanyal and M. P. Dudukovic, "Numerical Simulation of Bubble Columns Flows: Effect of Different Breakup and Coalescence Closures," Chemical Engineering Science 60, 1085-1101 (2005).

Clift, R., J. R. Grace and M. E. Weber, "Bubbles, Drops, and Particles," Academic Press, New York (1978).

Deckwer, W.-D., "Bubble Column Reactors," John Wiley and Sons, New York (1992).

Degaleesan, S., M. Dudukovic and Y. Pan, "Experimental Study of Gas-Induced Liquid-Flow Structures in Bubble Columns," AIChE Journal 47, 1913-1931 (2001).

Delnoij, E., J. A. M. Kuipers and W. P. M. van Swaaij, "Computational Fluid Dynamics Applied to Gas-Liquid Contactors," Chemical Engineering Science 52, 3623-3638 (1997).

Delnoij, E., J. A. M. Kuipers and W. P. M. van Swaaij, "A Three-Dimensional CFD Model for Gas-Liquid Bubble Columns," Chemical Engineering Science 54, 2217-2226 (1999).

Essemiani, K., G. Ducom, C. Cabassud and A. Line, "Spherical Cap Bubbles in a Flat Sheet Nanofiltration Module: Experiments and Numerical Simulation," Chemical Engineering Science 56, 6321-6327 (2001).

Hikita, H., S. Asai, K. Tanigawa, K. Segawa and K. Kitao, "Holdup in Bubble Columns," Chemical Engineering Journal 20, 59 (1980).

Hirt, C. W. and B. D. Nichols, "Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries," Journal of Computational Physics 39, 201-225 (1981).

Joshi, J. B. and M. Sharma, "A Circulation Cell Model for Bubble Columns," Trans Inst. Chem. Eng. 57, 244-251 (1979).

Karamanev, D. G., "Rise of Gas Bubbles in Quiescent Liquids," AIChE Journal 40, 1418-1421 (1994).

Krishna, R. and J. M. van Baten, "Rise Characteristics of Gas Bubbles in a 2D Rectangular Column: VOF Simulations vs Experiments," International Communications in Heat and Mass Transfer 26, 965-974 (1999).

Krishna, R. and J. M. van Baten, "CFD Simulations of a Bubble Column Operating in the Homogeneous and Heterogeneous Flow Regimes," Chemical Engineering Technology 25, 1081-1086 (2002).

Lain, S., D. Broder and M. Sommerfeld, "Experimental and Numerical Studies of the Hydrodynamics in a Bubble Column," Chemical Engineering Science 54, 4913-4920 (1999).

Lehr, F., M. Millies and D. Mewes, "Bubble-Size Distribution and Flow Fields in Bubble Columns," AIChE Journal 48, 2426-2443 (2002).

Li, Y., J. Zhang and L.-S. Fan, "Discrete-Phase Simulation of Single Bubble Rise Behavior at Elevated Pressures in a Bubble Column," Chemical Engineering Science 55, 4597-4609 (2000).

Liu, Z., Y. Zheng, L. Jia and Q. Zhang, "Study of Bubble Induced Flow Structure using PIV," Chemical Engineering Science 60, 3537-3552 (2005a).

Liu, Z., Y. Zheng and Q. Zhang, "PIV Study of Bubble Rising Behavior," in "7th World Congress of Chemical Engineering," IChemE, Glasgow (2005b).

Olmos, E., C. Gentric, C. Vial, G. Wild and N. Midoux, "Numerical Simulation of Multiphase Fflow in Bubble Column Reactors. Influence of Bubble Coalescence and Break-up," Chemical Engineering Science 56, 6359-6365 (2001).

Ranade, V. V., "Computational Flow Modeling for Chemical Reactor Engineering," Academic Press, New York (2002).

Rider, W. J. and D. B. Kothe, "Reconstructing Volume Tracking," Journal of Computational Physics 141, 112-152 (1998).

Sankaranarayanan, K. and S. Sundaresan, "Lift Force in Bubble Suspensions," Chemical Engineering Science 57, 3521-3542 (2002).

Sokolichin, A., G. Eigenberger and A. Laptin, "Simulation of Buoyancy Driven Bubbly Flow: Established Simplifications and Open Questions," AIChE Journal 50, 24-45 (2004).

Tomiyama, A., I. Zun, A. Sou and T. Sckaguchi, "Numerical Analysis of Bubble Motion with the VOF Method," Nuclear Engineering and Design 141, 69-82 (1993).

Van Sint Annaland, M., N. G. Deen and J. A. M. Kuipers, "Numerical Simulation of Gas Bubbles Behaviour using Three-Dimensional Volume of Fluid Method," Chemical Engineering Science 60, 2999-3011 (2005a).

Van Sint Annaland, M., N. G. Deen and J. A. M. Kuipers, "Numerical Simulation of Gas-Liquid-Solid Flows using a Combined Front Tracking and Discrete Particle Method," Chemical Engineering Science 60, 6188-6198 (2005b).

Van Wachem, B. G. M. and A. E. Almstedt, "Methods for Multiphase Computational Fluid Dynamics," Chemical Engineering Journal 96, 81-98 (2003).

Zhang, X. and G. Ahmadi, "Eulerian-Lagrangian Simulations of Liquid-Gas-Solid Flows in Three-Phase Slurry Reactors," Chemical Engineering Science 60, 5089-5104 (2005).

M. Abid Akhtar, Moses Tade and Vishnu Pareek *

Process Systems Computations Laboratory, Department of Chemical Engineering, Curtin University of Technology, GPO Box U1987, Perth, Western Australia

* Author to whom correspondence may be addressed.

E-mail address: [email protected]

Table 1. Various configurations of distributors used for simulation
2-D 3-D

Hole size Ug (cm/s) Mesh size No of Time (h)
D (cm) (mm) mesh nodes
10 1,5,10 0.46 9 900 8
5 0.25,1,5,10 0.44 10 020 9
2 0.01,1,5 0.4 12 050 10

Hole size Hole size Ug (cm/s) Mesh size No of Time (h)
D (cm) D (cm) (mm) mesh nodes

10 10 1,5,10 0.5 33 440 25
5 5 1,510 0.25 90 984 72
2 2 1,510 0.2 97 512 78