Abstract: Computational Fluid Dynamics (CFD) simulation is highly helpful for the study of large bubble column hydrodynamics because a validated model of a small reactor can be easily scaled up. Until recently CFD had not been used to study nozzle diameter effect on the hydrodynamics in bubble column. So, an attempt is made to simulate the flow pattern numerically so that scaling of numerical model can be done easily without depending on experiments. This research demonstrates the numerical simulation for bubble column using standard k-ε turbulence model. The problem is solved using Algebraic Slip Mixture Model (ASMM). The ASMM model solves conservation equations for only mixture rather than considering equations for both the phases individually. The simulation allows us to determine gas holdup, a key variable that need to be maximized to have a high reaction rates. A computer simulation of gas holdup is a complicated problem because the flow is generally turbulent and it requires tracking the transient movement of each bubble. The simulation predicts the flow pattern with different nozzle diameters and the use of sparger. The flow pattern was observed in simulations by considering physics of the problem.

INTRODUCTION

In chemical, biochemical industries for different process like absorption, fermentations, coal liquefaction and waste water treatment, etc., bubble columns are important as a simple contactor of two and three phase reactor. It offers many advantages over other conventional multiphase reactor due to their simple construction, low operating cost, high energy efficiency and good mass and heat transfer. In bubble column operation with low gas flowrate demonstrates the homogeneous bubbly flow while at relatively higher gas flow rates shows the heterogeneous flow regime with non-uniform bubble size. In homogeneous flow bubble size distribution and gas holdup are most desirable design parameters for efficient mass transfer.

Other hydrodynamic characteristics like distributor diameter, energy distribution and the flow characteristics in the column are also important in its design which are mostly depends on the column geometry, operating conditions, physical properties of the two phases and type of gas distributor. Since, there are different variables affect the design parameter in bubble column and due to the complexity of their combined effect on the design parameter various empirical or semi empirical correlations are developed to predict the design parameter and to scale up the reactor. The scale up of bubble columns is still poorly understood due to the complexity of flow patterns and their unknown behavior under different sets of design parameters. To establish the relationship between adjustable design and operating parameters correlation model and experimental observations are key points in multiphase reactor (Ranade, 2002). Other then correlation model various researchers studied the hydrodynamics in the bubble column from the theoretical model and their simulation of the bubble column gave insight of the column to interpret the various hydrodynamic characteristics in the column. During the last 4 decades, effort on understanding the complex flow behavior in bubble column is going on. Rafique et al. (2004) and Jakobsen et al. (1997) reported the extensive reviews on modeling and simulation of bubble columns. The complex flow pattern in bubble columns has inhibited the development of design of ideal bubble column. Therefore, understanding the relationship between the flow pattern and the design parameters such as pressure drop, fractional gas hold-up, energy distribution, flow behavior, etc., are still required in this area. Sokolichin and Eigenberger (1999) reported that a 2-dimensional representation of the turbulent flow in such columns k-ε turbulence model is stationary due to the stabilizing action of the turbulent eddy viscosity. Eddy viscosity plays an important role on the hydrodynamics like interfacial drag and the shear stress phenomena. Based on shear-induced and one bubble-induced eddy viscosity model, the eddy viscosity is one order of magnitude smaller than the contribution from the shear-induced turbulence (Pan et al., 2000). Lain et al. (1999) and Annaland et al. (2005) have successfully installed eddy simulation turbulence model in combination with the Euler-Lagrange model and obtained good agreement with experimental data in a square and cylindrical bubble column. Pfleger et al. (1999) and Sokolichin and Eigenberger (1999) enunciated the 2-dimensions and 3-dimensions simulations based on the principle of differences in effective viscosity. Buwa and Ranade (2002) and Oey et al. (2003) reported the same mechanism of the effect of viscosity on the simulation of bubble column.

Other various researchers reported their extensive research on studies on different hydrodynamic characteristics by CFD with different operating conditions (Buwa and Ranade, 2002; Jakobsen et al., 1997; Deen et al., 2001; Bove et al., 2004). In the present researchers, it is attempted to simulate the flow of air-water in a bubble column of height 2 and 0.05 m diameter for upflow problems to enunciate briefly the effect of gas distributor diameter on the gas holdup, axial liquid velocity and the energy distribution in bubble column.

MATERIALS AND METHODS

About 2 different theoretical models are generally used to simulate the hydrodynamics in bubble column reactor: One is Eulerian multiphase model and the other is Algebraic Slip Mixture Model (ASMM). In this study the Algebraic Slip Mixture Model (ASMM) were used for simulation due to homogeneity behavior of bubble column and incorporated in Fluent CFD software.

The algebraic slip mixture model: The details of mixture model have been reported by different researchers (Manninen et al., 1996; Blazej et al., 2004; Chen et al., 2005). The standard k-ε model is used to model the mixture turbulence. This model describes the flow regime as a single-phase pseudo-continuous mixture of the gas and liquid phases. This means that a single set of continuity and momentum equations can be used to model the flow phenomena. The principal assumptions in this formulation are as: the phases as 2 interpenetrating continua with the probability of existence of each phase at a point in the computational domain given by its respective phase holdup. A single equation is solved for continuity of the mixture and a single equation is solved for the momentum of the mixture. In any control volume motion of each phase relative to the center of mass of the mixture is viewed as a diffusion of that phase. For the mixture, the turbulent stress term in the mixture equation is closed by solving a k-ε model. The basic conservations equations for ASMM are formulated based on the above assumptions as follows: Continuity equation for the mixture:

(1)

Momentum equation for the mixture:

(2)

where, ρ_mix is mixture density given by:

(3)

and mixture viscosity:

(4)

The momentum equations including the interactions between each phase as a drift or slip velocity effect. Drift velocity:

(5)

Slip velocity equation:

(6)

Mass averaged velocity:

(7)

Friction factor: As per Blazej et al. (2004) the friction factor is calculated by Eq. 8:

(8)

where, Reynolds number is defined by:

(9)

The slip velocity effect depends on the volume fraction of the dispersed and continuous phases which is control by a volume fraction equation (Blazej et al., 2004) as given by Eq. 10.

Volume fraction equation:

(10)

Turbulence modeling: A fluctuating velocity field characterizes the turbulent flows in the bubble column. These fluctuations no matter how small mix the momentum, energy and species concentration. These small scale fluctuations with high frequency are computationally expensive to simulate the engineering problems. To remove the small scales, the exact governing equations can be time-averaged, ensemble-averaged or otherwise manipulated to result in a modified set of equations to easily solve less expensively. In this regard, the modified equations contain additional unknown variables which are to be determined in terms of known quantities. Here, we have used standard k-ε turbulence model because of its widespread applicability and simplicity in computing complex problems.

Standard k-ε turbulence model: The standard k-ε model is a semi-empirical model based on model transport equations for the turbulence kinetic energy (k) and its dissipation rate (ε). In the derivation of the k-ε model, it is assumed that the flow is fully turbulent and the effects of molecular viscosity are negligible. The standard k-ε model is therefore, valid only for fully turbulent flows. The turbulence kinetic energy, k and its rate of dissipation, ε are obtained from the following transport equations:

(11)

(12)

The generation of turbulence kinetic energy due to the mean velocity gradients and is given by:

(13)

The generation of turbulence kinetic energy due to buoyancy, Gb is represented as:

(14)

Where:

For the standard k-ε models, the default value of Pr_t is 0.85. The coefficient of thermal expansion, β is defined as:

(15)

The turbulent (or eddy) viscosity, μ_t is computed by combining k and ε as follows:

(16)

The buoyancy effects on e are neglected simply by setting G_b to zero in the transport equation for ε. The degree to which ε is affected by the buoyancy is determined by the constant C_3ε and is given by:

(17)

In the present case, it is taken as 1.0 for fluid flow along the axis. The effect of convective heat and mass transfer is assumed negligible in the present research. In the present problem for air-water system the parameters are taken as follows:

(18)

These values are recommended by Launder and Spalding (1974) for a dispersed multiphase system.

Problem setup
Dimension of the BCR considered: In the present researchers, the dimensions of the bubble column considered here is shown in Table 1. The air flowrate used for the computation is 44 lph. The sparger is considered as perforated type of plate which includes 14 inlet nozzles evenly distributed on both sides of the axis of diameter of 2 mm each. The symmetry of the domain about the axis of the column is utilized and computational geometry is thus reduced to only half of the original. The grid generated for discretization was made in GAMBIT software.

Table 1:	The dimensions of the bubble column used

The meshing scheme used here was Map type and elements were of Quad type. Total numbers of elements created are 50x800 (25 in radial direction and 800 in axial direction). As the flow is symmetrical in radial direction about the central axis, the column was divided into half and symmetry condition was given at the right wall so as to reduce the computational effort. The spacing between the grid points was kept 0.002 m in radial and 0.0025 m in axial direction. The meshing applied is uniform and structured in order to avoid any error due to non-uniformity and unstructured nature of the mesh. Also, this very fine mesh size eliminated the need to refine the mesh elements near boundaries or near sharp gradients.

Simulation procedure: The governing equations are solved sequentially Using segregated solution method (i.e., segregated from one another). A check for convergence of the equation set is made. Using discretization method, the discrete value of the scalars in each cell is stored at the cell centers. However, face values are required for the convection terms. So, the values must be interpolated from the cell center values. This is accomplished using an upwind scheme. In the present research, the first order upwind scheme has been used. In first order upwind scheme, quantities at cell faces are determined by assuming that the cell-center values of any field variable represent a cell-average value and hold throughout the entire cell. The face quantities are identical to the cell quantities. The face value is set equal to the cell-center value in the upstream cell. In first order upwind scheme, the face value φ_f is computed using:

(19)

where, φ and ∇φ are the cell-centered value and its gradient in the upstream cell and Δs is the displacement vector from the upstream cell centroid to the face centroid. Here, the gradient is computed using the divergence theorem as:

(20)

Because of the nonlinearity of the equations used here, it is necessary to control the change of φ. In a simple form, the new value of the variable φ within a cell depends upon the old value, φ_old, the computed change in φ, Δφ and the under-relaxation factor, α as follows:

Table 2:	Superficial gas velocities for various nozzle diameters used

(21)

Boundary conditions: Symmetrical conditions are applied at the center axis for all the variables and no-slip condition is applied on all the column walls. The wall roughness constant was set to 1/2.

Inlet and outlet boundary conditions: The inlet was divided into two parts the nozzle part and the wall part. The nozzle part was having the air coming at 44.0 lph into the column. For all the diameters constant air flow rate was maintained to study the hydrodynamics comparatively. Thus, inlet jet velocity, V_n = Q/A_n. For inlet through nozzle, superficial gas velocities considered are shown in Table 2. Since only air was inlet through the nozzle, the volume fraction of water was set to 0. The outlet was set to pressure outlet condition with gauge pressure set to 0, backflow turbulent kinetic energy to 1 m² sec^-2 and also backflow turbulent dissipation rate to 1 m² sec^-3.

Initial conditions: Initially the column was filled fully with water and thus the dull fluid grid points inside the column were initialized volume fraction of water as 1. The velocities in both axial and radial direction for water are initialized to zero and for air as 0 and 5 m sec^-1, respectively.

Parameters used for simulation: Lower values of under-relaxation factors can be used as shown in Table 3 for better stability in simulation but it would take a lot time for the simulations to go on. Similarly, it is obvious that we did not take higher order schemes for most of the variables in order to reduce the computation time (Table 4). In the present simulation, the residual for all of the variables is set at 10^-3 and hence, convergence is obtained with significant accuracy.

RESULTS AND DISCUSSIONS

From the simulation, it is found that the increase in nozzle diameter decreases the number of iterations it needs to attain convergence. This shows that as the nozzle diameter increases, it takes less time to get into a stable state but there is a possibility of weeping so, it is required to optimize in such a way that there is enough gas flow rate for weeping not to occur.

Table 3:	Relaxation factor for the variables

Table 4:	Discretization schemes for the variables

Effect of distributor size on gas holdup: The gas distributor plays an important role to affect the bubble characteristics which in turn, affect the gas holdup, mass transfer and other parameters. In the porous type gas distributors, the relatively small pores cause high pressure drops across the pores. In the gas distributor, small orifices enable the formation of smaller bubbles. Bouaifi et al. (2001) stated that the gas holdup values were higher with a small orifice gas distributor. Lin et al. (2004) studied the influence of different gas distributors. They observed that bubble sizes are much smaller and better distribution in case of porous sinter plate than in case of perforated. Their results also show that the radial profile of the gas holdup in case of sintered plate is much more uniform than that in perforated plate.

Formation of homogeneous regime and the gas holdup values are significantly influenced by the gas distributor. Different nozzle diameters show a decrease in the highest gas holdup point as the diameter increases. The variations of gas hold at different location of the bubble column for different nozzle size of the gas distributor are shown in Fig. 1-3. By using sparger one can see uniform distribution of gas hold-up across the column. Non-uniform distribution of gas bubble effects the gas holdup distribution in the bubble column.

The bubble size distribution becomes wider and the Sauter mean bubble diameter becomes larger as the superficial gas velocity increase (Chen et al., 2005). The superficial gas velocity when the flow is deep into the churn-turbulent regime has little effect on gas holdup of small bubbles. This may be due to the effect of the liquid circulation in the column.

The circulation effect is due to the pressure difference at radial direction due to radial velocity distribution. In the low gas velocity region, the size of bubbles generated from the single nozzle and sparger is larger when generated homogeneous small bubbles. In the high gas velocity region, small bubbles in the bulk liquid are formed by bubble breakup.

The bubble breakup increases the bubble population and the turbulence of the phases. In the high velocity region, bubbles are getting breaked up due to high kinetic energy of the liquid phase supplied by high different nozzles. But as per Andou et al., (1996), the size and distribution of bubbles in the bulk liquid became almost independent of the sparger types. Lower the nozzle diameter greater the kinetic energy distribution which enhance the breakup of bubbles and formation of finer bubbles. But the finer bubbles without coalesesence have a tendency to redistribute due to the circulation of the liquid in the vicinity of the bottom of the column. Therefore, different gas holdup distribution resulted in the vicinity just above the bottom of the bubble column as shown in Fig. 1.

Li et al. (2009) reported from their simulation results that the configuration of gas distributor has an important impact on liquid velocity and local gas hold-up in the vicinity of the gas distributor.

They observed that an increase in the number of gas sparging pipes used in gas distributors is beneficial in improving the gas hold-up but is disadvantageous in reducing bubble size due to a decrease in turbulent kinetic dissipation. Lau et al. (2009) studied the effects of gas distributor on hydrodynamics like overall gas holdup, bubble size distributions and bubble rise velocity distributions in an air-water shallow bubble column. They reported that single nozzle is not suitable for shallow bed operation. While in the absence of solids perforated plate and porous plate distributors have comparable behavior in their operation. The perforated plate distributor enhances the bubble coalescence while porous plate inhibits the bubble coalescence. This may results in different behavior to encounter the flow field in the bubble bed for the two distributors. Also the physical properties of the liquid in the column play an important role for the bubble size distribution based on the kinetic energy distribution inside the column by the liquid jet for different size of nozzle and sparger hole. Dhotre and Joshi (2007) found that the gas chamber configuration has an effect on the uniformity of gas distribution particularly in the sparger region of bubble column reactors. The homogeneous regime results in relatively uniform radial profile of the gas holdup wheras heterogeneous regime results in more parabolic (Hills, 1974; Shaikh and Al-Dahhan, 2005).

Fig. 1:	Gas hold-up vs radial vector plots of different nozzle diameter just above bottom

Fig. 2:	Gas hold-up vs radial vector plots of different nozzle diameter at 0.8 m above the bottom

This may be due to lateral forces like transverse lift force, turbulent dispersion force and wall lubrication force acting on the bubble surface which depends on the bubble size and flow regime. This may cause the variation of gas holdup with radial position.

Effect of distributor size on axial liquid velocity: Recently, Krishna et al. (1999) proposed a Computational Fluid Dynamics (CFD) based model for prediction of holdup profile and axial velocity profile. The axial liquid velocity profile depends on the distribution of gas in the column. The circulation of the gas inside the column is the main factor to variation of the liquid velocity profile in the column.

Fig. 3:	Gas hold-up vs radial vector plots of different nozzle diameter at 2 m above the bottom, i.e., at outlet

The intensity of the circulation along the column axis varies due to the variation of energy distribution. The energy distribution of course is a function of the phase physical properties. Along the core region of the column, the velocity is higher than the other radial position due to the dilution of interfacial stress between wall of the column and the liquid. Also the nozzle diameter has its effect to distribute the intensity of energy. Lower size of the nozzle will distribute the more energy at a particular liquid and gas velocity. At the vicinity of the nozzle or gas distributor, the energy distribution depends on the intensity of the circulation of the fluid. As per (Fig. 4, 5) the maximum axial liquid velocity decreases with the increase of nozzle diameter. The variations of the gas holdup at different locations for a particular gas distributor as shown in Fig. 1-3 may influence the pressure distribution along the column.

At the centre of the column the axial liquid velocity varies compared to other radial position of the column. It is zero at the wall of the bubble column. Due to this variation of the liquid velocity, the fluid gets circulates at different intensity along the column. Also due to the variation of the circulation velocity the momentum exchange results in distribution in kinetic energy for which the flow pattern and the gas holdup will change radically. The sparger axial liquid velocity suddenly changes its direction at some height and rapidly increases. The velocity distribution from sparger is shown in Fig. 6. The radial profiles of the gas holdup and liquid velocity are closely related to each other. As the superficial gas velocity increases the center-line liquid velocity increases. Also the variation of the liquid velocity with respect to the superficial gas velocity depended on both the flow regime and column diameter and physical properties (Degaleesan et al., 2001; Krishna and van Baten, 1999). Forret et al. (2003) stated that the column diameter significantly results in variation on the radial profile of the liquid velocity. The centerline liquid velocity increases from 0.44-1.04 m sec^-1 in column of diameter 0.15-1.0 m. This is due to liquid circulation effect with column diameter.

Effect of nozzle diameter on turbulent kinetic energy: The turbulent kinetic energy peak height decreases with increasing diameter as shown in the Fig. 7, 8 but the intensity increases with height. The turbulent kinetic energy explanation given in the first 3 cases with different nozzle diameters is based on k-ε model. The liquid velocity affects the global recirculation patterns for the gas and liquid phases. The gas phase recirculation size is much narrower being restricted to a thin layer close to the wall. The radial lift force acts to drive the bubbles against the liquid phase vertical velocity gradients which resulted in different kinetic energy distribution based on different distributor size. The kinetic energy variation depends on the gradients and the fluctuation in the liquid velocity at the vicinity of the distributor. For bubbly flow regime due to suppressed turbulence the effect is not very significant (Chen et al., 2005). Reported that the liquid turbulent kinetic energy is practically flat as a function of radius.

Fig. 4:	Axial liquid velocity profiles for columns of nozzle diameters (a) 4 mm (b) 8 mm

Fig. 5:	Axial liquid velocity profiles for columns of nozzle diameters (c) 12 mm and (d) sparger

Fig. 6:	Velocity distribution of air from sparger

Fig. 7:	Turbulent kinetic energy for columns of nozzle diameters (a) 4 mm (b) 8 mm

Fig. 8:	Turbulent kinetic energy for columns of nozzle diameters (a) 12 mm (b) sparger

The turbulence microscale governs the turbulent kinetic energy based on 2-phase k-ε formulation. They also reported that the 2D axisymmetric simulation indicates all the turbulence time scales smaller than the total averaging time to contribute to the turbulent kinetic energy at the microscale.

The kinetic energy distribution also depends on the bubble mean fluctuating velocity. The mean horizontal fluctuation of the bubbles is considerably higher than the vertical component for lower gas hold-up and is a result of the zigzag or helical motion of particular size bubble specially ellipsoidal bubble (Sommerfeld and Broder, 2009). The horizontal 2 components of the mean fluctuating velocity of the liquid phase are lower than the vertical one. This results in an isotropic turbulence in the bubble column. The higher vertical component is the result of the flow being mainly driven by the dominant vertical bubble rise due to buoyancy. Sommerfeld and Broder (2009) stated that the average turbulent kinetic energy first increases from the distributor to a maximum at about one-third of the column height and then continuously decreases toward the surface of the bubble column.

The turbulent kinetic energy decreases toward the surface because of dispersion of bubbles over the entire cross-section of the bubble column. As a result a smaller momentum transfer occurs from the bubble to the liquid. The level of turbulence increases with the increase in gas flow rate due to increase in agitation by the bubble. The increase in gas flow rate increase in bubble size which further enhance the agitation. They also claimed that for the larger bubbles the turbulent kinetic energy is unusually higher compared to the cases with smaller bubble size. Also higher level of bubble fluctuation energy exhibits due to increase of gas flow rate. Therefore, the kinetic energy distribution not only depends on the bubble size distribution but also distribution of gas holdup. Akhtar et al. (2007) studied the effect of distributor hole size on bubble size distribution and hydrodynamics at constant superficial gas velocity of 1.0 cm sec^-1.

They have considered 3 distributors with hole size ranging from 20-100 mm for 2 and 3 dimensional simulations. They reported that decrease in rise velocity with decrease in hole-size which becomes apparent by comparing bubble rise trajectories for 50 and 20 mm size distributor. Based on this they concluded that the bubble size has similar dependency on hole size as superficial gas velocity.

They stated that small bubbles having low rise velocity usually formed from small size distributor. They also observed that relatively bigger size bubbles are formed with 100 mm hole when compared with 50 and 20 mm with their less rise time to the same height of column. They showed from their simulation that at relatively low superficial bubbles exhibited rectilinear trajectory. They claimed that their results are in well agreement with the research of Glover and Generalis (2004) for identical problem by using algebraic slip mixture model.

CONCLUSION

In summary, the importance of the computational fluid dynamics in predicting flows provides good description about the column hydrodynamics in steady state. Here in this study, the hydrodynamics like gas holdup, axial liquid velocity and the kinetic energy distribution enunciated with the geometric parameters, different nozzle diameters and with spargers. The sparger provided good results with high mixing of gas and liquid. To accurately model the motion of gas and liquid phases in bubble columns, the prediction of optimum nozzle diameter is a must. Further, the nozzle diameter controls the bubble diameter and if the studies of nozzle diameter and the bubble coalescence and break-up theories when combined would provide highly accurate results with mass transfer and heat transfer rates because the effective surface area depends on diameter of bubbles.

The present research would provide a basis for such further research. The accuracy of the solution is enhanced by taking finer mesh elements in the regions closer to wall. It is also recommended to take higher resolutions of the mesh in the region where there are large gradients of gas velocities. The gas holdup and the axial liquid velocity are the background of enhancing the turbulence in the bubble column which gives the average kinetic energy distribution in the column. The distribution of kinetic energy is due to the bubble distribution along the axis of the column and the bubble fluctuation energy around the bubble with liquid due to the various lateral forces acting on it. The distributor size which causes the intensity of homogeneity and heterogeneity effects the size distribution of the bubble around the distributor. Also the radial distribution of the gas hold up and the axial liquid velocity depends on the column diameter as well as the physical properties of the liquid.

NOMENCLATURE