## Info

where kL and kT are characteristic laminar and turbulent resistance coefficients of the sparger, and pS is the pressure below the sparger (plenum pressure).

For column walls, which are impermeable to fluids, standard wall boundary conditions may be specified. Whether the full column should be considered in the solution domain or symmetry or other boundary conditions may be invoked to reduce the extent of the solution domain, essentially depends on the objective and the proposed mathematical model. If the interest is in estimating long-time-averaged flow characteristics (as done by Ranade, 1997), invoking symmetry is often useful and can facilitate rapid results. However, when the interest is in capturing inherently unsteady flow characteristics, which are not symmetrical, it is essential to consider the whole column as the solution domain (as done by Ranade and Tayalia, 2001). Overall flow can be modeled using an axis-symmetric assumption, if and only if, the governing equations are derived in such a way as to represent the influence of local dynamic flow structures on time-averaged flow via additional terms. Ranade (1997) developed such a model based on two empirical coefficients to simulate long-term-averaged flow patterns. The two parameters were estimated from a data set of Yao et al. (1991). The same values of these parameters were used for subsequent simulations. Results show good agreement between predicted results and long-time-averaged experimental results of Hillmer et al. (1994), Grienberger and Hofmann (1992), Menzel et al. (1990) and Hills (1974). However, when the objective is to simulate inherently unsteady flow, which is asymmetric, it will be necessary to avoid imposing symmetry boundary condition. For a two-dimensional column, therefore, it will be necessary to consider the whole domain. For three-dimensional cylindrical columns, either a body-fitted grid should be used, which will avoid specification of boundary condition at the axis, or conventional axis boundary conditions should be modified to allow flow through the axis (the value of any variable at the axis location is set to average of all the computational cells surrounding it).

Boundary conditions for the top horizontal surface of the bubble column requires special attention. Following the experimental practice of keeping the height of the gas-liquid dispersion smaller than actual column height, the top surface of the column can be modeled as an outlet for gas and liquid phase. It is expected that the solution of the model equations will determine the height of the gas-liquid dispersion and only gas will exit from the inlet. In other words, it is expected that the model will predict a gas volume fraction of unity above the gas-liquid interface. Such an approach, however, has to ensure that the governing equations are capable of handling change in the prevailing continuous phase (liquid below the gas-liquid dispersion height and gas above it). This is seldom done and despite this, solution of the conventional two-fluid model using the top surface as an outlet is attempted. As expected, such an attempt always leads to non-physical velocity values at the region of the gas-liquid interface and encounters severe convergence difficulties. In many cases with high superficial gas velocity, the liquid mass of the column is lifted out of the column. In order to enhance convergence behavior, Padial et al. (2000) used a smoothly varying continuous phase density across the gas-liquid interface:

Such an empirically adjusted smooth profile of continuous phase density helps to maintain a fairly stable interface. Some fraction of the liquid, however, may still escape the column during initial phases of solution development.

In most reactor engineering applications, it may not be necessary to include the gas-liquid interface in the solution domain. In an alternative approach, the solution domain is restricted to the height of gas-liquid dispersion. Of course, an exact value of dispersion height is not known a priori. However, in most cases, overall gas volume fractions and therefore, height of the gas-liquid dispersion can be estimated. Even if there is 40% error in the prediction of overall gas volume fraction, it will result in only 10% error in the estimation of height of gas-liquid dispersion, if the volume fraction is about 25%. Except for very shallow bubble columns, fluid dynamics is not sensitive to the 10% error in the height of gas-liquid dispersion. Thus, it is always possible to estimate a reasonable solution domain height to model gas-liquid flows in bubble columns. The top surface of the solution domain may then be assumed to coincide with the free surface of dispersion. This free surface may or may not be assumed to be flat. The normal liquid phase velocity, the tangential shear stress and the normal fluxes k, e and <p are set to zero at the free surface. The gas bubbles are free to escape from the top surface. If source code is accessible, one can modify the code to implement these boundary conditions. In most commercial CFD codes, user-defined routines with options such as 'patch boundary conditions' may be used to implement this boundary condition. When direct implementation is not possible, Ranade (1998, 2000) proposed two approximate alternatives.

In the first alternative, if the terminal rise velocity of gas bubbles is known (or can be estimated with confidence), the top surface of the dispersion may be defined as an 'inlet'. Normal liquid velocity may be set to zero while normal gas velocity may be set to terminal rise velocity. The implicit assumption here is that gas bubbles escape the dispersion with terminal rise velocity. It should be noted that even after defining the top surface as an inlet, gas volume fraction at the top surface is a free variable. There is no implicit forcing of gas volume fraction distribution. Alternatively, the top surface of the dispersion can be modeled as a no shear wall. This will automatically set normal liquid velocity to zero. It will also set normal gas velocity to zero. In order to represent escaping gas bubbles, an appropriate sink may be defined for all the computational cells attached to the top surface (Figure 11.7):

where AB is the area of the bottom surface of the computational cell attached to the top surface, WGB and aGB are the normal velocity of gas bubbles and gas volume fraction of the computational cells lying below the computational cell attached to the top surface. Such formulations of top surface avoid handling sharp gradients of gas volume fractions at the gas-liquid interface and are much more stable numerically.

To illustrate application of the Eulerian-Eulerian approach, some results of two-dimensional bubble columns are discussed here. Three-dimensional bubble columns and other reactor engineering applications are discussed in Section 11.3.

Gas-liquid Interface

No shear wall

Gas-liquid Interface

Top row of ■ computational cells

Define sink for gas phase at the top row using Eq. (11.24)

FIGURE 11.7 Top boundary condition for bubble column reactor.

Brief review of recent simulations of bubble columns

Rigorous experimental data is available for apparently two-dimensional bubble columns (Becker et al., 1994, 1999). It is, therefore, useful to simulate the fluid dynamics of such a column to validate the underlying mathematical model. Pfleger et al. (1999) recently reported comparison of their simulated results with experimental data (Fig. 11.3). Sokolichin and Eigenberger (1999) also reported good agreement between predicted results and experimental data. To illustrate such simulations, we reproduce some of the results obtained by Ranade (2000). He simulated gas-liquid flow in a two-dimensional bubble column (0.2 m width, 0.45 m dispersion height and 0.04 m depth) having the same geometrical column configuration as used in the experiments by Becker et al. (1999). A solution domain and computational grid (51 x 90 x 11 computational cells for width x height x depth) is shown in Fig. 11.8. A two-fluid model was used to simulate gas-liquid flow in such a column. A standard k-e model was used to simulate turbulence. A QUICK discretization scheme was used with a SUPERBEE limiter functions (see Chapter 6 for more details on discretization schemes). Modeled interphase coupling terms comprised drag force and virtual mass terms. For estimation of drag coefficient, a correlation of Schiller and Naumann (1935) was used. Governing equations were solved in a transient manner with a time step of 0.008 s. For each time step, 100 internal iterations were carried

FIGURE 11.8 Solution domain and grid.

out. This number of internal iterations was sufficient to ensure an adequate degree of convergence of all the equations. Further refinement of time step did not affect the simulation results significantly.

Simulation results clearly showed the inherently unsteady flow characteristics of gas-liquid flows. The predicted instantaneous flow field is shown in Figs 11.9 and 11.10 at three different planes. It can be seen that gas bubbles rise in the column in a meandering way. Experiments by Becker et al. (1999) indicate that the bubble swarm moves laterally with a period of about 16 s. Eulerian-Eulerian simulations carried out by Ranade (2000) also show a meandering effect. The time history of the simulated velocity field at four locations within the column clearly shows the oscillatory nature of the flow (Fig. 11.11). It can be seen that simulation results somewhat underpredict the period of oscillations (~ 12 s). However, overall characteristics are in good agreement with experimental observations. Recent simulations by Pfleger et al. (1999) also show some underprediction of the period of oscillation. Despite this, time-averaged predictions agree quite well with experimental data. The predicted results shown in Figs 11.9 and 11.10 indicate the three-dimensional nature of the flow even in

FIGURE 11.10 Predicted gas volume fraction distribution (from Ranade, 2000).

an apparent two-dimensional bubble column. If such three-dimensional effects are not considered, the meandering motion of the bubble swarm is not captured with adequate accuracy because of the overprediction of turbulent kinetic energy. Volume integration of the predicted turbulent energy dissipation rate indicates that about 30% of the input energy (g(Vo}) is dissipated in the form of turbulent energy. The remainder of the energy must be dissipating at the gas-liquid interface. Predicted overall gas volume fraction is 0.63%. If we assume that average slip velocity is about 0.23 ms-1, the amount of energy dissipated at the gas-liquid interface may be estimated as the product of the gravitational constant, gas volume fraction and slip velocity. This is about 68% of the total input energy. This value agrees quite well with the predicted overall turbulent energy dissipation (assuming that the laminar dissipation is negligible compared to turbulent dissipation). The simulated results can be examined in a variety of ways to understand the dynamic characteristics of gas-liquid flows in bubble columns. For example, three-dimensional iso-surface plots

FIGURE 11.11 Transient velocity traces at different locations in rectangular 2D column (from Ranade, 2000).

may give more information about the prevalent flow structures. The iso-surface of gas volume fraction may also clearly show the meandering motion of the bubble swarm (not shown here). In order to identify streamwise vortices, it is useful to examine values of normalized helicity (which is a scalar quantity representing the cosine of the angle between the velocity and vorticity vectors). Near the vortex cores, the magnitude of normalized helicity will be near unity with its sign depending on the orientation of the velocity vector to the vorticity vector. A typical helicity iso-surface plot is shown in Fig. 11.12 and shows several regions high magnitude normalized helicity indicating the presence of vortex cores. Thus, Eulerian-Eulerian simulations are able to capture the inherently unsteady flow characteristics of gas-liquid flows in bubble columns.

For bubble columns with relatively low gas volume fractions, bubble size distribution is fairly narrow. As gas velocity and therefore, gas volume fraction increases, a heterogeneous or churn-turbulent regime sets in with much wider bubble size distribution than the homogeneous regime. With such a wide bubble size distribution, it is important to develop appropriate averaging methods and corresponding closure models. Use of governing equations, derived based on the assumption of a single bubble size, generally lead to significant overprediction of gas volume fraction, though comparison of liquid phase mean velocity is not bad (see Kumar et al., 1994, for example). Recently, Krishna et al. (2000a) proposed use of a three-phase model to simulate the churn-turbulent regime of bubble columns. Basic concepts are shown schematically in Fig. 11.13. Based on experimental observations, the gas phase is divided into two separate phases, one with small bubbles of diameter of the order of a few millimeters, and the other with large spherical cap bubbles with diameter of a few centimeters. Large bubbles were introduced only in the central core, in their simulations. Drag coefficient values were calculated based on empirical correlations of observed bubble

FIGURE 11.12 Iso-surface of normalized helicity. Normalized helicity = 0.5 (dark surface) and -0.5 (light surface) (from Ranade, 2000).

slip velocities. Krishna et al. (2000a) had introduced four adjustable parameters to tune the drag coefficient and diameter of the large bubbles. With these parameters, they were successful in simulating liquid velocity profiles and gas volume fractions as observed by Hills (1974). Treatment of the outlet boundary conditions was not discussed in their paper. Comparison of their simulated results and experimental data is shown in Fig. 11.14. The comparison looks quite adequate for most engineering applications. It must, however, be noted that these results may not be grid independent and the values of adjustable parameters may depend on grid size.

Instead of arbitrarily considering two bubble classes, it may be useful to incorporate a coalescence break-up model based on the population balance framework in the CFD model (see for example, Carrica et al., 1999). Such a model will simulate the evolution of bubble size distribution within the column and will be a logical extension of previously discussed models to simulate flow in bubble columns with wide bubble size distribution. Incorporation of coalescence break-up models, however, increases computational requirements by an order of magnitude. For example, a two-fluid model with a single bubble size generally requires solution of ten equations (six momentum, pressure, dispersed phase continuity and two turbulence characteristics). A ten-bubble class model requires solution of 46 (33 momentum, pressure,

Phase 1:

liquid

FIGURE 11.13 Three phase simulation of bubble columns (from Krishna et al., 2000a).

Phase 1:

liquid

FIGURE 11.13 Three phase simulation of bubble columns (from Krishna et al., 2000a).

10 continuity and two turbulence characteristics) equations. In such a model, each bubble class is considered as a separate phase and momentum and continuity equations are solved for each bubble class. If the bubble size distribution is relatively narrow and can be represented by a single velocity field, it is possible to reduce the computational requirements by solving the momentum equations for only one bubble class and assigning the same velocity field to all bubble classes (Lo, 2000). Separate continuity equations comprising terms pertaining to coalescence and break-up need to be solved for each bubble class. For such a case with the ten-bubble class model, the number of governing equations reduces from 46 to 19. This approach will be quite useful to simulate gas-liquid interfacial area and mass transfer with moderate increase in computational demands. Recently, Buwa and Ranade (2000) employed this approach to simulate evolution of bubble size distribution in a two-dimensional bubble column. Their model is briefly discussed in Appendix 11.1.

Bubble column reactors are also widely used as slurry reactors. The simplest way to simulate the fluid dynamics of slurry bubble column reactors is to treat the solid and liquid phases as a single slurry phase with effective slurry properties. This approach was used by Torvik and Svendsen (1990) and may give reasonable predictions for low solid volume fractions. The variation of solid volume fraction within the column is ignored in such an approach. For denser slurry applications, however, solid volume fraction may vary significantly within the column and solid and liquid phases need to be considered separately. The overall fluid dynamics can then be simulated using a three-fluid model (gas, liquid and solid). It is necessary to formulate appropriate interphase exchange terms for these three phases. When the solid volume fraction in slurry bubble column reactors increases even further and is more than 10%, it is necessary to use granular flow models to simulate flow of solid particles. Details of granular flow models are discussed in the next chapter when describing the modeling of fluidized

O Experimental data of Hills (1974) - 3D simulations ---2D axi-symmetric

(a) U=0.019m/s (b) U=0.038m/s (c) U=0.064m/s (d) U= 0.095 m/s (e) U=0.169 m/s

0 0