ENGINEERING DEVELOPMENT OF

SLURRY BUBBLE COLUMN REACTOR

(SBCR) TECHNOLOGY

Contract No. DE-FC-22-95 PC 95051
Monthly Report, Budget Year 6
Reporting Period: July 1-31, 2001
(For the 26th Quarterly Period: July 1 to September 30, 2001)
from
Chemical Reaction Engineering Laboratory (CREL)
Washington University

TO:Dr. Bernard Toseland

DOE Contract Program Manager

Air Products and Chemicals, Inc.

P. O. Box 25780

Lehigh Valley, PA 18007

FROM:Dr. Milorad P. Dudukovic'

The Laura and William Jens Professor and Chair

Director, Chemical Reaction Engineering Laboratory

Washington University

One Brookings Drive

Campus Box 1198

St. Louis, MO 63130

Cc:R. Klippstein, Air Products and Chemicals, Inc.

M. Phillips, Air Products and Chemicals, Inc.

L. S. Fan, Ohio State University

K. Shollenberger, Sandia National Laboratory

1

BUBBLE DYNAMICS IN BUBBLE COLUMNS

INTRODUCTION

During the month of July 2001 work was in progress on a variety of topics in order to complete the tasks set for the current budget period. In this monthly report we address research progress made in:

Computational fluid dynamics (CFD) of bubble column flows, where we examined the effect of various closures for interfacial terms on the computed flow fields;

Radioactive particle tracking in providing data at high gas velocities to assess the effect of operating conditions and of gas distributor on the flow field;

Calibration of the optical probe for bubble size distribution measurements;

Bubble population balance modeling, which is undertaken to predict the mean bubble size and remove the uncertainty regarding bubble size that was present in our previous CFD models.

The accomplishments in the above four areas are briefly outlined below.

1.EFFECT OF CLOSURES IN CFD MODELS ON BUBBLE COLUMN FLOWS

Highlights

  • Various closures for interfacial drag and virtual mass force used in the Euler-Euler simulation of two phase gas liquid bubby flows have been identified and coded for incorporation in the CFDLib code.
  • The mesh system has been developed for simulation of a two dimensional bubble column at different operating conditions for which experimental data exists.
  • CFDLIB Euler-Euler simulations have been conducted using different expressions for the drag and virtual mass force. The effect of these terms on the computed time averaged velocity profiles will be examined and reported next month.

Background

The bubble-driven flow in a bubble column is a typical two-phase system in which the flow is driven by buoyancy. From the viewpoint of physical modeling, two-phase bubbly flow is commonly defined as a flow pattern in which the gas phase is distributed within a liquid continuum in the form of discrete bubbles much smaller than the characteristic dimension of the container (column diameter). Numerical simulation of such two-phase flows either treats the system in an average sense as interpenetrating continua, or explicitly moves the discrete bubbles due to the body forces and interactions arising with the continuous fluid. In some cases feedback of the effect of the bubbles on the carrier fluid motion is modeled through various means. In the first approach, called Eulerian-Eulerian method, the two-fluid model is developed (Drew (1983)) to describe the motion for each of the phases in the Eulerian frame of reference. In the second approach, called Eulerian-Lagrangian method, the carrier phase is still described in the Eulerian frame while the discrete phase is treated as is, that is, each discrete particle (bubble) is tracked by solving its equation of motion in the continuous phase.

In order to simulate gas-liquid bubble column flows in large columns, usually a Eulerian-Eulerian approach is applied (Sokolichin and Eigenberger (1994), Pan et al., (1999) and Pan et al., (2000)). The main advantage of this approach resides in the fact that it is computationally less expensive as compared to Eulerian-Lagrangean approach or Direct Numerical Simulations. However, this approach requires specific consideration of the interface exchange terms in momentum balance equations, which arise due to the averaging process. The physical phenomena contributing to the momentum exchange at the interface include: (i) the momentum exchange due to drag between the phases, (ii) a virtual force experienced by the discrete phase due to the acceleration or deceleration of the carrier phase, (iii) a lateral force caused by the presence of vorticity in the continuous phase, (iv) force due to the pressure gradient, (v) force due to the history effects and so on. The major contributions to the inter-phase momentum exchange come from the drag force, virtual mass force and the lift force. As the inter-phase exchange phenomenon is a function of gas holdup, relative velocity between the phases, bubble size and fluid properties, the mathematical modeling of these forces must take into account these factors. Various studies reported in the literature use different drag coefficient and virtual mass coefficient correlations that take into account these parameters (e.g., gas holdup etc.) to a different extent. Most studies claim some limited agreement with the experimental data. However a systematic study that compares the various correlations and their effect on the predicted flow quantities with the experimental data is still lacking. For scale-up purposes of bubble column or slurry bubble column reactors, it is necessary to identify the sets of correlations which can be used in a wide range of parameters. This study is aimed to fill this gap and provide a comparison of the most commonly used correlations with each other and with the experimental data. There are several empirical/semi-empirical correlations available in the literature (Fan et al. (1990), Sankaranarayanan et al., (2000), Biesheuvel et al. (1989) and Felderhof (1991)) for the closure of different interface terms. The range of validity of each of these correlations depends not only on the flow regimes but also strongly on the flow parameters like: gas holdup, bobble size, superficial gas velocity etc.

In this study, using Eulerian-Eulerian approach, the mass and momentum balance equations are solved in the framework of Los Alamos National Laboratory’s code CFDLIB. Different correlations for drag force and virtual mass force are implemented in the code library. The average flow characteristics are obtained by ensemble averaging of instantaneous flow variables.

Ensemble averaged equations for two-phase flow

The ensemble-averaged equations for incompressible two-phase flow are given by Drew (1983). The continuity equations are:

/ (1)
/ (2)

The momentum equations for each phase can be written as:

/ (3)
/ (4)

Where l and g are the liquid and gas phase stress tensors.

The tensor represents the stress inside the bubble. One can ignore the stress inside a bubble due to small gas viscosity as compared to the liquid. The term() is the stress experienced by the liquid (gas) phase due to the presence of the gas (liquid) phase. In the case of very low gas superficial velocity (low gas holdup), which is the case in the bubbly flow regime, the direct contributions due to these terms can be ignored. However, the presence of the second phase can be taken into account by introducing the concept of Einstein’s viscosity correction for dilute suspensions, which implies that , where m is the modified viscosity of the first phase containing a dilute suspension of second phase. The turbulence in the liquid phase is taken care of using the single-phase k- model.

Interfacial Closures

The inter-phase momentum exchange due to drag force, Md and due to added mass force, Mvm are defined as:

/ (5)
/ (6)

For qualitative analysis, it is necessary not only to derive the correct form of these equations but also to obtain reliable estimates of the various averaged quantities, such as the viscous drag and added mass coefficients. The drag coefficient CD depends on the flow regime, gas holdup and the liquid properties.

Drag Force

The drag force arises in the multiphase flows due to the friction between the phases. Most of the studies reported in the literature consider the drag force as a force exerted on a single spherical particle moving in a large expanse of a continuous fluid at rest (see Equation 5). However, in real systems as the dispersed phase volume fraction increases from zero, the multi-body effects increase drag. Different correlations reported in the literature take account of this fact to different extents. The different drag force formulations used in this study are listed in Table 1. The correlations listed in Table 1 for Dr=1, 2 and 4 are based on the single sphere concept (of a mean diameter dp) and do not take into account the gas holdup while for Dr=3, 5 and 6 include the effect of local gas holdup.

Table 1: Drag coefficient correlations

Dr=1 / , where (CFDLIB)
Dr=2 / ; Schiller and Naumann (1935)
Dr=3 / , where ,; where Eo is the Eötvös number. Tsuchiya et al. (1997)
Dr=4 / Morsi and Alexander (1972)

Dr=5 / , Sankaranarayanan et al. (2000)
Assuming , one gets where is the Richardson-Zaki exponent, , , Re and Eo are bubble Reynolds number and Eötvös number, respectively, as defined earlier. Mo is the Morton number and  is the surface tension while is the bubble Reynolds number based on the terminal velocity of an isolated bubble given by Fan and Tsuchiya (1990).
Dr=6 / The drag coefficient Dr=1 with Richardson-Zaki correction for multi-body effects, that is

Virtual Mass Force

In unsteady flows, the inertia of the fluid effects the acceleration of the particle, which gives rise to the concept of added/virtual mass force. This force depends on the flow regime and the fluid properties. The virtual mass force correlations analyzed in this study are listed in Table 2.

Table 2: The virtual mass coefficient correlations

Vm=1 / (CFDLIB)
Vm=2 / ;
Biesheuvel and Spoelstra (1989)
Vm=3 / ; Sankaranarayanan et al. (2000)

Numerical Simulation

The flow dynamics in a two-dimensional bubble column (15.2x170 cm diameter x height) is chosen for simulation using different interfacial closures specified in Table 1 and Table 2. The liquid phase (water) is in batch mode with initial static height of 152 cm (static height/diameter = 10). The gas phase (air) enters the column through three nozzles (each of size 0.1 cm) at a superficial velocity of 1 cm/s over the section. This system is chosen for simulation for several reasons. First, extensive experimental data exists for the two dimensional column described above and two papers have treated this column already (Mudde et al., 1997; Pan et al., 2000) thus providing good background for comparison. Second, two dimensional computation can be executed fast so that feedback on the effect of various closures can be provided rapidly in contrast to three-dimensional simulations that may require a month to execute.

The conservation equations (1) through (4) are solved in the framework of Los Alamos National Laboratory code CFDLIB (Kashiwa and Rauenzahn, 1994) along with different closures for the drag force (see Table: 1) and the virtual mass force (see Table: 2). CFDLIB is a Fortran code library for multiphase flow simulations. It uses a cell-centered finite-volume method applied to the time-dependent equations. The methodology and the numerical scheme are described elsewhere (Cranfill 1983 and Kashiwa et al. 1993). To adapt the code for the present study, some parts of the code related to the inter-phase momentum exchange are modified according to the interfacial closures discussed in the previous section.

Figure 1 shows the mesh system for the column. The horizontal grid consists of 36 points, distributed as 6, 2, 10, 2, 10, 2 and 6 over the segments of 0.4, 0.1, 0.49, 0.1, 0.49, 0.1 and 0.4 cm. The vertical grid consists of 170 cells each of 1cm. At the inlet the gas phase has a fixed gas superficial velocity while the liquid phase has no-slip conditions. Both phases obey the no-slip velocity conditions at the vertical walls and the outflow boundary condition at the top of the column. All the simulations start from a static initial condition where the column is filled with liquid only up to the static height. Above the static height, there is gas only. The simulations are then performed until a quasi-steady-state is reached. Once the quasi-steady-state is reached, the velocity and gas/liquid holdup fields are stored for next 110 seconds with a time step of 0.1 second.

Figure 1: Geometry and Mesh System

The computations have been completed. The results are being analyzed and will be reported in the next monthly report.

References:

Bel F’Dhila R. and Simonin O., 1992, “Eulerian prediction of a turbulent bubbly flow down-stream a sudden pipe expansion”, 6th workshop on two-phase flow predictions, Erlangen, Germany.

Biesheuvel A. and Spoelstra S., 1989, “The added mass coefficient of a dispersion of spherical gas bubbles in liquid”, Int. J. Multiphase Flows, 15, 911-924.

Cranfill C. W., 1983, “EOSPAC: A subroutine package for accessing the Los Alamos Sesame EOS Data Library”, Los Alamos National Laboratory Report LA-9728-M, UC-32, Los Alamos, NM.

Drew D. A., 1983, “Mathematical modeling of two-phase flow”, Ann. Rev, Fluid Mech., 15.

Fan L. –S. and Tsuchiya K., 1990, “Bubble wake dynamics in liquids and liquid-solid suspensions”, Butterworth-Heinemann, Stoneham, MA.

Felderhof B. U., 1991, “Virtual mass and drag in two-phase flow”, J. Fluid Mech., 25, 177-196.

Kashiwa B. A., Padial N. T., Rauenzahn R. M. and Vanden-Heyden W. B., 1993, “A cell centered ICE method for multiphase flow simulations”, Los Alamos National Laboratory Report LA-UR-93-3922.

Krishna R., Urseanu M. I., Van Baten J. M. and Ellenberger J., 1999, “Rise velocity of a swarm of large bubbles in liquids”, Chem. Eng. Sci., 54, 171-183.

Mudde R. F., Lee D. J., Reese J. and Fan L, -S, 1997, “Role of coherent structures on Reynolds stresses in a 2-D bubble column”, AIChe. J., 43.

Pan Y., Dudukovic M. P. and Chang M., 1999, “Dynamic simulation of bubbly flow in bubble columns”, Chem. Eng. Sci., 54, 2481-2489.

Pan Y., Dudukovic M. P. and Chang M., 2000, “Numerical investigation of gas-driven flow in 2D bubble columns”, Fluid Mech. Transport Phenomena, Vol. 46, No. 3, 434-449.

Sankaranarayanan K., Shan X., Kevrekidis I. G. and Sundaresan S., 2000, “Analysis of drag and virtual mass forces in bubbly suspensions using an implicit formulation of the Lattice Boltzmann Method”.

Schiller L. and Naumann Z., 1935, Z. Ver. Deutsch. Ing., 77:318, 1935

2.CARPT DATA PROCESSING

Processing of CARPT data to obtain ensemble averaged liquid veloicities and various measures of liquid mixing is still in progress. The conditions summarized below in Table 3 have been covered by the CARPT runs which were completed in June 2001.

Table 3.Conditions of Recent CARPT Experiments

No. / Pressure (bar) / Superficial Gas Velocity (Ug (cm/s) /
Distributor Used
1 / 1 / 30 / D2, Cross Sparger (4 holes, 0.1 % porosity)
2 / 1 / 30 / D4, Perforated Plate (163 holes, 0.15 % porosity)
3 / 1 / 30 / D6, Perforated Plate (163 holes, 1 % porosity)
4 / 4 / 30 / D3, Single Nozzle (0.1 % porosity)
5 / 4 / 30 / D4, Perforated Plate (163 holes, 0.15 % porosity)
6 / 4 / 30 / D6, Perforated Plate (163 holes, 1 % porosity)
7 / 10 / 30 / D6, Perforated Plate (163 holes, 1 % porosity)
8 / 1 / 45 / D4, Perforated Plate (163 holes, 0.15 % porosity)
9 / 4 / 45 / D4, Perforated Plate (163 holes, 0.15 % porosity)
The designations D1, D2, etc. for the distributors refer to Figure 7 of the previous monthly report

The processed data from the above runs should enable us to assess the effects of i) distributor type, ii) pressure and iii) superficial gas velocity on liquid recirculation rate and liquid mixing and turbulence.

3.CALIBRATION OF OPTICAL PROBE

Our work with the four point optical probe continues. It was determined (see Table 4 below) that front edge of the liquid provides better estimates of bubble rise velocity. Further intricate details of the probe operation are being carefully examined.

Table 4.Comparison of Optical Probe Determined and Camera Determined Bubble Rise Velocity

Camera / Probe
(front edge method) / Probe
(back edge method)
Rise velocity (cm/s) / 26.5 / 26.8 / 30.8
Relative Error / 1.2% / 16.2%

4.INCORPORATION OF THE POPULATION BALANCE INTO THE CFD MODEL

Highlights:

The effect on the computed flow field of preserving either the surface balance or the number balance when reassigning bubbles in a numerically simulated coalescense – redispersion process to the sizes whose evolution is followed by the code is minimal. Hence, either of the above choices can be made in the code.

Explanation:

In the previous monthly report, we implemented the bubble population balance equation (7) into the CFD framework using FLUENT.

(7)

The basic idea when discretizating the above equation is that bubbles in a size range, say Ri, are assigned to a pivotal size xi. However, breakup and coalescence processes may produce particles that are between such pivotal sizes (except in the case of a uniform linear grid, i.e., xi = ivmin) and must be reassigned to the pivots. The reassignment must be done carefully to preserve the accurate calculation of selected moments of the probability density function.

Kumar and Ramkrishna(1996a) proposed the following way to preserve any selected moments. For xi v xi+1, let the fraction of bubbles of size v assigned to xi be denoted by , and a fraction of be assigned to size xi+1. The reassignment will preserve the rth moment provided

for r = r1, r2(2)

These two equations above (i.e. Eq. 8 for r = r1 and r = r2) yield a unique solution for the quantity . In order to preserve the mass balance, we must set r1 = 1. Then we can choose to preserve either the surface (r2 = 2/3) or number (r2 = 0) of bubbles upon reassignment since we must maintain the mass balance (r1 =1). In the present work, we studied the effect of the choice of preserving either the number density (r2 =0) or of preserving the surface (r2=2/3).

The flow in bubble columns was modeled using the Algebraic Slip Mixture model (ASMM) incorporated in the FLUENT software. The breakup kernel as given by Martínez-Bazán et al. (1999) and coalescence kernel as given by Luo (1993) are used, respectively. Nine-class bubble population was considered in the simulation.

In Figure 2 the computed time-averaged liquid axial velocity profiles are compared against the time-averaged velocity profiles obtained by CARPT. The superficial velocity is Ug = 12 cm/sec and the column diameter is 8 inch. The results presented are at a height of 60 cm above the sparger, typical of the fully developed region of the flow. The centerline axial velocity is somewhat overpredicted by all simulations (within 14%). However, the general shape of the profile is well captured and the discrepancy in the model predicted and CARPT-measured time-averaged axial liquid velocity diminishes as one moves radially outwards in the column. Moreover, the crossover point (i.e., the radial location where the axial velocity component becomes zero) is well predicted. The difference between the two predicted results is much smaller than the difference between the CARPT data and the liquid velocity prediction.