- Email: [email protected]
Opaque multiphase flows: experiments and modeling
Experimental Thermal and Fluid Science 26 (2002) 747–761 www.elsevier.com/locate/etfs
Opaque multiphase flows: experiments and modeling M.P. Dudukovic
*
The Chemical Reaction Engineering Laboratory (CREL), Washington University, One Brookings Drive, Campus Box 1198, St. Louis, MO 63130-4899, USA Accepted 28 November 2001
Abstract This brief review deals with the characterization of multiphase flows needed in reaction engineering and describes the techniques and methods available to accomplish it. It emphasizes the fact that most multiphase systems of industrial interest (e.g. stirred tanks, slurry bubble columns, gas–solid or liquid–solid risers and fluidized beds, pneumatic transport reactors, etc.) operate in churn turbulent flow at large volume fraction of the dispersed phase and are opaque. This inevitably leads to the use of the Euler–Euler, interpenetrating fluid model in describing their fluid dynamics, which raises issues of closure for the interfacial momentum terms and turbulence. Hence, computational fluid dynamic (CFD) models based on the Euler–Euler approach with selected closures require experimental validation. Due to the non-transparency of the systems involved, non-conventional experimental techniques are required. The paper briefly describes two such techniques which have been successfully implemented at the Chemical Reaction Engineering Laboratory for characterization of opaque, multiphase flows. One is the gamma ray computed tomography (CT), the other is the computer automated radioactive particle tracking (CARPT). Their use in characterizing diverse systems such as liquid–solid risers and bubble columns is illustrated as well as the utility of CARPT-CT data in validation of the Euler–Euler CFD models. The goal is to enable the reader to grasp the advantages and limitations of the described experimental techniques and computational fluid dynamic models. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Riser; Bubble column; Gas–liquid; Liquid–solid; Radioactive particle tracking; Computed tomography
1. Introduction Applications of multiphase technology are prevalent in petroleum refining, synthesis gas conversion to fuels and chemicals, bulk commodity chemicals, specialty chemicals, conversion of undesired byproducts to recyclable products, manufacture of polymers, biotechnology and pollution abatement. It is estimated that the value of shipments generated by industry that predominantly uses multiphase technology exceeds 640 billion dollars a year in the United States alone [1]. As a rule the multiphase reactor is the heart of any process in which chemical transformations take place, because reactor performance determines the number and size of needed separation units in front and after the reactor and, hence *
Corresponding author. Tel.: +1-314-935-6021; fax: +1-314-9354832. E-mail address: [email protected] (M.P. Dudukovic).
dictates the economics of the whole process. It is thus essential to quantify various measures of reactor performance, such as conversion, selectivity, production rate, etc., as a function of the design and operating parameters which requires appropriate descriptions of micro- and macro-transport. The discipline of reaction engineering attempts precisely to do that by quantifying the transport-kinetic interactions on various scales. Starting with species continuity equations and the energy balance (Fig. 1), one recognizes the need to describe the phenomena that occur in the reactor on a myriad of scales. Molecular scale knowledge allows us to write the kinetic rate expressions, R, based on first principles. Understanding the interaction of flow and transport in a single eddy or on a single catalyst particle scale enables us to modify the kinetic forms by local transport effects, g, and describe more precisely the source/sink terms in the species and energy balance equations as well as to describe transport between phases. The left hand side
0894-1777/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 8 9 4 - 1 7 7 7 ( 0 2 ) 0 0 1 8 5 - 1
748
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
Nomenclature C Drr Dzz E F g k khs ks Ksf p Q r R T t u, U z
concentration of scalar quantity (like tracer), kg/m3 radial eddy diffusivity, m2 /s axial eddy diffusivity, m2 /s E-curve, exit age density function, s1 force, N gravitational acceleration, m/s2 mass exchange coefficient, m/s coefficient of granulary conductivity, kg/m s kinetic energy per unit mass of solids, m2 /s2 solids–fluid exchange coefficient, kg/m3 s pressure, Pa volumetric flow rate, m3 /s radial position, m reaction rate, mol/m3 s temperature, K time, s velocity, m/s axial distance, m
Greek symbols DH heat of reaction, J/mol
chs e g h q s /fs
collisional dissipation term, kg/m s3 void fraction, dimensionless effectiveness factor, dimensionless granular temperature, m2 /s2 density, kg/m3 shear stress, N/m2 solids energy dissipation by correlations between solids and liquid velocity fluctuations, kg/m s3
Subscripts b bulk fluid f fluid phase g gas phase l liquid phase s solid phase k kth phase 1, 2 phase 1, 2 Other Symbols h i average 0 fluctuating quantity
Fig. 1. Reaction engineering models.
of the balance equations, containing the input–output terms due to convection, diffusion and interphase exchange, rests on the description of the flow field and its approximation. Rarely is in reactor engineering the momentum balance solved coupled with the species and energy balance to provide for the exact flow field and for the exact species spatial and temporal distribution. In recent years there has been an increased drive towards the use of fundamentals in reaction engineering, pushing the industrial practice from strictly empirical to increasingly fundamental, mechanism based, descriptions of kinetics. At the same time direct numerical simulations (DNS) allow at certain conditions an increasingly accurate description of flow fields and transport around small groups of particles helping steadily
advance the understanding on the single eddy/particle scale. In contrast, the description of reactor scale flows and mixing has not progressed much and in the industrial applications it is at best at the level of ideal reactor assumptions of plug flow or perfect mixing. In a plug flow reactor (PFR) a piston type one dimensional (1D) flow with flat velocity profile is assumed. In a perfectly mixed continuous flow stirred tank reactor (CSTR) it is assumed that the feed looses identity instantly and that all of the reactor content is at a uniform composition and temperature equal to those of the outlet stream. The problem of mixing and flow is avoided by assuming that the time constants for those phenomena are orders of magnitude shorter than the characteristic reaction time. This, however, is hardly ever the case! Typical models
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
749
Fig. 2. Flow pattern models for reactors with two moving phases.
for a two phase reactor with two flowing phases are sketched in Fig. 2. Either plug flow or perfect mixing is assumed in each phase. Plug flow and perfect mixing combination (not shown) is also used [2]. On occasion an empirical 1D axial dispersion model (ADM) is used in each phase to fill the gap between PFR and CSTR. The dispersion coefficient seldom can be scaled up. In each of these models the interphase exchange is lumped into an overall exchange coefficient that needs to be multiplied by an average driving force to describe the transport. A much more desirable approach would be to rely on phenomenological models that capture the essence of the flow field. Such phenomenological models can be developed based on experimental observations coupled with CFD model calculations. Ultimately, CFD calculations should be used to assess the flow field of each phase. Unfortunately, for multiphase flow, CFD calculations are still subject to uncertainty, especially since DNS cannot be used for simulation of large reactors. It seems that for simulation of pilot plant and commercial scale multiphase reactors the Euler–Euler k-fluid model is the best suited. This, of course, requires various closures for the interphase interaction forces and intraphase turbulence. Thus, the results produced by the multiphase CFD codes need to be experimentally validated. Since the systems to be interrogated are opaque, and we cannot ‘‘see’’ into them, most of the non-invasive techniques that are based on optical methods or use lasers cannot be employed. Fortunately, as recent reviews illustrate [3,4], there are techniques which can provide us with the desired information. In this paper the focus is on introducing two of them: gamma ray computer tomography (computed tomography, CT), for measurement of holdup profiles, and computer aided radioactive particle tracking (computer automated radioactive particle tracking, CARPT), for measurement
of velocity profiles and mixing parameters. We show how these techniques can be used to obtain information in systems of industrial interest such as a liquid–solid riser, or a gas–liquid bubble column. The ability of the available CFD Euler–Euler models to correctly calculate the observed hydrodynamic quantities is also briefly discussed.
2. Gamma ray computed tomography (CT) To measure the phase holdup distribution at any desired cross-section of the column we implemented a gamma ray source based fan beam type CT unit (Fig. 3). A collimated hard source (100 mCi of Cs-137) is positioned opposite to eleven 200 sodium iodide detectors in a fan beam arrangement. The lead collimators in front of the detectors have manufactured slits and the lead assembly can move so as to allow repeated use of the same detectors for additional projections. A 360° scan can be executed at any desired axial location and columns from 100 to 1800 in diameter can be scanned in the current configuration. The principle of computed tomography (CT) is simple. From the measured attenuation of the beams of radiation through the two phase mixture (projections) we calculate due to the different attenuation by each phase, the distribution of phases in the cross-section that was scanned. In our Chemical Reaction Engineering Laboratory (CREL) CT unit we achieve 3465–4000 projections and in the original unit we had a spatial resolution of 3 mm and density resolution of 0.02 g/cm2 . Of course we need a long time to scan the cross-section (about 45 min) and, hence, can only determine the time averaged density distribution. Details of the CT unit at CREL are described elsewhere [5–7]. Among the suggested techniques for reconstruction, e.g. convolution or filtered back projection, algebraic
750
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
Fig. 3. CT setup (a) schematic (b) unit in CREL.
reconstruction and estimation–maximization algorithms (E–M), we have found the E–M algorithm to be the best [6]. The CT hardware and software used in CREL and described above is quite suitable to provide us with reliable images of time averaged cross-sectional density distribution of the opaque system of interest at any desired elevation. It accomplishes that in a non-invasive manner and provides us with the ability to assess the spatial holdup distribution in two phase systems with good resolution compared to the vessel scale. It cannot, at present, yield ‘‘instantaneous images’’ and capture the dynamics of the system, nor can it provide very fine
resolution of the order of a single catalyst particle or a single small bubble. Since it is the reactor scale flow description that we are interested in, our CT unit provides valuable information for that purpose as illustrated above.
3. Computer automated radioactive particle tracking (CARPT) We have adopted CARPT (Fig. 4), the technique based on following the motion of a single tracer particle, as a method of mapping the velocity field in the whole
Fig. 4. CARPT setup (a) schematic (b) unit in CREL.
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
system. This Lagrangian technique was introduced for monitoring the motion of solids in fluidized beds by Lin et al. [8] and was adapted for measurement of liquid velocities in bubble columns by Devanathan and coworkers [9,10]. A single radioactive particle dynamically similar to the phase to be traced is introduced into the system. When solids motion is monitored in ebullated beds, fluidized beds or slurries, a particle of the same size and density as the solids in the system is prepared. For monitoring of liquid motion a neutrally buoyant particle is made [9]. For example, a polypropylene bead of about 2 mm is hollowed out, a small amount of Scandium 46 is inserted and a polypropylene plug is put in place, so that the density of the composite particle consisting of polypropylene, Scandium and air gap equals that of the liquid. Application of a thin film metallic coating assures that bubbles do not preferentially adhere to the particle. An array of scintillation detectors is located round the column. In our case up to 32 NaI 200 detectors are used. The detectors are calibrated in situ with the tracer particle to be used to get the counts–distance maps. Conventionally, CARPT calibration was done by positioning the tracer particle (usually Sc-46 of 250 lCi strength was used) at about 1000 known locations and recording the counts obtained at each detector. Various ways of doing this have been described [9,11]. An alternative way of constructing distance–count maps is via modeling of particle emission of photons and transmission and subsequent detection at the detectors. The Monte Carlo method [12,13] in which the photon histories are tracked in their flight from the source, through the attenuating medium and their final detection (or lack of it) at the detector can be used for this purpose. Thus, both the geometry and radiation effects are accounted for in estimation of the detector efficiencies in capturing and recording the photons. This involves calculation of complex three dimensional (3D) integrals which are evaluated using the Monte Carlo approach by sampling modeled photon histories over the main directions of their flight from the source. Once calibration is complete, the tracer particle is let loose in the system and the operating conditions are controlled and kept constant for many hours while the particle is tracked. A regression, least-squares method is used to evaluate the position of the particle. Sampling frequency is adjusted between 50 Hz and 500 Hz to assure the required accuracy. Typically, it is selected at 50 Hz for bubble columns since the finite size particle used as tracer cannot capture motion of frequencies above 25 Hz [11]. A wavelet based filtering algorithm [11] is employed to remove the noise in position readings created by the statistical nature of gamma radiation. The processing of data proceeds along the following lines. Filtered particle positions at each sampling time establish the tracer particle trajectory throughout the system. In systems like bubble columns with batch liquid the par-
751
ticle trajectory is confined to the liquid phase and the tracer particle stays in the vessel at all times. In risers, or bubble columns with flowing liquid, the tracer particle leaves and re-enters the system repeatedly. In either case from filtered particle positions at subsequent sampling times the instantaneous particle velocity is calculated and assigned to the column compartment (a grid is preestablished for the column) into which the mid-point falls. Thus, the instantaneous velocities are obtained repeatedly for each compartment. The run continues until adequate statistics are obtained. Then for each compartment (cell) ensemble average velocities are evaluated, which by ergodic hypothesis are assumed to be the time averaged velocities for each such location in the column. The differences between instantaneous and average velocities allow the estimation of the fluctuating velocities and all their auto- and cross-correlations leading to estimates of turbulent kinetic energy and Reynolds stresses. Most importantly, from the reaction engineering viewpoint, the appropriate Lagrangian correlations lead to the evaluation of the components of the eddy diffusivity tensor. The CARPT technique is very well suited for describing the reactor scale flows in opaque system by noninvasive means. The spatial resolution of 2–5 mm is quite adequate and frequencies up to 25 Hz can be captured with relatively large neutrally buoyant particle while higher frequencies are accessible with smaller particles. The technique cannot provide at present a fine spatial or temporal resolution. Its utility is illustrated below for liquid–solid and gas–liquid flows.
4. Liquid–solid riser flows Liquid–solid riser is one of the reactors of choice for alkylation reactions needed in production of higher octane fuels and linear alkylbenzene. This reactor with novel solid acid catalyst replaces the environmentally problematic processes that used HF and H2 SO4 as liquid catalysts. The solid acid catalyst, however, does deactivate and needs to be replaced. For that reason it is carried by liquid reactants through a riser, in which reaction takes place, and then is recycled through a regeneration unit. For proper reactor design and performance prediction it is necessary not only to know the solids holdup profile but also the solids mean velocity profile and measures of solids backmixing. Before our study the prevailing working assumption was that solids and liquid are in plug flow (piston flow). The average slip velocity between them defined, via Richardson– Zake correlation and the global drift flux model, the total solids holdup in the bed. This was insufficient information for proper analysis of reactor performance. We have provided now a more detailed description of the solids holdup and flow pattern [14,15].
752
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
4.1. Experimental results Our setup consisted of a 6 in. (15 cm) diameter 7 ft (213 cm) long plexiglass column. 2.5 mm glass beads of 2.5 density of water were used in water as the solid– liquid system. The solids mass flux was controlled by water flow through an eductor which was calibrated using a novel technique [16]. Tracer studies in the liquid phase (using KCl as tracer and conductivity probes) confirmed that liquid is in turbulent flow and close to plug flow [14]. The dimensionless variance of the measured liquid residence time distribution (RTD) indicates always more than 10 CSTRs in series at several superficial liquid velocities and three solids to liquid ratios. If liquid were flowing alone, the Reynolds number for the riser would have been between 23,000 and 35,000, so that the liquid flow is definitely turbulent. Solids holdup cross-sectional profiles were measured by gamma ray computer tomography using the CT equipment by Roy [14]. Relatively symmetric cross-sectional distributions are observed which allows us to azimuthally average the solids holdup. In Fig. 5 the radial profiles at three elevations and at three solids to liquid flux (S=L) ratios of 0.10, 0.15 and 0.20 at a liquid superficial velocity of 20 cm/s are presented. Clearly, the solids holdup profiles are quite similar at all elevations and the flow can be considered from the engineering standpoint as fully developed. At all conditions there are slightly more solids at the wall than in the middle of the riser but the radial profiles are not nearly as pronounced as in a gas–solid riser. The magnitude of measurement errors are indicated by bars in Fig. 5. Solids velocity was monitored by a radioactive tracer particle made of an aluminum ball of 2.5 mm in diameter with scandium inserted into it so as to match the density of the glass beads. The settling velocity of the tracer in water was determined to be 30:1 0:1 (cm/s) while the glass beads had a settling velocity of 30:2 0:1 cm/s. Sc-46 was activated to 290 lCi. The tracer particle behaves as the solids being traced. The sampling frequency was 50 Hz and the time of the tracer particle entry and exit of the riser was also carefully recorded providing the distribution of particle residence times and a probability density function (p.d.f.) of independent trajectories. The Lagrangian trace of a single particle trajectory during its sojourn through the riser is shown in Fig. 6 projected onto a horizontal plane and two vertical planes. The time series for the radial and axial particle position are also shown. Several thousands of such trajectories were sampled at each steady operating condition. Clearly the solid particle does not experience plug flow. From these trajectories the instantaneous and ultimately average tracer velocity is calculated at each location of the riser. The vector plots of the ensemble averaged solids velocity vectors projected on various
Fig. 5. Circumferentially averaged time averaged solids holdup distribution in the liquid–solid riser at U1 ¼ 15 cm/s. (a) S=L ¼ 0:10, (b) S=L ¼ 0:15, (c) S=L ¼ 0:20 (bars show standard deviation of azimuthal variation).
vertical planes, not shown here but available in Roy [14] and Roy and Dudukovic [15] reveal a fully developed flow pattern with solids rising in the middle and descending by the wall. Azimuthal averaging of these vector plots produces almost zero average solids radial velocities and radial profiles of the axial solids velocity that are almost independent of elevation from 50 to 150 cm at all solid to liquid (S=L) flux ratios studied as seen in Fig. 7. It is evident that the axial solids velocity profile in the radial direction is very pronounced with the solids in the time average sense rising in the middle of the riser and falling by the walls. Coupled with the
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
753
Fig. 6. Lagrangian trace of tracer particle during one sojourn in the riser (U1 ¼ 20 cm/s; S=L ¼ 0:15) and time series of radial and axial position.
solids holdup profile this provides the information on the extent of solids backflow. The combination of the solids holdup profile and axial velocity profile at all conditions satisfies the overall solids continuity usually within 5% with maximum discrepancy bounded by 10%. For reactor modeling we need to describe the extent of solids backmixing, i.e. departure from the assumed plug flow. The classical way of determining the extent of backmixing is to examine the RTD and its variance. In CARPT experiments the solids RTD is obtained directly as the p.d.f. of the tracer particle sojourn times since each time of particle entry and exit is detected. The dimensionless variance was found to vary with operating conditions used from 0.18 to 0.61 indicating modest to severe backmixing of the solids corresponding to 2–6 stirred tanks in series. This is in sharp contrast to the liquid which is essentially in plug flow with dimensionless variance of the liquid RTD always <0.1 (more than 10 CSTRs) [14]. The extent of the solids backmixing tends to increase with liquid velocity and, to a lesser extent, with S=L ratio. In more advanced reactor modeling we would like to have a quantitative description of the local dispersion
that is superimposed on the mean flow field. Typically, one introduces a bundle of particles and from their spreading in time calculates the eddy diffusivity. In our CARPT experiment we use the N independent tracer particle trajectories through the riser compartment of interest to calculate the local eddy diffusivity. Using ergodic hypothesis this ensemble of trajectories is viewed as generated by the swarm of particles released at the same time. The eddy diffusivities are obtained then by the appropriate integrals of the fluctuating Lagrangian velocity auto-correlation coefficients [14,15]. Typically two types of plots are seen [14]. The radial auto-correlation oscillates in a damped manner indicating the presence of the reflective column wall (i.e. the particle disperses back into the column upon reaching it). Such behavior is characteristic of wall bounded flows. The axial autocorrelation decays monotonically as the particles become uncorrelated with the previous state. An example of the calculated local eddy diffusivities as functions of time is shown in Fig. 8. In the radial direction the diffusivities reach a peak value after short times (0.1 s) then decay and level off asymptotically after 0.3 s. The asymptotic values are less than half of peak values. In the axial direction the
754
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
Fig. 7. Azimuthally averaged time-averaged solids velocity (U1 ¼ 23 cm/s).
asymptotic values of eddy diffusivity are reached monotonically after about 0.15 s. These asymptotic values in the axial direction are an order of magnitude larger than the peak values of radial diffusivities. Clearly the asymptotic values of both radial and axial diffusivities are position (radial location) dependent. However, upon averaging the eddy diffusivities over a cross-section we can represent the riser with a single number for axial and radial eddy diffusivity which would likely be done for a model to be used in reactor analysis. Both axial and radial eddy diffusivities increase with the increase in liquid superficial velocity and, at fixed liquid velocity, increase mildly with solid to liquid flux ratio. 4.2. Modeling The information generated by CARPT-CT, and briefly described above, provides us with the previously
Fig. 8. Local solids diffusivities.
unavailable opportunities to model the riser reactor based on realistic flow patterns. For example, the ubiquitously used plug flow model for the liquid and solid can now be replaced at the minimum with the plug flow model for liquid and ADM for the solids where the axial dispersion coefficient is readily calculable from the observed solids RTDs. Or in a two dimensional (2D) convection–diffusion model one could use plug flow of liquid, the observed axial solids velocity profile and the CARPT determined axial and radial eddy diffusivity to describe the flow pattern. The question then arises whether the needed parameters for such reactor models can be obtained via CFD models rather than by experimentation. If they can, this would provide confidence in scale-up. To address this situation we have adopted the Euler– Euler interpenetrating two fluid model to describe the motion of the liquid and solid phase. The customary [17] ensemble averaged equations of continuity and mo-
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
mentum conservation for each phase are shown here in Table 1, and are written assuming negligible interphase exchange of mass. No reactions are considered, setting the right hand side of the continuity equations to zero. The hydrodynamic pressure is assumed shared by both phases. In the case of solids, collisions between individual particles contribute to the normal stress while in a continuum representation of solids this is reflected in the ‘‘solids pressure’’. The stress term in the liquid phase requires closure for the turbulent interaction of liquid eddies. The shear stress term in the solids phase captures via closure the shear resulting from collisions of solid particles, the normal component of which was accounted for in the solids pressure. The momentum exchange (drag term) represents the surface integral of all the forces that act in the interface between phases. Here it is assumed that it can be represented as a momentum exchange term multiplied by the slip velocity between phases. Our estimates [14,15] of the Bagnold number, which is the ratio of inertial to viscous forces, based on
CARPT data, indicates that it is of order 100. Hence, it is necessary to model explicitly both particle–particle and particle–fluid interactions. At Ba > 450 the particle motion is governed entirely by particle–particle interactions and at Ba < 40 by fluid viscosity, so our experiments are in the middle. Since rigorous theory that would account exactly for particle–liquid interactions due to mean velocity, particle interactions with fluctuating liquid velocity, mean velocity interactions with fluctuating ones for the solids and the liquid is not available, we have resorted to a judicious choice of empirical correlations. Solid phase fluctuations were modeled using the granular temperature concept. The governing equations are shown in Table 2. Partial justification for this is provided by the fact that experimentally observed p.d.fs of solids fluctuating velocity in various column compartments can be approximated with a Gaussian distribution [14]. The equation for the granular temperature is derived in the usual way by appropriate ensemble averaging and its right hand side
Table 1 Two fluid approach
755
Phases treated as interpenetrating continua Probability of occurrence of a given phase in multiple realizations of the flow is given by volume fraction of that phase Equations written for ensemble averaged realizations of each phase Both fluids are treated as incompressible o ðek qk Þ þ r ðek qk~ uk Þ ¼ 0 es þ ef ¼ 1; k ¼ f for fluid; k ¼ s for solids ot o ðef qf~ uf Þ þ r ðef qf~ uf ~ uf Þ ¼ ef rp þ r sf þ ef qf~ us ~ uf Þ þ ~ Ff g þ Ksf ð~ ot o ðes qs~ us Þ þ r ðes qs~ us ~ us Þ ¼ es rp rps þ r ss þ es qs~ uf ~ us Þ þ ~ Fs g þ Ksf ð~ ot
Table 2 Granular flow kinetic theory Inspired by kinetic theory of gases, phenomenological model for micro-scale interactions Thermodynamic temperature replaced by granular temperature, which is a flow property Closures for various terms in granular kinetic energy transport equation written in terms of granular temperature h ¼ 23ks
756
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
contains as the first term the production of kinetic energy by interaction of the normal and shear stress matrix with the mean velocity field by which energy is extracted from the mean field into the fluctuating field. The second term is the transport of kinetic energy by conduction, i.e. flow of energy caused by elastic and inelastic collisions between solid particles. The third term represents the dissipation of kinetic energy due to collisions while the fourth is the energy exchange due to interactions with liquid phase fluctuations. The solids pressure, shear viscosity and conductivity are calculated using the expressions of Chapman and Cowling [18] as modified by Ding and Gidaspow [19]. The radial distribution function used by Louge et al. [20] was employed. Collisional terms dominate the contribution for solids pressure and viscosity for the solids volume fraction encountered in our system, while both kinetic and collisional terms are important in determining solids conductivity. Energy dissipation through collisions was modeled as Chapman and Cowling [18] and the restitution coefficient is assumed very close to one (0.98) due to the damping effect of the liquid phase. The variation in this coefficient from 0.95 to 1 does not affect the results. The dissipation due to the friction of the fluctuating solids and liquid is modeled consistent to the drag formulation for the interaction of mean velocities. For liquid–solid drag we have used the correlation by Wen and Yu [21] based on systems with similar solids loadings. Finally to close the equations for liquid phase turbulence a modified k-e model is used that accounts for the presence of the dispersed phase [22]. Appropriate wall boundary conditions, explained in detail in the thesis of Roy [14] were used. Computation was performed in the framework of the Fluent two fluid model version 4. Both 2D axisymmetric and 3D simulations were performed. Best comparison with data was obtained with 3D simulations since the flow clearly has a 3D character. Such 3D simulations using 39,800 cells produce good agreement of calculated time averaged solids axial velocity profile and solids holdup profile with CARPT and CT data as shown in Fig. 9. Averaging is done from 25 to 100 s simulation time. Granular temperature predictions are also good except in the wall region. However the 3D simulation took more than a month on a DEC Alpha 500 MHz machine to execute. To provide faster input to reactor modeling, a 2D simulation coupled with a scalar tracer transport equation was executed and the results for the liquid and solids RTD are in agreement with experimental observation that liquid is close to plug flow and solids are much more backmixed (Fig. 10). The differences in the mean and variance between the computed and experimentally determined values are well within 10%. The following conclusions can be derived from this study: the CARPT-CT technique provides the needed
Fig. 9. Comparison of computed and measured values.
information on solids holdup, velocity and mixing for establishment of more accurate reactor models. While liquid in a liquid–solid riser is close to being in plug flow [14], solids are extensively backmixed. The Euler–Euler model is capable of computing reasonably well the overall features of the reactor scale liquid and solids flow via 2D axisymmetric simulation. The true transient nature of the flow is, however, only revealed by 3D simulations which produce results for time averaged velocity, holdup and kinetic energy in good agreement with experimental observations.
5. Bubble column flows Let us now examine a bubble column, which is the reactor of choice for conversion of natural gas and
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
Fig. 10. Computed RTD’s of the liquid and solids at UL ¼ 20 cm/s, S=L ¼ 0:15.
synthesis gas to clean fuels and is used in a great number of industrial applications. The flow is buoyancy driven as the gas superficial velocity exceeds by orders of magnitude the liquid superficial velocity. For very small catalyst particles the slurry phase can be treated as pseudohomogeneous. The key to proper modeling of the bubble column reactor is to understand and quantify the reactor scale flow pattern. Current models are based on either assumptions of ideal flow (e.g. plug flow for the gas and perfect mixing for the liquid) or on the ADM. However, no truly predictive correlations for the dispersion coefficient can be found. Also ignored is the fact that most columns need to operate in the churn turbulent regime at high gas velocities and high gas holdup. 5.1. Experimental data Using CT we have confirmed that gas holdup profiles are in the time averaged sense symmetric with more gas in the middle than by the walls while in bubbly flow the profiles are relatively flat. In churn turbulent flow the gas holdup profile becomes parabolic [7]. Using CARPT with a neutrally buoyant particle, we have confirmed the findings of Hills [23] and many others that there is a distinct time averaged liquid recirculation profile with liquid rising in the middle and falling by the walls. The time averaged liquid flow pattern is remarkably symmetric and forms a single recirculation cell [11]. However, in addition to the ensemble averaged velocity
757
profile, CARPT also provides eddy diffusivities in the radial and axial direction. The values in the axial direction are an order of magnitude larger than the values in radial direction [11]. Diffusivities increase with superficial gas velocity and with column diameter. Axial diffusivities are the highest in the region of downward time averaged liquid flow which can be explained by the tendency of bubbles to resist downward motion creating more backmixing of the liquid. As already mentioned, our CT and CARPT data reveal conclusively that a time averaged liquid recirculation profile exists through most of the column and is driven by a pronounced radial gas holdup profile [7,10,11]. The fact that even in churn turbulent flow the time averaged flow pattern is symmetric is clearly established. The liquid rises up in the middle and falls by the walls. The long term pattern is axisymmetric. However, at each point (cell, compartment) one can establish, based on CARPT data, a histogram of instantaneous velocities which yield zero means for radial and azimuthal velocities through most of the column [11]. Axial Lagrangian auto-correlation coefficient and axial–radial cross-correlation coefficient can also be established and are indispensable in evaluation of eddy diffusivities. Turbulent stresses for the 44 cm (1800 ) diameter column in churn turbulent flow at superficial gas velocity of 10 cm/s exhibit the usual trends for bubble columns [11,24]. The axial normal stresses are twice as high as the normal radial or azimuthal stresses. The only non-zero shear stress is the r–z one which peaks close to the velocity inversion point. It is important to note that turbulent stresses measured by CARPT are in good agreement with the values obtained by hot film anemometry and laser doppler anemometry under conditions of bubbly flow when the latter two techniques can be made to work [25,26]. Application of the R=S analysis to the CARPT data provides the Hurst exponents which indicate diffusion type mixing in the radial direction and strong convective effects in the axial direction [27]. Application of the Lagrangian and qualitative dynamics tools analysis to CARPT data clearly reveals the role of bubble wakes on the motion of the tracer particle and identifies and quantifies the chaotic fluid dynamics of the liquid phase. The effect of column diameter and operating conditions on the Kolmogorov entropy, strange attractors and kinematics of bubble wakes can then be established [28]. We have shown that the experimental information obtained by CARPT-CT can be utilized in an improved 2D phenomenological reactor model for the column presented schematically in Fig. 11 [29,30]. Based on rigorous arguments one can show that the time averaged velocity profile needed for the 2D scalar transport model is the ensemble averaged velocity determined from CARPT and that the needed liquid holdup profile is indeed the time averaged holdup determined by CT. The
758
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
Fig. 11. Suggested model for tracer spreading in the liquid phase of the bubble column.
assumption made is that locally the unknown crosscorrelation of the liquid velocity––tracer concentration fluctuations can be closed by the Bousinesq approximation i.e. by a product of the diffusivity tensor and the gradient of the tracer concentration. Moreover, it is assumed that diffusivities are the ones determined by CARPT. The usefulness of these assumptions is tested by feeding the CARPT measured velocity profile and CT measured holdup profile into the 2D convection diffusion model together with the CARPT measured diffusivities. One is now able to calculate the normalized exit mixing cup concentration for the impulse injection of tracer at the inlet and obtain the E-curve (exit age density function p.d.f. of the RTD). Comparison with the experimentally determined E-curve in response to the dye impulse injection is astonishingly good as shown in Fig. 12. This means that the model is capable of predicting reactor performance exactly for all linear kinetics. The success of this approach in being able to produce tracer response in agreement with measurement on a real hot pilot plant bubble column 1800 in diameter has been discussed elsewhere [29,30]. 5.2. Modeling We have qualitatively compared CARPT data (outline of particle trajectory and mean distance traveled) with the Euler–Lagrange simulation of bubbly flows [31]. More than 50,000 bubbles (assuming no bubble– bubble interactions) were simulated in a 3D simulation. Reasonable agreement between data and simulations was observed.
Since churn turbulent flow is of most current interest, Lagrange–Euler simulations are not practical, and the Euler–Euler interpenetrating two fluid model is a natural choice. The issues remain as to how to represent phase–phase interactions such as lift, added mass and drag and how to describe bubble induced turbulence. A simple approach that neglects lift, uses classical approach for added mass, models the drag force based on a sphere affected by presence of other spheres has been tried successfully in the Euler–Euler interpenetrating two fluid model to compute liquid velocities, gas holdup profiles, turbulent parameter and even diffusivities in bubble columns. The models were implemented and executed either in the CFDLIB code framework of Los Alamos or in Fluent 4.5 codes. Usually available closures, based on steady state expressions were used for drag and added mass. Liquid phase turbulence was accounted for either via the bubble induced turbulence [32,33] or via k–e model adjusted for the presence of the second phase [34]. Until recently an average bubble size was assumed for evaluation of the interfacial forces. At present we have successfully implemented the bubble population balance model with breakup coalescence models to predict the evolution of the bubble size distributions. Reasonable agreement between computed and experimentally observed results were reported for bubbly flow in quasi-2D columns [32,33] and for both bubbly and churn turbulent flow in 3D columns [34,35]. Pan et al. [33] clearly demonstrated that not only the time averaged liquid velocity profile was correctly calculated by the CFDLIB Euler–Euler code but that the right magnitude was predicted for liquid turbulent kinetic energies and stress profile in bubbly flows. Moreover, the calculated periods of enriched gas plume fluctuation agreed well with the experimentally observed ones. This provided confidence that the Euler–Euler two fluid model is capable of calculating liquid velocities in bubbly flow within engineering accuracy. In churn turbulent flow, 2D axisymmetric calculations produced liquid time average velocity in reasonable agreement with the experimental CARPT data [34,35]. Correct overall gas holdup magnitudes were calculated but computed gas holdup profiles did not exhibit the CT experimentally observed parabolic radial profile. 3D computations, which take a long time to execute, came closer to producing the observed holdups but substantial discrepancies still exist. Work is in progress to determine whether the incorporation of the population balance model and coalescence break up models can alter sufficiently the computed holdup profile. Most importantly, from the reaction engineering viewpoint, it has recently been shown [35] that the CFDLIB Euler–Euler model with the above discussed closures can, via 3D computations, produce the eddy diffusivities whose radial profile is in
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
759
Fig. 12. Prediction of non-volatile tracer impulse response.
close agreement with values obtained from CARPT experiments. Hence, CFD offers possibilities in calculating the parameters needed for reactor modeling, such as time averaged velocity profile and eddy diffusivities. As mentioned earlier a 2D convection–diffusion model that utilizes CFD calculated flow parameters can be quite useful in predicting tracer responses and reactor performance of industrial size reactors.
and in describing the liquid flow field in churn turbulent bubble columns is established. The capabilities of the Euler–Euler interpenetrating two fluid model in generating flow pattern information in agreement with CARPT-CT data that can form a basis for reactor modeling are also shown and are significant for industrial applications.
7. Concluding remarks 6. Practical significance/usefulness This paper describes the key features of and provides the key references for CARPT-CT. The techniques are unique in their capabilities of providing the flow field description in large enclosures for opaque systems. Their utility in obtaining the solids holdup distribution, velocity and kinetic energy profiles in a liquid–solid riser
The CARPT-CT techniques provide us with the unique opportunities to assess holdup and velocity fields in opaque two phase systems. The resolutions of both techniques allow reliable information to be generated for large scale motion at lower frequencies. While CT provides the time averaged density distribution, CARPT, in addition to time averaged velocities, allows estimation
760
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
of a number of turbulent quantities, e.g. kinetic energy, stresses, eddy diffusivities, etc. The wealth of information provided by CARPT-CT allows us to quantify the flow fields in systems like liquid–solid risers and bubble columns which leads to improved and more detailed reactor models. CFD computations, based on the Euler–Euler interpenetrating two fluid model have been shown to be capable of producing information in good agreement with data, which renders hope for CFD based scale-up in the future. It is fair to conclude that CARPT-CT are indeed a valuable measurement tool for opaque multiphase systems and that the Euler–Euler model has potential in simulation of reactor scale multiphase flows at conditions of industrial interest.
Acknowledgements At the end I would like to thank the many collaborators that have produced all the results presented here. Above all, Narsi Devanathan, Yubo Yang and Sailesh Kumar deserve credit for designing, constructing and implementing the hardware and software for CARPTCT in CREL. Shantanu Roy and Abdenour Kemoun executed the experiments in the liquid–solid riser discussed here, and Shantanu modeled them using Fluent with the help of Jay Sanyal (Fluent). Sue Degaleesan, Jinwen Chen, Booncheng Ong experimented with the bubble column, while Yu Pan and S. Roy solved CFD models with CFDLIB and Fluent. Muthanna Al-Dahhan offered helpful suggestions and Bernard Toseland of Air Products was instrumental in focussing us to meet industrial needs. Finally, we are thankful to our financial sponsors: Exxon-Mobil for donating the CARPT-CT elements, Department of Energy and all CREL industrial sponsor companies from five continents (ABB Lummus, Air Products, Bayer, Chevron, Conoco, Dow Chemicals, DuPont, Elf Atofina, ENI Technologies, Exxon-Mobil, IFP, Intevep, Mitsubishi, Praxair, Sasol, Shell, Solutia, Statoil, Synetix-ICI, UOP). Their financial support made this work possible.
References [1] P.L. Mills, Recent development and advances in laboratory-scale multiphase reactors, Plenary Lecture, First North American Symposium on Chemical Reaction Engineering (NASCRE1), Houston, TX, 2001. [2] G.F. Froment, K.B. Bischoff, Chemical Reactor Analysis and Design, second ed., Wiley, NY, 1990. [3] J. Chaouki, F. Larachi, M.P. Dudukovic, Non-invasive tomographic and velocimetric monitoring of multiphase flows, Ind. Eng. Chem. Res. 36 (11) (1997) 4476–4503.
[4] J. Chaouki, F. Larachi, M.P. Dudukovic (Eds.), Non-invasive monitoring of multiphase flows, Elsevier, Amsterdam, 1997. [5] S.B. Kumar, D. Moslemian, M.P. Dudukovic, A c-ray tomographic scanner for imaging of void distribution in two-phase flow systems, Flow Meas. Instrum. 6 (3) (1995) 61–73. [6] S.B. Kumar, M.P. Dudukovic, Computer assisted gamma and x-ray tomography: Application to multiphase flow systems, in: Chaouki et al. (Eds.), Non-Invasive Monitoring of Multiphase Flows, Elsevier, Amsterdam, 1997, pp. 47–103, Chapter 2. [7] S.B. Kumar, D. Moslemian, M.P. Dudukovic, Gas holdup measurements in bubble columns using computed tomography, AIChE J. 43 (6) (1997) 1414–1425. [8] J.S. Lin, M.M. Chen, B.T. Chao, A novel radioactive particle tracking facility for measurement of solids motion in gas fluidized beds, AIChE J. 31 (3) (1985) 465–473. [9] N. Devanathan, Investigation of liquid hydrodynamics in bubble columns via computer-automated radioactive particle tracking (CARPT), D.Sc. Thesis, Washington University, St. Louis, MO, 1991. [10] N. Devanathan, D. Moslemian, M.P. Dudukovic, Flow mapping in bubble columns using CARPT, Chem. Eng. Sci. 45 (8) (1990) 2285–2291. [11] S. Degaleesan, Fluid dynamic measurements and modeling of liquid mixing in bubble columns, D.Sc. Thesis, Washington University, St. Louis, MO, 1997. [12] F. Larachi, G. Kennedy, J. Chaouki, A c-ray detection system for 3-D particle tracking in multiphase reactors, Nucl. Instrum. Meth. A 338 (2-3) (1994) 568–576. [13] P. Gupta, Bubble column dynamics: experiments and models, D.Sc. Thesis, Washington University, St. Louis, MO, May, 2002. [14] S. Roy, Quantification of two phase flows in liquid–solid risers, D.Sc. Thesis, Washington University, St. Louis, MO, December, 2000. [15] S. Roy, M.P. Dudukovic, Flow mapping and modeling of liquidsolid risers, Ind. Eng. Chem. Res. 40 (23) (2001) 5440–5454. [16] S. Roy, A. Kemoun, M. Al-Dahhan, M.P. Dudukovic, A method for estimating the solids circulation rate in a closed-loop circulating fluidized bed, Powder Technol. 121 (2–3) (2001) 213–222. [17] D.A. Drew, Mathematical modeling of two phase flow, Ann. Rev. Fluid Mech. 15 (1983) 261–291. [18] S. Chapman, T.G. Cowling, Mathematical Theory and Nonuniform Gases, third ed., Cambridge University Press, Cambridge UK, 1990. [19] J. Ding, D. Gidaspow, A bubbling fluidization model using kinetic theory of granular flow, AIChE J. 36 (4) (1990) 523–538. [20] M.Y. Louge, E. Mastorakos, J.T. Jenkings, The role of particle collisions on pneumatic transport, J. Fluid Mech. 231 (1991) 345– 359. [21] C.Y. Wen, Y.H. Yu, Mechanics of fluidization, Chem. Eng. Progr. Symp. Sci. 62 (62) (1966) 100–111. [22] S.E. Elghobashi, T.W. Abou-Arab, A two-equation turbulence model for two-phase flows, Phys. Fluids 26 (4) (1983) 931–938. [23] J.H. Hills, Packed nonuniformity of velocity and voidage in a bubble column, Trans. Inst. Chem. Eng. 52 (1974) 1–9. [24] J. Chen, L. Fan, S. Degaleesan, P. Gupta, M.H. Al-Dahhan, M.P. Dudukovic, B.A. Toseland, Fluid dynamic parameters in bubble columns with internals, Chem. Eng. Sci. 54 (13/14) (1999) 2187– 2197. [25] T. Menzel, T. Inderweide, O. Standacher, O. Wein, U. Onken, Reynolds shear stress for modeling of bubble column reactors, Ind. Eng. Chem. Res. 29 (6) (1990) 988–994. [26] R.F. Mudde, J.S. Groen, H.E.A. van den Akker, Liquid velocity field in bubble columns: LDA experiments, Chem. Eng. Sci. 52 (21/22) (1997) 4217–4224. [27] Y.B. Yang, N. Devanathan, M.P. Dudukovic, Liquid backmixing in bubble columns via computer-automated radioactive particle tracking (CARPT), Exp. Fluids 16 (1) (1993) 1–9.
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761 [28] M. Cassanello, F. Larachi, A. Kemoun, M. Al-Dahhan, M.P. Dudukovic, Inferring liquid chaotic dynamics in bubble columns using CARPT, Chem. Eng. Sci. 56 (21/22) (2002) 6125–6134. [29] S. Degaleesan, S. Roy, S.B. Kumar, M.P. Dudukovic, Liquid mixing based on convection and turbulent dispersion in bubble columns, Chem. Eng. Sci. 51 (11) (1996) 2721–2725. [30] S. Degaleesan, M.P. Dudukovic, B.A. Toseland, B.L. Bhatt, A two compartment convective diffusion model for slurry bubble column reactors, Ind. Eng. Chem. Res. 36 (11) (1997) 4670–4680. [31] N. Devanathan, M.P. Dudukovic, A. Lapin, A. Lubbert, Chaotic flow in bubble column reactors, Chem. Eng. Sci. 50 (16) (1995) 2661–2667.
761
[32] Y. Pan, M.P. Dudukovic, M. Chang, Dynamic simulation of bubbly flow in bubble columns, Chem. Eng. Sci. 54 (13/14) (1999) 2481–2489. [33] Y. Pan, M.P. Dudukovic, M. Chang, Numerical investigation of gas driven flow in 2-D bubble columns, AIChE J. 46 (3) (2000) 434–449. [34] J. Sanyal, S. Vasquez, S. Roy, M.P. Dudukovic, Numerical simulation of gas–liquid dynamics in cylindrical bubble column reactors, Chem. Eng. Sci. 55 (21) (1999) 5071–5083. [35] Y. Pan, M.P. Dudukovic, CFD simulations of bubble columns, 6th World Congress of Chemical Engineering, Melbourne, Australia, September, 2001.
Opaque multiphase flows: experiments and modeling M.P. Dudukovic
*
The Chemical Reaction Engineering Laboratory (CREL), Washington University, One Brookings Drive, Campus Box 1198, St. Louis, MO 63130-4899, USA Accepted 28 November 2001
Abstract This brief review deals with the characterization of multiphase flows needed in reaction engineering and describes the techniques and methods available to accomplish it. It emphasizes the fact that most multiphase systems of industrial interest (e.g. stirred tanks, slurry bubble columns, gas–solid or liquid–solid risers and fluidized beds, pneumatic transport reactors, etc.) operate in churn turbulent flow at large volume fraction of the dispersed phase and are opaque. This inevitably leads to the use of the Euler–Euler, interpenetrating fluid model in describing their fluid dynamics, which raises issues of closure for the interfacial momentum terms and turbulence. Hence, computational fluid dynamic (CFD) models based on the Euler–Euler approach with selected closures require experimental validation. Due to the non-transparency of the systems involved, non-conventional experimental techniques are required. The paper briefly describes two such techniques which have been successfully implemented at the Chemical Reaction Engineering Laboratory for characterization of opaque, multiphase flows. One is the gamma ray computed tomography (CT), the other is the computer automated radioactive particle tracking (CARPT). Their use in characterizing diverse systems such as liquid–solid risers and bubble columns is illustrated as well as the utility of CARPT-CT data in validation of the Euler–Euler CFD models. The goal is to enable the reader to grasp the advantages and limitations of the described experimental techniques and computational fluid dynamic models. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Riser; Bubble column; Gas–liquid; Liquid–solid; Radioactive particle tracking; Computed tomography
1. Introduction Applications of multiphase technology are prevalent in petroleum refining, synthesis gas conversion to fuels and chemicals, bulk commodity chemicals, specialty chemicals, conversion of undesired byproducts to recyclable products, manufacture of polymers, biotechnology and pollution abatement. It is estimated that the value of shipments generated by industry that predominantly uses multiphase technology exceeds 640 billion dollars a year in the United States alone [1]. As a rule the multiphase reactor is the heart of any process in which chemical transformations take place, because reactor performance determines the number and size of needed separation units in front and after the reactor and, hence *
Corresponding author. Tel.: +1-314-935-6021; fax: +1-314-9354832. E-mail address: [email protected] (M.P. Dudukovic).
dictates the economics of the whole process. It is thus essential to quantify various measures of reactor performance, such as conversion, selectivity, production rate, etc., as a function of the design and operating parameters which requires appropriate descriptions of micro- and macro-transport. The discipline of reaction engineering attempts precisely to do that by quantifying the transport-kinetic interactions on various scales. Starting with species continuity equations and the energy balance (Fig. 1), one recognizes the need to describe the phenomena that occur in the reactor on a myriad of scales. Molecular scale knowledge allows us to write the kinetic rate expressions, R, based on first principles. Understanding the interaction of flow and transport in a single eddy or on a single catalyst particle scale enables us to modify the kinetic forms by local transport effects, g, and describe more precisely the source/sink terms in the species and energy balance equations as well as to describe transport between phases. The left hand side
0894-1777/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 8 9 4 - 1 7 7 7 ( 0 2 ) 0 0 1 8 5 - 1
748
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
Nomenclature C Drr Dzz E F g k khs ks Ksf p Q r R T t u, U z
concentration of scalar quantity (like tracer), kg/m3 radial eddy diffusivity, m2 /s axial eddy diffusivity, m2 /s E-curve, exit age density function, s1 force, N gravitational acceleration, m/s2 mass exchange coefficient, m/s coefficient of granulary conductivity, kg/m s kinetic energy per unit mass of solids, m2 /s2 solids–fluid exchange coefficient, kg/m3 s pressure, Pa volumetric flow rate, m3 /s radial position, m reaction rate, mol/m3 s temperature, K time, s velocity, m/s axial distance, m
Greek symbols DH heat of reaction, J/mol
chs e g h q s /fs
collisional dissipation term, kg/m s3 void fraction, dimensionless effectiveness factor, dimensionless granular temperature, m2 /s2 density, kg/m3 shear stress, N/m2 solids energy dissipation by correlations between solids and liquid velocity fluctuations, kg/m s3
Subscripts b bulk fluid f fluid phase g gas phase l liquid phase s solid phase k kth phase 1, 2 phase 1, 2 Other Symbols h i average 0 fluctuating quantity
Fig. 1. Reaction engineering models.
of the balance equations, containing the input–output terms due to convection, diffusion and interphase exchange, rests on the description of the flow field and its approximation. Rarely is in reactor engineering the momentum balance solved coupled with the species and energy balance to provide for the exact flow field and for the exact species spatial and temporal distribution. In recent years there has been an increased drive towards the use of fundamentals in reaction engineering, pushing the industrial practice from strictly empirical to increasingly fundamental, mechanism based, descriptions of kinetics. At the same time direct numerical simulations (DNS) allow at certain conditions an increasingly accurate description of flow fields and transport around small groups of particles helping steadily
advance the understanding on the single eddy/particle scale. In contrast, the description of reactor scale flows and mixing has not progressed much and in the industrial applications it is at best at the level of ideal reactor assumptions of plug flow or perfect mixing. In a plug flow reactor (PFR) a piston type one dimensional (1D) flow with flat velocity profile is assumed. In a perfectly mixed continuous flow stirred tank reactor (CSTR) it is assumed that the feed looses identity instantly and that all of the reactor content is at a uniform composition and temperature equal to those of the outlet stream. The problem of mixing and flow is avoided by assuming that the time constants for those phenomena are orders of magnitude shorter than the characteristic reaction time. This, however, is hardly ever the case! Typical models
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
749
Fig. 2. Flow pattern models for reactors with two moving phases.
for a two phase reactor with two flowing phases are sketched in Fig. 2. Either plug flow or perfect mixing is assumed in each phase. Plug flow and perfect mixing combination (not shown) is also used [2]. On occasion an empirical 1D axial dispersion model (ADM) is used in each phase to fill the gap between PFR and CSTR. The dispersion coefficient seldom can be scaled up. In each of these models the interphase exchange is lumped into an overall exchange coefficient that needs to be multiplied by an average driving force to describe the transport. A much more desirable approach would be to rely on phenomenological models that capture the essence of the flow field. Such phenomenological models can be developed based on experimental observations coupled with CFD model calculations. Ultimately, CFD calculations should be used to assess the flow field of each phase. Unfortunately, for multiphase flow, CFD calculations are still subject to uncertainty, especially since DNS cannot be used for simulation of large reactors. It seems that for simulation of pilot plant and commercial scale multiphase reactors the Euler–Euler k-fluid model is the best suited. This, of course, requires various closures for the interphase interaction forces and intraphase turbulence. Thus, the results produced by the multiphase CFD codes need to be experimentally validated. Since the systems to be interrogated are opaque, and we cannot ‘‘see’’ into them, most of the non-invasive techniques that are based on optical methods or use lasers cannot be employed. Fortunately, as recent reviews illustrate [3,4], there are techniques which can provide us with the desired information. In this paper the focus is on introducing two of them: gamma ray computer tomography (computed tomography, CT), for measurement of holdup profiles, and computer aided radioactive particle tracking (computer automated radioactive particle tracking, CARPT), for measurement
of velocity profiles and mixing parameters. We show how these techniques can be used to obtain information in systems of industrial interest such as a liquid–solid riser, or a gas–liquid bubble column. The ability of the available CFD Euler–Euler models to correctly calculate the observed hydrodynamic quantities is also briefly discussed.
2. Gamma ray computed tomography (CT) To measure the phase holdup distribution at any desired cross-section of the column we implemented a gamma ray source based fan beam type CT unit (Fig. 3). A collimated hard source (100 mCi of Cs-137) is positioned opposite to eleven 200 sodium iodide detectors in a fan beam arrangement. The lead collimators in front of the detectors have manufactured slits and the lead assembly can move so as to allow repeated use of the same detectors for additional projections. A 360° scan can be executed at any desired axial location and columns from 100 to 1800 in diameter can be scanned in the current configuration. The principle of computed tomography (CT) is simple. From the measured attenuation of the beams of radiation through the two phase mixture (projections) we calculate due to the different attenuation by each phase, the distribution of phases in the cross-section that was scanned. In our Chemical Reaction Engineering Laboratory (CREL) CT unit we achieve 3465–4000 projections and in the original unit we had a spatial resolution of 3 mm and density resolution of 0.02 g/cm2 . Of course we need a long time to scan the cross-section (about 45 min) and, hence, can only determine the time averaged density distribution. Details of the CT unit at CREL are described elsewhere [5–7]. Among the suggested techniques for reconstruction, e.g. convolution or filtered back projection, algebraic
750
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
Fig. 3. CT setup (a) schematic (b) unit in CREL.
reconstruction and estimation–maximization algorithms (E–M), we have found the E–M algorithm to be the best [6]. The CT hardware and software used in CREL and described above is quite suitable to provide us with reliable images of time averaged cross-sectional density distribution of the opaque system of interest at any desired elevation. It accomplishes that in a non-invasive manner and provides us with the ability to assess the spatial holdup distribution in two phase systems with good resolution compared to the vessel scale. It cannot, at present, yield ‘‘instantaneous images’’ and capture the dynamics of the system, nor can it provide very fine
resolution of the order of a single catalyst particle or a single small bubble. Since it is the reactor scale flow description that we are interested in, our CT unit provides valuable information for that purpose as illustrated above.
3. Computer automated radioactive particle tracking (CARPT) We have adopted CARPT (Fig. 4), the technique based on following the motion of a single tracer particle, as a method of mapping the velocity field in the whole
Fig. 4. CARPT setup (a) schematic (b) unit in CREL.
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
system. This Lagrangian technique was introduced for monitoring the motion of solids in fluidized beds by Lin et al. [8] and was adapted for measurement of liquid velocities in bubble columns by Devanathan and coworkers [9,10]. A single radioactive particle dynamically similar to the phase to be traced is introduced into the system. When solids motion is monitored in ebullated beds, fluidized beds or slurries, a particle of the same size and density as the solids in the system is prepared. For monitoring of liquid motion a neutrally buoyant particle is made [9]. For example, a polypropylene bead of about 2 mm is hollowed out, a small amount of Scandium 46 is inserted and a polypropylene plug is put in place, so that the density of the composite particle consisting of polypropylene, Scandium and air gap equals that of the liquid. Application of a thin film metallic coating assures that bubbles do not preferentially adhere to the particle. An array of scintillation detectors is located round the column. In our case up to 32 NaI 200 detectors are used. The detectors are calibrated in situ with the tracer particle to be used to get the counts–distance maps. Conventionally, CARPT calibration was done by positioning the tracer particle (usually Sc-46 of 250 lCi strength was used) at about 1000 known locations and recording the counts obtained at each detector. Various ways of doing this have been described [9,11]. An alternative way of constructing distance–count maps is via modeling of particle emission of photons and transmission and subsequent detection at the detectors. The Monte Carlo method [12,13] in which the photon histories are tracked in their flight from the source, through the attenuating medium and their final detection (or lack of it) at the detector can be used for this purpose. Thus, both the geometry and radiation effects are accounted for in estimation of the detector efficiencies in capturing and recording the photons. This involves calculation of complex three dimensional (3D) integrals which are evaluated using the Monte Carlo approach by sampling modeled photon histories over the main directions of their flight from the source. Once calibration is complete, the tracer particle is let loose in the system and the operating conditions are controlled and kept constant for many hours while the particle is tracked. A regression, least-squares method is used to evaluate the position of the particle. Sampling frequency is adjusted between 50 Hz and 500 Hz to assure the required accuracy. Typically, it is selected at 50 Hz for bubble columns since the finite size particle used as tracer cannot capture motion of frequencies above 25 Hz [11]. A wavelet based filtering algorithm [11] is employed to remove the noise in position readings created by the statistical nature of gamma radiation. The processing of data proceeds along the following lines. Filtered particle positions at each sampling time establish the tracer particle trajectory throughout the system. In systems like bubble columns with batch liquid the par-
751
ticle trajectory is confined to the liquid phase and the tracer particle stays in the vessel at all times. In risers, or bubble columns with flowing liquid, the tracer particle leaves and re-enters the system repeatedly. In either case from filtered particle positions at subsequent sampling times the instantaneous particle velocity is calculated and assigned to the column compartment (a grid is preestablished for the column) into which the mid-point falls. Thus, the instantaneous velocities are obtained repeatedly for each compartment. The run continues until adequate statistics are obtained. Then for each compartment (cell) ensemble average velocities are evaluated, which by ergodic hypothesis are assumed to be the time averaged velocities for each such location in the column. The differences between instantaneous and average velocities allow the estimation of the fluctuating velocities and all their auto- and cross-correlations leading to estimates of turbulent kinetic energy and Reynolds stresses. Most importantly, from the reaction engineering viewpoint, the appropriate Lagrangian correlations lead to the evaluation of the components of the eddy diffusivity tensor. The CARPT technique is very well suited for describing the reactor scale flows in opaque system by noninvasive means. The spatial resolution of 2–5 mm is quite adequate and frequencies up to 25 Hz can be captured with relatively large neutrally buoyant particle while higher frequencies are accessible with smaller particles. The technique cannot provide at present a fine spatial or temporal resolution. Its utility is illustrated below for liquid–solid and gas–liquid flows.
4. Liquid–solid riser flows Liquid–solid riser is one of the reactors of choice for alkylation reactions needed in production of higher octane fuels and linear alkylbenzene. This reactor with novel solid acid catalyst replaces the environmentally problematic processes that used HF and H2 SO4 as liquid catalysts. The solid acid catalyst, however, does deactivate and needs to be replaced. For that reason it is carried by liquid reactants through a riser, in which reaction takes place, and then is recycled through a regeneration unit. For proper reactor design and performance prediction it is necessary not only to know the solids holdup profile but also the solids mean velocity profile and measures of solids backmixing. Before our study the prevailing working assumption was that solids and liquid are in plug flow (piston flow). The average slip velocity between them defined, via Richardson– Zake correlation and the global drift flux model, the total solids holdup in the bed. This was insufficient information for proper analysis of reactor performance. We have provided now a more detailed description of the solids holdup and flow pattern [14,15].
752
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
4.1. Experimental results Our setup consisted of a 6 in. (15 cm) diameter 7 ft (213 cm) long plexiglass column. 2.5 mm glass beads of 2.5 density of water were used in water as the solid– liquid system. The solids mass flux was controlled by water flow through an eductor which was calibrated using a novel technique [16]. Tracer studies in the liquid phase (using KCl as tracer and conductivity probes) confirmed that liquid is in turbulent flow and close to plug flow [14]. The dimensionless variance of the measured liquid residence time distribution (RTD) indicates always more than 10 CSTRs in series at several superficial liquid velocities and three solids to liquid ratios. If liquid were flowing alone, the Reynolds number for the riser would have been between 23,000 and 35,000, so that the liquid flow is definitely turbulent. Solids holdup cross-sectional profiles were measured by gamma ray computer tomography using the CT equipment by Roy [14]. Relatively symmetric cross-sectional distributions are observed which allows us to azimuthally average the solids holdup. In Fig. 5 the radial profiles at three elevations and at three solids to liquid flux (S=L) ratios of 0.10, 0.15 and 0.20 at a liquid superficial velocity of 20 cm/s are presented. Clearly, the solids holdup profiles are quite similar at all elevations and the flow can be considered from the engineering standpoint as fully developed. At all conditions there are slightly more solids at the wall than in the middle of the riser but the radial profiles are not nearly as pronounced as in a gas–solid riser. The magnitude of measurement errors are indicated by bars in Fig. 5. Solids velocity was monitored by a radioactive tracer particle made of an aluminum ball of 2.5 mm in diameter with scandium inserted into it so as to match the density of the glass beads. The settling velocity of the tracer in water was determined to be 30:1 0:1 (cm/s) while the glass beads had a settling velocity of 30:2 0:1 cm/s. Sc-46 was activated to 290 lCi. The tracer particle behaves as the solids being traced. The sampling frequency was 50 Hz and the time of the tracer particle entry and exit of the riser was also carefully recorded providing the distribution of particle residence times and a probability density function (p.d.f.) of independent trajectories. The Lagrangian trace of a single particle trajectory during its sojourn through the riser is shown in Fig. 6 projected onto a horizontal plane and two vertical planes. The time series for the radial and axial particle position are also shown. Several thousands of such trajectories were sampled at each steady operating condition. Clearly the solid particle does not experience plug flow. From these trajectories the instantaneous and ultimately average tracer velocity is calculated at each location of the riser. The vector plots of the ensemble averaged solids velocity vectors projected on various
Fig. 5. Circumferentially averaged time averaged solids holdup distribution in the liquid–solid riser at U1 ¼ 15 cm/s. (a) S=L ¼ 0:10, (b) S=L ¼ 0:15, (c) S=L ¼ 0:20 (bars show standard deviation of azimuthal variation).
vertical planes, not shown here but available in Roy [14] and Roy and Dudukovic [15] reveal a fully developed flow pattern with solids rising in the middle and descending by the wall. Azimuthal averaging of these vector plots produces almost zero average solids radial velocities and radial profiles of the axial solids velocity that are almost independent of elevation from 50 to 150 cm at all solid to liquid (S=L) flux ratios studied as seen in Fig. 7. It is evident that the axial solids velocity profile in the radial direction is very pronounced with the solids in the time average sense rising in the middle of the riser and falling by the walls. Coupled with the
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
753
Fig. 6. Lagrangian trace of tracer particle during one sojourn in the riser (U1 ¼ 20 cm/s; S=L ¼ 0:15) and time series of radial and axial position.
solids holdup profile this provides the information on the extent of solids backflow. The combination of the solids holdup profile and axial velocity profile at all conditions satisfies the overall solids continuity usually within 5% with maximum discrepancy bounded by 10%. For reactor modeling we need to describe the extent of solids backmixing, i.e. departure from the assumed plug flow. The classical way of determining the extent of backmixing is to examine the RTD and its variance. In CARPT experiments the solids RTD is obtained directly as the p.d.f. of the tracer particle sojourn times since each time of particle entry and exit is detected. The dimensionless variance was found to vary with operating conditions used from 0.18 to 0.61 indicating modest to severe backmixing of the solids corresponding to 2–6 stirred tanks in series. This is in sharp contrast to the liquid which is essentially in plug flow with dimensionless variance of the liquid RTD always <0.1 (more than 10 CSTRs) [14]. The extent of the solids backmixing tends to increase with liquid velocity and, to a lesser extent, with S=L ratio. In more advanced reactor modeling we would like to have a quantitative description of the local dispersion
that is superimposed on the mean flow field. Typically, one introduces a bundle of particles and from their spreading in time calculates the eddy diffusivity. In our CARPT experiment we use the N independent tracer particle trajectories through the riser compartment of interest to calculate the local eddy diffusivity. Using ergodic hypothesis this ensemble of trajectories is viewed as generated by the swarm of particles released at the same time. The eddy diffusivities are obtained then by the appropriate integrals of the fluctuating Lagrangian velocity auto-correlation coefficients [14,15]. Typically two types of plots are seen [14]. The radial auto-correlation oscillates in a damped manner indicating the presence of the reflective column wall (i.e. the particle disperses back into the column upon reaching it). Such behavior is characteristic of wall bounded flows. The axial autocorrelation decays monotonically as the particles become uncorrelated with the previous state. An example of the calculated local eddy diffusivities as functions of time is shown in Fig. 8. In the radial direction the diffusivities reach a peak value after short times (0.1 s) then decay and level off asymptotically after 0.3 s. The asymptotic values are less than half of peak values. In the axial direction the
754
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
Fig. 7. Azimuthally averaged time-averaged solids velocity (U1 ¼ 23 cm/s).
asymptotic values of eddy diffusivity are reached monotonically after about 0.15 s. These asymptotic values in the axial direction are an order of magnitude larger than the peak values of radial diffusivities. Clearly the asymptotic values of both radial and axial diffusivities are position (radial location) dependent. However, upon averaging the eddy diffusivities over a cross-section we can represent the riser with a single number for axial and radial eddy diffusivity which would likely be done for a model to be used in reactor analysis. Both axial and radial eddy diffusivities increase with the increase in liquid superficial velocity and, at fixed liquid velocity, increase mildly with solid to liquid flux ratio. 4.2. Modeling The information generated by CARPT-CT, and briefly described above, provides us with the previously
Fig. 8. Local solids diffusivities.
unavailable opportunities to model the riser reactor based on realistic flow patterns. For example, the ubiquitously used plug flow model for the liquid and solid can now be replaced at the minimum with the plug flow model for liquid and ADM for the solids where the axial dispersion coefficient is readily calculable from the observed solids RTDs. Or in a two dimensional (2D) convection–diffusion model one could use plug flow of liquid, the observed axial solids velocity profile and the CARPT determined axial and radial eddy diffusivity to describe the flow pattern. The question then arises whether the needed parameters for such reactor models can be obtained via CFD models rather than by experimentation. If they can, this would provide confidence in scale-up. To address this situation we have adopted the Euler– Euler interpenetrating two fluid model to describe the motion of the liquid and solid phase. The customary [17] ensemble averaged equations of continuity and mo-
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
mentum conservation for each phase are shown here in Table 1, and are written assuming negligible interphase exchange of mass. No reactions are considered, setting the right hand side of the continuity equations to zero. The hydrodynamic pressure is assumed shared by both phases. In the case of solids, collisions between individual particles contribute to the normal stress while in a continuum representation of solids this is reflected in the ‘‘solids pressure’’. The stress term in the liquid phase requires closure for the turbulent interaction of liquid eddies. The shear stress term in the solids phase captures via closure the shear resulting from collisions of solid particles, the normal component of which was accounted for in the solids pressure. The momentum exchange (drag term) represents the surface integral of all the forces that act in the interface between phases. Here it is assumed that it can be represented as a momentum exchange term multiplied by the slip velocity between phases. Our estimates [14,15] of the Bagnold number, which is the ratio of inertial to viscous forces, based on
CARPT data, indicates that it is of order 100. Hence, it is necessary to model explicitly both particle–particle and particle–fluid interactions. At Ba > 450 the particle motion is governed entirely by particle–particle interactions and at Ba < 40 by fluid viscosity, so our experiments are in the middle. Since rigorous theory that would account exactly for particle–liquid interactions due to mean velocity, particle interactions with fluctuating liquid velocity, mean velocity interactions with fluctuating ones for the solids and the liquid is not available, we have resorted to a judicious choice of empirical correlations. Solid phase fluctuations were modeled using the granular temperature concept. The governing equations are shown in Table 2. Partial justification for this is provided by the fact that experimentally observed p.d.fs of solids fluctuating velocity in various column compartments can be approximated with a Gaussian distribution [14]. The equation for the granular temperature is derived in the usual way by appropriate ensemble averaging and its right hand side
Table 1 Two fluid approach
755
Phases treated as interpenetrating continua Probability of occurrence of a given phase in multiple realizations of the flow is given by volume fraction of that phase Equations written for ensemble averaged realizations of each phase Both fluids are treated as incompressible o ðek qk Þ þ r ðek qk~ uk Þ ¼ 0 es þ ef ¼ 1; k ¼ f for fluid; k ¼ s for solids ot o ðef qf~ uf Þ þ r ðef qf~ uf ~ uf Þ ¼ ef rp þ r sf þ ef qf~ us ~ uf Þ þ ~ Ff g þ Ksf ð~ ot o ðes qs~ us Þ þ r ðes qs~ us ~ us Þ ¼ es rp rps þ r ss þ es qs~ uf ~ us Þ þ ~ Fs g þ Ksf ð~ ot
Table 2 Granular flow kinetic theory Inspired by kinetic theory of gases, phenomenological model for micro-scale interactions Thermodynamic temperature replaced by granular temperature, which is a flow property Closures for various terms in granular kinetic energy transport equation written in terms of granular temperature h ¼ 23ks
756
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
contains as the first term the production of kinetic energy by interaction of the normal and shear stress matrix with the mean velocity field by which energy is extracted from the mean field into the fluctuating field. The second term is the transport of kinetic energy by conduction, i.e. flow of energy caused by elastic and inelastic collisions between solid particles. The third term represents the dissipation of kinetic energy due to collisions while the fourth is the energy exchange due to interactions with liquid phase fluctuations. The solids pressure, shear viscosity and conductivity are calculated using the expressions of Chapman and Cowling [18] as modified by Ding and Gidaspow [19]. The radial distribution function used by Louge et al. [20] was employed. Collisional terms dominate the contribution for solids pressure and viscosity for the solids volume fraction encountered in our system, while both kinetic and collisional terms are important in determining solids conductivity. Energy dissipation through collisions was modeled as Chapman and Cowling [18] and the restitution coefficient is assumed very close to one (0.98) due to the damping effect of the liquid phase. The variation in this coefficient from 0.95 to 1 does not affect the results. The dissipation due to the friction of the fluctuating solids and liquid is modeled consistent to the drag formulation for the interaction of mean velocities. For liquid–solid drag we have used the correlation by Wen and Yu [21] based on systems with similar solids loadings. Finally to close the equations for liquid phase turbulence a modified k-e model is used that accounts for the presence of the dispersed phase [22]. Appropriate wall boundary conditions, explained in detail in the thesis of Roy [14] were used. Computation was performed in the framework of the Fluent two fluid model version 4. Both 2D axisymmetric and 3D simulations were performed. Best comparison with data was obtained with 3D simulations since the flow clearly has a 3D character. Such 3D simulations using 39,800 cells produce good agreement of calculated time averaged solids axial velocity profile and solids holdup profile with CARPT and CT data as shown in Fig. 9. Averaging is done from 25 to 100 s simulation time. Granular temperature predictions are also good except in the wall region. However the 3D simulation took more than a month on a DEC Alpha 500 MHz machine to execute. To provide faster input to reactor modeling, a 2D simulation coupled with a scalar tracer transport equation was executed and the results for the liquid and solids RTD are in agreement with experimental observation that liquid is close to plug flow and solids are much more backmixed (Fig. 10). The differences in the mean and variance between the computed and experimentally determined values are well within 10%. The following conclusions can be derived from this study: the CARPT-CT technique provides the needed
Fig. 9. Comparison of computed and measured values.
information on solids holdup, velocity and mixing for establishment of more accurate reactor models. While liquid in a liquid–solid riser is close to being in plug flow [14], solids are extensively backmixed. The Euler–Euler model is capable of computing reasonably well the overall features of the reactor scale liquid and solids flow via 2D axisymmetric simulation. The true transient nature of the flow is, however, only revealed by 3D simulations which produce results for time averaged velocity, holdup and kinetic energy in good agreement with experimental observations.
5. Bubble column flows Let us now examine a bubble column, which is the reactor of choice for conversion of natural gas and
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
Fig. 10. Computed RTD’s of the liquid and solids at UL ¼ 20 cm/s, S=L ¼ 0:15.
synthesis gas to clean fuels and is used in a great number of industrial applications. The flow is buoyancy driven as the gas superficial velocity exceeds by orders of magnitude the liquid superficial velocity. For very small catalyst particles the slurry phase can be treated as pseudohomogeneous. The key to proper modeling of the bubble column reactor is to understand and quantify the reactor scale flow pattern. Current models are based on either assumptions of ideal flow (e.g. plug flow for the gas and perfect mixing for the liquid) or on the ADM. However, no truly predictive correlations for the dispersion coefficient can be found. Also ignored is the fact that most columns need to operate in the churn turbulent regime at high gas velocities and high gas holdup. 5.1. Experimental data Using CT we have confirmed that gas holdup profiles are in the time averaged sense symmetric with more gas in the middle than by the walls while in bubbly flow the profiles are relatively flat. In churn turbulent flow the gas holdup profile becomes parabolic [7]. Using CARPT with a neutrally buoyant particle, we have confirmed the findings of Hills [23] and many others that there is a distinct time averaged liquid recirculation profile with liquid rising in the middle and falling by the walls. The time averaged liquid flow pattern is remarkably symmetric and forms a single recirculation cell [11]. However, in addition to the ensemble averaged velocity
757
profile, CARPT also provides eddy diffusivities in the radial and axial direction. The values in the axial direction are an order of magnitude larger than the values in radial direction [11]. Diffusivities increase with superficial gas velocity and with column diameter. Axial diffusivities are the highest in the region of downward time averaged liquid flow which can be explained by the tendency of bubbles to resist downward motion creating more backmixing of the liquid. As already mentioned, our CT and CARPT data reveal conclusively that a time averaged liquid recirculation profile exists through most of the column and is driven by a pronounced radial gas holdup profile [7,10,11]. The fact that even in churn turbulent flow the time averaged flow pattern is symmetric is clearly established. The liquid rises up in the middle and falls by the walls. The long term pattern is axisymmetric. However, at each point (cell, compartment) one can establish, based on CARPT data, a histogram of instantaneous velocities which yield zero means for radial and azimuthal velocities through most of the column [11]. Axial Lagrangian auto-correlation coefficient and axial–radial cross-correlation coefficient can also be established and are indispensable in evaluation of eddy diffusivities. Turbulent stresses for the 44 cm (1800 ) diameter column in churn turbulent flow at superficial gas velocity of 10 cm/s exhibit the usual trends for bubble columns [11,24]. The axial normal stresses are twice as high as the normal radial or azimuthal stresses. The only non-zero shear stress is the r–z one which peaks close to the velocity inversion point. It is important to note that turbulent stresses measured by CARPT are in good agreement with the values obtained by hot film anemometry and laser doppler anemometry under conditions of bubbly flow when the latter two techniques can be made to work [25,26]. Application of the R=S analysis to the CARPT data provides the Hurst exponents which indicate diffusion type mixing in the radial direction and strong convective effects in the axial direction [27]. Application of the Lagrangian and qualitative dynamics tools analysis to CARPT data clearly reveals the role of bubble wakes on the motion of the tracer particle and identifies and quantifies the chaotic fluid dynamics of the liquid phase. The effect of column diameter and operating conditions on the Kolmogorov entropy, strange attractors and kinematics of bubble wakes can then be established [28]. We have shown that the experimental information obtained by CARPT-CT can be utilized in an improved 2D phenomenological reactor model for the column presented schematically in Fig. 11 [29,30]. Based on rigorous arguments one can show that the time averaged velocity profile needed for the 2D scalar transport model is the ensemble averaged velocity determined from CARPT and that the needed liquid holdup profile is indeed the time averaged holdup determined by CT. The
758
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
Fig. 11. Suggested model for tracer spreading in the liquid phase of the bubble column.
assumption made is that locally the unknown crosscorrelation of the liquid velocity––tracer concentration fluctuations can be closed by the Bousinesq approximation i.e. by a product of the diffusivity tensor and the gradient of the tracer concentration. Moreover, it is assumed that diffusivities are the ones determined by CARPT. The usefulness of these assumptions is tested by feeding the CARPT measured velocity profile and CT measured holdup profile into the 2D convection diffusion model together with the CARPT measured diffusivities. One is now able to calculate the normalized exit mixing cup concentration for the impulse injection of tracer at the inlet and obtain the E-curve (exit age density function p.d.f. of the RTD). Comparison with the experimentally determined E-curve in response to the dye impulse injection is astonishingly good as shown in Fig. 12. This means that the model is capable of predicting reactor performance exactly for all linear kinetics. The success of this approach in being able to produce tracer response in agreement with measurement on a real hot pilot plant bubble column 1800 in diameter has been discussed elsewhere [29,30]. 5.2. Modeling We have qualitatively compared CARPT data (outline of particle trajectory and mean distance traveled) with the Euler–Lagrange simulation of bubbly flows [31]. More than 50,000 bubbles (assuming no bubble– bubble interactions) were simulated in a 3D simulation. Reasonable agreement between data and simulations was observed.
Since churn turbulent flow is of most current interest, Lagrange–Euler simulations are not practical, and the Euler–Euler interpenetrating two fluid model is a natural choice. The issues remain as to how to represent phase–phase interactions such as lift, added mass and drag and how to describe bubble induced turbulence. A simple approach that neglects lift, uses classical approach for added mass, models the drag force based on a sphere affected by presence of other spheres has been tried successfully in the Euler–Euler interpenetrating two fluid model to compute liquid velocities, gas holdup profiles, turbulent parameter and even diffusivities in bubble columns. The models were implemented and executed either in the CFDLIB code framework of Los Alamos or in Fluent 4.5 codes. Usually available closures, based on steady state expressions were used for drag and added mass. Liquid phase turbulence was accounted for either via the bubble induced turbulence [32,33] or via k–e model adjusted for the presence of the second phase [34]. Until recently an average bubble size was assumed for evaluation of the interfacial forces. At present we have successfully implemented the bubble population balance model with breakup coalescence models to predict the evolution of the bubble size distributions. Reasonable agreement between computed and experimentally observed results were reported for bubbly flow in quasi-2D columns [32,33] and for both bubbly and churn turbulent flow in 3D columns [34,35]. Pan et al. [33] clearly demonstrated that not only the time averaged liquid velocity profile was correctly calculated by the CFDLIB Euler–Euler code but that the right magnitude was predicted for liquid turbulent kinetic energies and stress profile in bubbly flows. Moreover, the calculated periods of enriched gas plume fluctuation agreed well with the experimentally observed ones. This provided confidence that the Euler–Euler two fluid model is capable of calculating liquid velocities in bubbly flow within engineering accuracy. In churn turbulent flow, 2D axisymmetric calculations produced liquid time average velocity in reasonable agreement with the experimental CARPT data [34,35]. Correct overall gas holdup magnitudes were calculated but computed gas holdup profiles did not exhibit the CT experimentally observed parabolic radial profile. 3D computations, which take a long time to execute, came closer to producing the observed holdups but substantial discrepancies still exist. Work is in progress to determine whether the incorporation of the population balance model and coalescence break up models can alter sufficiently the computed holdup profile. Most importantly, from the reaction engineering viewpoint, it has recently been shown [35] that the CFDLIB Euler–Euler model with the above discussed closures can, via 3D computations, produce the eddy diffusivities whose radial profile is in
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
759
Fig. 12. Prediction of non-volatile tracer impulse response.
close agreement with values obtained from CARPT experiments. Hence, CFD offers possibilities in calculating the parameters needed for reactor modeling, such as time averaged velocity profile and eddy diffusivities. As mentioned earlier a 2D convection–diffusion model that utilizes CFD calculated flow parameters can be quite useful in predicting tracer responses and reactor performance of industrial size reactors.
and in describing the liquid flow field in churn turbulent bubble columns is established. The capabilities of the Euler–Euler interpenetrating two fluid model in generating flow pattern information in agreement with CARPT-CT data that can form a basis for reactor modeling are also shown and are significant for industrial applications.
7. Concluding remarks 6. Practical significance/usefulness This paper describes the key features of and provides the key references for CARPT-CT. The techniques are unique in their capabilities of providing the flow field description in large enclosures for opaque systems. Their utility in obtaining the solids holdup distribution, velocity and kinetic energy profiles in a liquid–solid riser
The CARPT-CT techniques provide us with the unique opportunities to assess holdup and velocity fields in opaque two phase systems. The resolutions of both techniques allow reliable information to be generated for large scale motion at lower frequencies. While CT provides the time averaged density distribution, CARPT, in addition to time averaged velocities, allows estimation
760
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761
of a number of turbulent quantities, e.g. kinetic energy, stresses, eddy diffusivities, etc. The wealth of information provided by CARPT-CT allows us to quantify the flow fields in systems like liquid–solid risers and bubble columns which leads to improved and more detailed reactor models. CFD computations, based on the Euler–Euler interpenetrating two fluid model have been shown to be capable of producing information in good agreement with data, which renders hope for CFD based scale-up in the future. It is fair to conclude that CARPT-CT are indeed a valuable measurement tool for opaque multiphase systems and that the Euler–Euler model has potential in simulation of reactor scale multiphase flows at conditions of industrial interest.
Acknowledgements At the end I would like to thank the many collaborators that have produced all the results presented here. Above all, Narsi Devanathan, Yubo Yang and Sailesh Kumar deserve credit for designing, constructing and implementing the hardware and software for CARPTCT in CREL. Shantanu Roy and Abdenour Kemoun executed the experiments in the liquid–solid riser discussed here, and Shantanu modeled them using Fluent with the help of Jay Sanyal (Fluent). Sue Degaleesan, Jinwen Chen, Booncheng Ong experimented with the bubble column, while Yu Pan and S. Roy solved CFD models with CFDLIB and Fluent. Muthanna Al-Dahhan offered helpful suggestions and Bernard Toseland of Air Products was instrumental in focussing us to meet industrial needs. Finally, we are thankful to our financial sponsors: Exxon-Mobil for donating the CARPT-CT elements, Department of Energy and all CREL industrial sponsor companies from five continents (ABB Lummus, Air Products, Bayer, Chevron, Conoco, Dow Chemicals, DuPont, Elf Atofina, ENI Technologies, Exxon-Mobil, IFP, Intevep, Mitsubishi, Praxair, Sasol, Shell, Solutia, Statoil, Synetix-ICI, UOP). Their financial support made this work possible.
References [1] P.L. Mills, Recent development and advances in laboratory-scale multiphase reactors, Plenary Lecture, First North American Symposium on Chemical Reaction Engineering (NASCRE1), Houston, TX, 2001. [2] G.F. Froment, K.B. Bischoff, Chemical Reactor Analysis and Design, second ed., Wiley, NY, 1990. [3] J. Chaouki, F. Larachi, M.P. Dudukovic, Non-invasive tomographic and velocimetric monitoring of multiphase flows, Ind. Eng. Chem. Res. 36 (11) (1997) 4476–4503.
[4] J. Chaouki, F. Larachi, M.P. Dudukovic (Eds.), Non-invasive monitoring of multiphase flows, Elsevier, Amsterdam, 1997. [5] S.B. Kumar, D. Moslemian, M.P. Dudukovic, A c-ray tomographic scanner for imaging of void distribution in two-phase flow systems, Flow Meas. Instrum. 6 (3) (1995) 61–73. [6] S.B. Kumar, M.P. Dudukovic, Computer assisted gamma and x-ray tomography: Application to multiphase flow systems, in: Chaouki et al. (Eds.), Non-Invasive Monitoring of Multiphase Flows, Elsevier, Amsterdam, 1997, pp. 47–103, Chapter 2. [7] S.B. Kumar, D. Moslemian, M.P. Dudukovic, Gas holdup measurements in bubble columns using computed tomography, AIChE J. 43 (6) (1997) 1414–1425. [8] J.S. Lin, M.M. Chen, B.T. Chao, A novel radioactive particle tracking facility for measurement of solids motion in gas fluidized beds, AIChE J. 31 (3) (1985) 465–473. [9] N. Devanathan, Investigation of liquid hydrodynamics in bubble columns via computer-automated radioactive particle tracking (CARPT), D.Sc. Thesis, Washington University, St. Louis, MO, 1991. [10] N. Devanathan, D. Moslemian, M.P. Dudukovic, Flow mapping in bubble columns using CARPT, Chem. Eng. Sci. 45 (8) (1990) 2285–2291. [11] S. Degaleesan, Fluid dynamic measurements and modeling of liquid mixing in bubble columns, D.Sc. Thesis, Washington University, St. Louis, MO, 1997. [12] F. Larachi, G. Kennedy, J. Chaouki, A c-ray detection system for 3-D particle tracking in multiphase reactors, Nucl. Instrum. Meth. A 338 (2-3) (1994) 568–576. [13] P. Gupta, Bubble column dynamics: experiments and models, D.Sc. Thesis, Washington University, St. Louis, MO, May, 2002. [14] S. Roy, Quantification of two phase flows in liquid–solid risers, D.Sc. Thesis, Washington University, St. Louis, MO, December, 2000. [15] S. Roy, M.P. Dudukovic, Flow mapping and modeling of liquidsolid risers, Ind. Eng. Chem. Res. 40 (23) (2001) 5440–5454. [16] S. Roy, A. Kemoun, M. Al-Dahhan, M.P. Dudukovic, A method for estimating the solids circulation rate in a closed-loop circulating fluidized bed, Powder Technol. 121 (2–3) (2001) 213–222. [17] D.A. Drew, Mathematical modeling of two phase flow, Ann. Rev. Fluid Mech. 15 (1983) 261–291. [18] S. Chapman, T.G. Cowling, Mathematical Theory and Nonuniform Gases, third ed., Cambridge University Press, Cambridge UK, 1990. [19] J. Ding, D. Gidaspow, A bubbling fluidization model using kinetic theory of granular flow, AIChE J. 36 (4) (1990) 523–538. [20] M.Y. Louge, E. Mastorakos, J.T. Jenkings, The role of particle collisions on pneumatic transport, J. Fluid Mech. 231 (1991) 345– 359. [21] C.Y. Wen, Y.H. Yu, Mechanics of fluidization, Chem. Eng. Progr. Symp. Sci. 62 (62) (1966) 100–111. [22] S.E. Elghobashi, T.W. Abou-Arab, A two-equation turbulence model for two-phase flows, Phys. Fluids 26 (4) (1983) 931–938. [23] J.H. Hills, Packed nonuniformity of velocity and voidage in a bubble column, Trans. Inst. Chem. Eng. 52 (1974) 1–9. [24] J. Chen, L. Fan, S. Degaleesan, P. Gupta, M.H. Al-Dahhan, M.P. Dudukovic, B.A. Toseland, Fluid dynamic parameters in bubble columns with internals, Chem. Eng. Sci. 54 (13/14) (1999) 2187– 2197. [25] T. Menzel, T. Inderweide, O. Standacher, O. Wein, U. Onken, Reynolds shear stress for modeling of bubble column reactors, Ind. Eng. Chem. Res. 29 (6) (1990) 988–994. [26] R.F. Mudde, J.S. Groen, H.E.A. van den Akker, Liquid velocity field in bubble columns: LDA experiments, Chem. Eng. Sci. 52 (21/22) (1997) 4217–4224. [27] Y.B. Yang, N. Devanathan, M.P. Dudukovic, Liquid backmixing in bubble columns via computer-automated radioactive particle tracking (CARPT), Exp. Fluids 16 (1) (1993) 1–9.
M.P. Dudukovic / Experimental Thermal and Fluid Science 26 (2002) 747–761 [28] M. Cassanello, F. Larachi, A. Kemoun, M. Al-Dahhan, M.P. Dudukovic, Inferring liquid chaotic dynamics in bubble columns using CARPT, Chem. Eng. Sci. 56 (21/22) (2002) 6125–6134. [29] S. Degaleesan, S. Roy, S.B. Kumar, M.P. Dudukovic, Liquid mixing based on convection and turbulent dispersion in bubble columns, Chem. Eng. Sci. 51 (11) (1996) 2721–2725. [30] S. Degaleesan, M.P. Dudukovic, B.A. Toseland, B.L. Bhatt, A two compartment convective diffusion model for slurry bubble column reactors, Ind. Eng. Chem. Res. 36 (11) (1997) 4670–4680. [31] N. Devanathan, M.P. Dudukovic, A. Lapin, A. Lubbert, Chaotic flow in bubble column reactors, Chem. Eng. Sci. 50 (16) (1995) 2661–2667.
761
[32] Y. Pan, M.P. Dudukovic, M. Chang, Dynamic simulation of bubbly flow in bubble columns, Chem. Eng. Sci. 54 (13/14) (1999) 2481–2489. [33] Y. Pan, M.P. Dudukovic, M. Chang, Numerical investigation of gas driven flow in 2-D bubble columns, AIChE J. 46 (3) (2000) 434–449. [34] J. Sanyal, S. Vasquez, S. Roy, M.P. Dudukovic, Numerical simulation of gas–liquid dynamics in cylindrical bubble column reactors, Chem. Eng. Sci. 55 (21) (1999) 5071–5083. [35] Y. Pan, M.P. Dudukovic, CFD simulations of bubble columns, 6th World Congress of Chemical Engineering, Melbourne, Australia, September, 2001.