A hybrid CFD—reaction engineering framework for multiphase reactor modelling: basic concept and application to bubble column reactors

Chemical Engineering Science 58 (2003) 3077 – 3089

www.elsevier.com/locate/ces

A hybrid CFD—reaction engineering framework for multiphase reactor modelling: basic concept and application to bubble column reactors Stelios Rigopoulos, Alan Jones∗ Department of Chemical Engineering, University College London, Torrington Place, London WC1E 7JE, UK Received 6 November 2002; received in revised form 2 April 2003; accepted 3 April 2003

Abstract The coupling of turbulent mixing and chemical phenomena lies at the heart of multiphase reaction engineering, but direct CFD approaches are usually confronted with excessive computational demands. In this hybrid approach, the quanti5cation of mixing is accomplished through averaging the 6ow and concentration pro5les resulting from a CFD 6ow 5eld calculation and a computational (“virtual”) tracer experiment. Based on these results, we construct a mapping of the CFD grid into a generalised compartmental model where the chemistry calculations can be e:ciently carried out. In contrast to the empirical models used in the residence time distribution (RTD) approach, the compartmental model in this methodology, owning to its CFD origins, retains the essential features of the equipment geometry and 6ow 5eld. A procedure for extracting the mixing information from k– based CFD codes is outlined, but the main concept of the approach is not restricted to any particular type of turbulence modelling, and will therefore bene5t from future developments. A phenomenological model of mass transfer and chemical reaction, based on the penetration theory, is employed to simulate the interfacial phenomena in gas–liquid reactors, and a study of CO2 absorption into alkali solution is presented to demonstrate the method. ? 2003 Elsevier Ltd. All rights reserved. Keywords: Mixing; Reaction engineering; Computational 6uid dynamics

1. Present challenges in multiphase reactor modelling The performance of gas–liquid reactors relies on the combined outcome of multiphase 6uid dynamics, interfacial mass transfer and chemical reactions. The detailed, simultaneous analysis of the above phenomena is out of reach given the currently available computing power, partly because they take place at di?erent scales. Furthermore, most of these phenomena are non-linear, and empirical models cannot be extrapolated too far from the set of design parameters and operating conditions under which they were extracted. As a result, multiphase reaction engineering is still confronted with major challenges from a modelling point of view. The quantitative description of mixing is a major issue of reaction engineering. Early research attempted to quantify mixing without knowledge of the 6ow 5eld, and a major step in this direction was taken when Danckwerts (1953, 1958) introduced the concept of residence time ∗

Corresponding author. Present address: Department of Mechanical Engineering, Imperial College London, Exhibition Road, London SW7 2BX, UK. E-mail address: [email protected] (A. Jones). 0009-2509/03/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0009-2509(03)00179-9

distribution (RTD). The RTD arising from tracer experiments can be used to identify situations that lie between the two extremes of plug 6ow and complete mixing, and a variety of compartmental models, such as the tank-in-series and the cell model with back6ow (or its continuous equivalent, the axial dispersion model) have been proposed to accommodate them (Mecklenburgh & Hartland, 1975; Nauman & Bu?ham, 1983). Also notable is the “network-of-zones” approach of Mann and Hackett (1988), a semi-empirical alternative to CFD. More fundamental information about mixing can be obtained through the Navier–Stokes equations which, though, require an excessive amount of computation to resolve all turbulence scales by direct numerical simulation (DNS). Shortcuts with variable levels of detail have been developed, ranging from the commonly used in engineering turbulent viscosity-based approaches such as the k– (Launder & Spalding, 1972) to the much more accurate Large Eddy Simulation (LES, see e.g. Pope, 2000), which, though sparingly used today, may become the industry standard in the future. The coupling of 6uid dynamics with complex and often non-linear chemical reaction systems, however, results in a steep increase of computational time, that

3078

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089

often prohibits the use of this approach for reactor scale-up purposes. In the last decades, the reaction engineering community has started to explore the prospect of employing CFD in the modelling of multiphase reactors. Torvik and Svendsen (1990) made one of the 5rst attempts to model mixing and reaction in a bubble column with CFD, based on the k– model, while mass transfer and chemical reaction were still dealt with by assuming constant reaction enhancement. Ranade (1995), Delnoij, Kuipers, and van Swaaij (1997), Van den Akker (1997), Kuipers and van Swaaij (1997), and Dudukovic, Larachi, and Mills (1999) among others, have reviewed the current state of the art regarding the application of CFD to chemical reaction engineering. The reviews of Jakobsen, Sannaes, Grevskott, and Svendsen (1997) and Joshi (2001) focus on CFD modelling of gas–liquid 6ows, with relevance to bubble column reactors. The direct application of CFD to chemical processes faces several obstacles, however. Even in single phase reactors, chemical reactions are described by highly non-linear terms that often cause numerical instabilities. More complex phenomena, such as interfacial mass transfer with chemical reaction, are encountered in multiphase or catalytic reactors, while in cases such as precipitation and biochemical engineering additional equations, such as the population balance, must be solved alongside the Navier–Stokes equations. The computational resources required are often prohibitive, while the complexity of the problems that arise from the coupling of the 6uid dynamics with the chemical phenomena means that the systems must be treated case-by-case. Recently, “hybrid” approaches have emerged as an alternative. In those CFD is employed only for hydrodynamic simulation, while the chemical phenomena are resolved in a custom-built compartmental model. Although this decoupling cannot be applied to cases where the coupling between hydrodynamics and chemistry is very strong, such as in combustion, many chemical reactors with a liquid bulk phase fall into this category. Bauer and Eigenberger (1999) used a “zone model” to study a bubble column reactor; Bezzo, Macchietto and Pantelides (2000) developed an interface of communication between the gPROMS modelling software and a commercial CFD code; Zauner and Jones (2000) used a segregated feed model in conjunction with CFD to study precipitation in a stirred tank. Though these approaches represent essentially a compromise, their lifetime is likely to be quite long, as the full CFD solution would require a radical breakthrough in computer technology. The approach proposed here belongs to this class. In its context, the detailed hydrodynamic information on mesoand macro-mixing that CFD provides is mapped into a more tractable compartmental model, suitable for e:cient calculation of chemical reactions and other non-linear phenomena (e.g. precipitation), while still incorporating the essential features of the geometry and operating conditions of the reactor. The main concept of the approach, the mapping of the CFD grid into the generalised compartmental model,

is introduced in Section 2. This is followed by a procedure for linking it to k–-based CFD, and an application to bubble column reactors by superimposing a phenomenological model of mass transfer with chemical reaction. The basic methodology is not limited to these, however, and may be used in conjunction with more advanced turbulence models, or other phenomenological reaction engineering models. 2. The generalised compartmental model 2.1. Basic concepts In this work, we shall distinguish between mixing due to the mean (time-averaged) 6ow 5eld, and turbulent local mixing, i.e. any localised mixing e?ects, caused by turbulence at scales higher than the Kolmogorov. To describe the former, we shall use the mean 6ow information that is yielded by a CFD code supplied with a turbulence model based on time averaging, such as the k–. To account for the later we shall employ a “virtual” tracer experiment, i.e. compute the dispersion of a passive scalar using the above CFD code. The 6ow and scalar concentration 5elds that result from CFD will subsequently be used to derive a compartmental model, on which the chemical reactions can be superimposed provided that they have no signi5cant feedback on the hydrodynamics. For this purpose, we shall introduce the generalised compartmental model, a network of compartments with many degrees of freedom that may be tailored to accommodate the main features of the 6ow 5eld and turbulent mixing. The outcome of an actual tracer experiment is a response curve that, though indicative of processes that take place (e.g. bypassing, short-circuiting), does not yield insight into the features of the geometry and 6ow that cause them; different systems can produce the same RTD. On the contrary, the objective of the “virtual” tracer experiment proposed here is to reproduce the detailed concentration distribution that would result from the mixing conditions that prevail during the actual operation of the reactor. The information thus yielded is of a much more mechanistic nature, and the compartmental model derived on its basis has the main features of the 6ow 5eld and turbulent exchange inherent in it. 2.2. Model topology The generalised compartmental model is a projective mapping of the CFD grid into a coarser network of fully mixed compartments. The exact shape of the compartments and their boundaries is of no essence, since they are spatially uniform; instead, all that matters is the way they are connected. The concept is best illustrated through the example of Fig. 1, where 5ve subdomains of the CFD grid are mapped into a compartmental model that retains only the 6uxes between the original subdomains. The compartmentalisation must be based on physical reasoning, and will be discussed in more detail later; at the moment it will be

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089

2

1

3

5

1

consequence, both 2D and 3D CFD grids may be mapped in the same way. To describe multi-phase reactors, separate compartments are employed for each phase. In this work the structure of the network will be the same for all phases, to allow the implementation of mass exchange between the phases via the penetration model. Thus the building block of a two-phase compartmental network comprises the two compartments of each phase, in contact with up to L other compartments (Fig. 2). The following notation is introduced:

2

3

Cikp

4 4

5

Vkp in Qkp Fig. 1. Mapping of a CFD grid into a compartmental network.

out Qkp

pointed out that in Fig. 1 compartment 5, entirely enclosed within 4, may be devised to contain a vortex that results in dead space. It is important to stress the di?erence in the physical interpretation of a compartment between the RTD models and the present approach. Whereas in RTD models a compartment arises from the 5tting of an RTD curve but does not necessarily correspond to a physical part of the reactor, in the present approach the compartments are generated from the CFD grid, and are thus tightly bound to the geometry of the particular piece of equipment. Let K be the number of compartments in the network. Each compartment will be identi5ed by the compartment index, k: k ∈ {1; 2; : : : ; K};

(1)

Each compartment is allowed to have boundaries with a maximum of L neighbours, the indices of which are established with the mapping: I : (k; l) → I (k; l):

(2)

With each compartment k we associate a maximum of L connectivity 6owrates, denoting 6ow towards its lth neighbour. Consequently, a pair of opposite direction 6owrates is associated with each boundary. This two-way connectivity provides a means of introducing local turbulent exchange. The inverse map of I (k; l) yields the indices of the L (at most) compartments that have an out6ow towards k: I
: (k; l) → I
(k; l):

3079

(3)

The one-dimensional array of compartment indices, together with the maps I (k; l) and I
(k; l) establish the topology of the network, i.e. the way in which the cells are connected. Only I (k; l) is required to specify the network; I
(k; l) may be easily determined by the program. The dimensionality of the CFD grid space is not retained in the compartmental model: each compartment is simply aware of its neighbours’ identities (via I (k; l) and I
(k; l)), regardless of whether they are in the x, y or z direction. As a

I (k; l) Qkp

QIk
(k; l); p

concentrations of chemical species i, spatially averaged over compartment k at phase p volume of compartment k at phase p in6ow into the system (i.e. the whole reactor) at compartment k, at phase p out6ow from the system at compartment k, at phase p 6ow from compartment k, at phase p, towards its lth neighbour, the compartment I (k; l) 6ow into compartment k, at phase p, from its lth neighbour, the compartment I
(k; l)

The total number of parameters thus are: • number of cells; • volumes of cells (K); • connectivity 6owrates, maximum [K × L]. The total number of parameters that may be arbitrarily speci5ed (i.e. degrees of freedom) is however constrained by the mass balances, as it will be shown in the next section. 2.3. Constraints All of the parameters mentioned above may not be arbitrary speci5ed. To begin with, the number of cells is inferred from the speci5cation of the cell volumes. Moreover, the cell volumes in a phase must sum up to the total volume of the system in that phase: Vp =

K


Vkp :

(4)

k=1

For every cell in a phase the mass balance must be satis5ed, resulting in K further constraints: L
l=1

I (k; l) out Qkp + Qkp =

L


in QIk
(k; l)p + Qkp :

(5)

l=1

2.4. Mass balances A compartment k has 6owrates to and from L neighbouring cells, and there may also be input or output from the whole system at a location in the vicinity of k (Fig. 2). In the following, i stands for the species index. Mass balances

3080

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089 out

Qk1 : Inflow into the reactor at compartment k

out

Qk2

Compartment k

I(k,l) Qk1

L flowrates : from compartment k to its l-th neighbour

I(k,l)

Phase 1

k

L flowrates QI’(k)1 : from the l-th neighbour to compartment k

Phase 2

Qk1

k

QI’(k)2

in

Qk1 : Inflow into the reactor at compartment k

in

Qk2

Fig. 2. The building block of the multi-phase generalised compartmental model (two-phase).

take the form:



dCikp (t) 1 out = −Qkp Cikp (t) − dt Vkp in in Cikp (t) + Qkp

+

L


L
l=1

is usually the case, as long as the reactions are not highly exothermal and the reaction product does not a?ect the mixture rheology (e.g. polymerisation, foaminess).

I (k; l) Qkp Cikp (t)



3.1. Overall tuning algorithm

QIk
(k; l)p CiI
(k; l)p (t)

l=1

+fi (C1 (t); C2 (t):::CN (t));

(6)

where fi is the reaction term. This balance may be written more concisely as dCikp (t) = Wout Cikp (t) + Win dt +fi (C1 (t); C2 (t):::CN (t)); where 1 Wout = − Vklp

 L
l=1

(7)

 I (k; l) Qkp

+

out Qkp

;

L 1
k 1 in Win = + QI (k; l)p CiI (k; l)p (t) + C in (t) Qkp : Vkp Vklp ikp l=1

Note that, in the case of the gas phase, the reaction term may include the consumption due to chemical absorption. As this is a system of ordinary di?erential equations, initial conditions must be speci5ed for its solution: 0 t = 0 → Cikp = Cikp :

(8)

3. Link with CFD We shall now describe how the generalised compartmental model may be reduced to a speci5c lumped model by using the results of a CFD simulation. It must be noted that the implicit assumption in the present approach is that the coupling between the hydrodynamics and the chemical reactions is one-way, i.e. the reactions have no signi5cant feedback on the hydrodynamics. In gas–liquid reactors this

The overall process of constructing and 5ne-tuning the generalised compartmental model can be summarised by the following algorithm: Step 1: Carry out the 6ow 5eld computation and virtual tracer experiment with CFD. Step 2: Determine the network structure and estimate connectivity 6owrates based on mean 6ow. Step 3: Introduce two-way 6owrates to account for local turbulent exchange. Step 4: Evaluate the convergence between the reaction engineering model and the CFD tracer experiment. If poor, return to step 3. If not possible at all, return to step 2 and modify the structure. Step 5: Having established a network that adequately describes the dispersion, superimpose the reaction engineering models to account for complex reaction systems, gas–liquid mass transfer, micromixing, etc. 3.2. Extraction of mixing information from k– based CFD In principle, the results from any CFD simulation may be averaged to parametrise the generalised compartmental model. The k– model (Launder & Spalding, 1972) belongs to the family of turbulence models based on the Reynolds-averaged Navier–Stokes equations, whose objective is to produce an accurate prediction of the mean 6ow 5eld. This 5eld encompasses the gross features of the 6ow such as recirculation zones, large-scale vortices giving rise to dead space, bypassing routes etc, that give rise to macromixing. Further, localised mixing e?ects are revealed when one computes the 5eld of a scalar quantity.

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089

In order to establish a correspondence between the hydrodynamic information yielded by the k– and the mixing concepts used in reaction engineering, we shall distinguish between mixing due to the mean (time-averaged) 6ow 5eld, and local turbulent exchange, i.e. any localised mixing effects, lateral or axial, caused by turbulence at scales higher than the Kolmogorov. The former is fully described by the mean 6ow 5eld, while the later can be revealed through CFD by the computational equivalent of a tracer experiment. No mention is made to micromixing, which is not at all captured by the k– model. Its incorporation would require further phenomenological models, such as those of Baldyga, Bourne, and Hearn (1997). 3.3. The virtual tracer experiment In RTD studies the tracer is introduced as an impulse, step or periodic feed. In contrast, in the context of this approach the virtual tracer experiment must emulate the actual feeding policy of the reactor. The reason is to maintain the structure of the compartmental model as simple as possible, so that it is su:cient to account for the e?ects that determine reactor performance while attaining a maximum reduction in computational cost. Not all details of the 6ow and turbulence are necessary for that purpose. A point inlet, for instance, would introduce steep gradients around it, requiring several small compartments to accommodate it. A step feed, on the other hand, would cause a uniform concentration over a larger region. By introducing the tracer in a way that emulates the actual feeding policy, we ensure that the compartmentalisation process will result in the smallest number of compartments necessary to describe adequately the macromixing and turbulent exchange during actual reactor operation. In the case of a bubble column, the feeding policy can be emulated by introducing the passive scalar continuously throughout the column, at a value proportional to the local interfacial area, which in turn can be computed by the local volume fraction and mean bubble diameter (db ): [V of cell] S˙ ∼ [interfacial area per bubble]: (9) [N of bubbles] It must be noted here that this methodology employs mean 6ow information for the estimation of the inter-compartmental 6ows, in contrast to earlier hybrid models that relied on predicted turbulent quantities such as turbulent kinetic energy and dissipation rate. As Montante, Lee, Brucato, and Yianneskis (2001) have shown, these quantities are severely underestimated by k–-based simulations. 3.4. Determination of network structure The mapping of the CFD grid into the compartmental model can be made on the basis of a qualitative inspection of the 6ow and concentration 5elds. A potential objective for future work could be the automation of this stage, e.g. by

3081

the use of arti5cial intelligence tools. A few heuristic rules can be suggested: • The 6ow should not change direction along a boundary. • Streamlines should cross similar distances within a certain compartment. • Regions of uniform concentration may be accommodated by a single compartment. The implication of the 5rst rule is that the all main 6ow features giving rise to macromixing, such as circulation loops, must be accounted for by the initial structure (apart from small vortices, which may be fully contained within a single compartment). The second rule suggests that all 6uid elements entering a compartment will have similar residence times at that compartment. The third rule suggests that turbulent exchange does not have to be accounted for by the inter-compartmental 6owrates alone, but may be inherent in the network itself. Tailoring the network to the steady–state concentration 5eld is likely to lead to a structure that is insu:cient to account for the dynamic behaviour during start-up, which might require a more complex structure. If start-up is of importance, the designer may still initialise the structure according to the steady state, and subsequently re5ne it by comparing its dynamic performance to the CFD. The 6owrate at any compartmental boundary can be calculated from the CFD results by summing up the 6ows through the CFD cell faces that comprise that boundary. Such information is readily available in most CFD codes, as the 6ux term (ui dAi ) (where Ai is the area multiplied by the unit vector in the ith direction, and i a su:x of Cartesian coordinates) must be calculated to form the convective part of the coe:cient matrix. The network thus generated features one-way connectivity; the two-way is reserved for the introduction of turbulent mixing. Before proceeding, however, it is useful to simulate the virtual tracer experiment with the one-way network. Acceptable convergence would signify that turbulent local mixing is fully contained within the boundaries of the compartments. The introduction of two-way connectivity to account for turbulent exchange may be guided by the distribution of turbulent kinetic energy and dissipation rate. High turbulence around an intercompartmental boundary would signify that further material interchange, not accounted for by the mean 6ow, is bound to take place between the compartments involved. Several re5nements of the turbulent exchange 6owrates may be required until the virtual tracer experiment can be adequately simulated with the compartmental model. If that is not at all possible, it signi5es a fundamental inadequacy of the compartmental structure, and that the network should be constructed anew. 3.5. Superposition of reaction engineering models Once a compartmental network that properly accounts for mixing has been established, various mechanistic reaction

3082

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089

engineering models can be superimposed on it. Although the present paper focuses on bubble column reactors, the main concept lends itself to other potential reaction engineering applications, such as: • Single-phase reactors with complex non-linear reaction systems, whose direct coupling with CFD would pose serious computational problems. • Stirred tank reactors with fast reactions, where micromixing models can be superimposed in the compartments located at the inlets. • Precipitation, where additional population balance equations describing crystal growth, agglomeration or breakage must be solved. • Multiphase non-catalytic reactors (e.g. bubble columns), where interfacial mass transfer and chemical reaction models can be accommodated within the compartments. These models must be resolved at the interfacial scale, and therefore cannot be directly incorporated in a CFD code. • Catalytic and biochemical reactors, where adsorption models or cell population balances may be needed. Some of the reaction engineering models required for the above applications require extra hydrodynamic information: e.g. a gas–liquid mass transfer model requires knowledge of gas hold-up or a particle agglomeration model (in precipitation) may be a function of the energy dissipation rate. This extra information may be obtained from the CFD simulations in the same way that the mixing information was extracted, i.e. by averaging the distributed information yielded by CFD within the compartments. The present work focuses on bubble column reactors, and a way of introducing interfacial mass transfer and chemical reactions in the network according to the penetration theory will be demonstrated. This way, the absorption enhancement due to a fast, but not instantaneous, chemical reaction system can be properly taken into account. 4. Application to bubble column reactors 4.1. Mixing in bubble column reactors The mixing events that will be simulated by the network are liquid circulation due to forced or natural convection (linked to the mean 6owrates), and turbulent exchange between regions with high gas holdup (where the reactant is introduced) with regions with low holdup (linked to the virtual tracer experiment). Furthermore, a phenomenological model of mass transfer with chemical reaction will be superimposed, to provide the link of the compartments of the liquid phase with those of the gas phase. 4.2. Model of interfacial mass transfer and chemical reaction Unless the time scales of a gas–liquid reaction and mass transfer di?er by orders of magnitude (instantaneous or

slow reaction regime), these two phenomena are intricately coupled. The consumption of absorbent at the interface increases the driving force for mass transfer, a phenomenon known as reaction enhancement, and the amount of absorption can only be resolved by detailed solution of the di?usion–reaction equations that result from a phenomenological model such as the 5lm, penetration and surface renewal. Steady-state processes can be characterised by an enhancement factor, which may be calculated or extracted experimentally, but in a dynamic process the enhancement is bound to vary, depending on the concentrations at every instance. Recently several dynamic models have appeared, based on either the 5lm (Romanainen, 1997) or the penetration theory (HostomskWy & Jones, 1995; Van Elk, Borman, Kuipers, & Versteeg, 1999; Rigopoulos & Jones, 2001). Among the established absorption models the penetration theory (Higbie, 1935) will be employed here, in a way similar to the works cited above. It is conceptually suited to the mechanism of bulk–interface interaction in bubbly 6ows: a bubble moving with a slip velocity relative to the bulk is in contact with a surrounding element of bulk 6uid for a short amount of time, and subsequently overcomes it. This process is emulated through the following scheme: the liquid element surrounding the bubble and the bulk are considered as two separate dynamic reactors that operate independent of each other and interact at discrete time intervals. In the beginning of the contact time, the interface is detached from the bulk. Subsequently, the equations describing the bulk and interface are solved in the bulk and interfacial reactors, respectively. When overcome by the bubble, the liquid element returns to the bulk and is presumed to get mixed instantaneously. Not all the liquid elements within a cell will contact the gas simultaneously, but since the contact time is very small compared to the process time the situation is equivalent to having a very large interfacial element emerging from the bulk, returning to be mixed and then re-emerging in the surface. The new concentrations that arise from the mixing are the initial conditions for the next contact time integration. 4.3. Interface equations Penetration theory assumes unsteady–state di?usion without convection on one spatial dimension (note that the notation of the compartmental model, such as indices referring to compartments and phases, does not apply to this section— the interfacial element is treated in isolation): @Ci (x; t) @2 Ci (x; t) =D @t @x2 +fi (C1 (x; t); (C2 (x; t); : : : ; (CN (x; t));

(10)

where fi stands for the non-linear reaction term. At the gas side of the interface (x=0) the concentration of the absorbed component can be explicitly speci5ed by an equilibrium relation such as Henry’s law, while zero 6ux is the boundary condition for the non-volatile components:

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089

Volatile species: x = 0;

t ¿ 0 ⇒ Ci (x; t) =

yP : H

(11)

Non-volatile species: dCi (x; t) x = 0; t ¿ 0 ⇒ = 0: (12) dx Regarding the boundary conditions at the bulk side we adopt the concept of Toor and Marchello (1958), according to which they can be speci5ed at a 5xed distance (x = ) from the interface, in contrast to the cumbersome (x → ∞) proposed by the original Higbie model. This is valid on the condition that concentration gradients are zero approaching the bulk at the given distance, according to the criterion ¡ 2 =D, easily satis5ed for 1–10 mm bubbles with a terminal slip velocity of 20 –30 cm s−1 , and an interface depth in the range of 10 –100 m (di?usion is in the order of 10−9 m2 s−1 ). Note that the contact time, can be calculated as the ratio of the bubble diameter and slip velocity. The initial conditions and the boundary conditions at the bulk side are speci5ed as the bulk concentrations at the beginning of the contact time: x = ;

t ¿ 0 ⇒ Ci = Ci; bulk ;

(13)

t = 0;

x ¿ 0 ⇒ Ci = Ci; bulk :

(14)

The concentration distribution in the interface, at the end of the contact time can be calculated as   = Ci (x; ) d x: (15) [Ci ]int tot 0

From this, and taking into account the stoichiometry of the reaction, the amount of absorption can be deduced. 5. Computational implementation 5.1. Generalised compartmental model In our implementation, the set of ODEs (Eqs. (7)) that describe the bulk of each compartment are integrated with the LSODE set of subroutines, employing an Adams method. LSODE is available in the public domain via www.netlib.org, and is discussed by Hindmarsh (1983). Successful tuning of the generalised compartmental model to the CFD requires, as it will be shown in the next section, consideration of many di?erent structures for successive re5nement. Hence the success and usability of the method relies on an e:cient computational implementation that facilitates the generation of structures as much as possible. In our code, all network characteristics are de5ned through the input 5les, in a comprehensive way that allows quick and easy modi5cation. Each structure is checked for satisfaction of the mass balance (Eq. (5)) before proceeding. The code has been programmed in FORTRAN 90, which allows elementary object oriented programming through the

3083

use of derived data types (structures). A structure is constructed to describe the network compartments, while all properties associated with them (volumes, 6owrates, concentrations, etc.) are de5ned as elements of the structure, corresponding to the matrix notation employed here. Only a few seconds are required for simulating tracer experiments, as these involve just convection of material without reaction terms. Once the network has been tuned to the CFD results, the reactions can be superimposed. Computations with reaction are still very fast; about 30 compartments were employed for the case study presented here compared with thousands of CFD cells. 5.2. Solution of interfacial PDEs Computational methods for solving interfacial mass transfer and reaction models were reviewed by Glasscock and Rochelle (1989) and Van Swaaij and Versteeg (1992). In this work, the discretised equations are generated in a manner similar to the 5nite volume method (Patankar, 1980), i.e. by integrating Eq. (10) over an one-dimensional control volume. Di?usion terms are approximated with central di?erences. The grid involves two neighbouring cells at the previous time point, and the cell in question at the previous and new time point (implicit method). The non-linear term f is linearised according to Newton’s method, with numerical calculation of the partial derivatives by perturbating f around Ci at the previous time point. 6. Case study: Parallel–consecutive reactions in a draft-tube bubble column In the example that follows, the usage of this methodology will be demonstrated via a gas–liquid reaction system (absorption of CO2 into alkali solution) taking place in a draft tube bubble column. It must be stressed, however, that the particular CFD implementation described below does not constitute a part of the methodology; in principle any CFD model can be used in conjunction with the generalised compartmental model. 6.1. CFD model Before integrating with the reaction engineering model, the hydrodynamic predictions of the CFD simulation will be compared with the experimental data of Wachi, Jones, and Elson (1991), who reported gas hold-up in the riser and downcomer sections of a draft-tube bubble column measured using the piezometric method. Di?erent liquids, gas 6owrates and internal draft tubes were tested in that work, but the data presented here for comparison are for the air– water system and a 0:12 m diameter draft tube. More details of the experiments may be found in that work. Recent reviews of CFD models of gas–liquid 6ows such as those encountered in bubble columns have

3084

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089

gas hold-up in riser, %

10 riser, model riser, experiment downcomer, model downcomer, experiment

8 6 4 2 0 0

2

4

6 3

8

10

-4

gas flowrate, m /s (x10 ) Fig. 3. Simulated and experimental gas hold-up in the riser and downcomer (experimental results from Wachi et al., 1991).

been contributed by Jakobsen et al. (1997) and Joshi (2001). Among the main kinds of models that can be employed (volume tracking, Eulerian–Lagrangian, Eulerian–Eulerian) it is generally accepted that the Eulerian–Eulerian approach is best suited for resolving the main 6ow characteristics in large scale equipment within a reasonable computational time (for its derivation via averaging see, e.g. Ishii, 1975 or Drew, 1983). With respect to turbulence modelling, for engineering— scale gas—liquid 6ow computations so far only the k– model has been applied and tested, though it is unlikely to be the most suitable (Jakobsen et al., 1997). However, the fact that the methodology is not linked to a speci5c CFD model means that it will bene5t from any future developments in the 5eld. The drag force for a bubble swarm is the most important empirical correlation required to use this model. The approach of Schwarz and Turner (1988), used subsequently with comparative success in several studies of vertical bubbly 6ows such as Sokolichin and Eigenberger (1995), is adopted for the purposes of this study. FD = 5:0 × 104 G usl :

(16)

For the virtual tracer experiment, the feeding of reactant in a bubble column is emulated by introducing the passive scalar continuously throughout the column at a value proportional to the local interfacial area (Eq. (9)). Our simulations were conducted with the CFX 4.3 code (AEA Technology, Harwell, UK), which solves the Eulerian–Eulerian model with the Inter-phase Slip Algorithm—IPSA (Spalding, 1982). The gas in6ow was introduced as a mass source in the bottom cells. To account for the initial expansion of the column during the introduction of the gas in the bottom, the computational domain was extended to include certain cells above the free surface, initially containing only gas. Fig. 3 shows a comparison between experimental and simulated gas hold-up at 5ve di?erent 6owrates. The shortcomings of the Eulerian–Eulerian approach with respect to predicting 6ow details and radial phase dis-

tribution have been pointed out by many authors, such as Jakobsen (2000). In the present framework, however, CFD is employed for calculating averaged properties, and in this context the model seems to be adequate. In the riser, agreement between CFD and experiment is excellent and may be considered quantitative at the low gas 6owrates tested, which explore the whole range of homogeneous bubbly 6ow. A slight deviation at the highest 6owrates may indicate the beginning of the transition to the heterogeneous regime. Gas entrapment in the downcomer is, however, underpredicted, as it is a complex phenomenon related to size classi5cation of the bubbles, whose proper modelling is not possible at this level (i.e. single db ). The CFD model would be more accurate if applied to external loop reactors, where gas is totally separated at the top. In this case, however, holdup in the riser is more important, as the entrapped bubbles in the downcomer are soon depleted from the reactant component and play no important role in the reaction engineering. The compartmentalisation is to be carried out for the steady state. For this reason, the concentration di?erence between two points in the column (one in the riser and one in the downcomer) was monitored until a pseudo-steady state was reached. The virtual tracer experiment was then carried out in a manner resembling the feeding policy of bubble columns (Eq. (9)). 6.2. Network structure and mean ?owrates Fig. 4 shows the tracer concentration, 6ow 5eld and volume fraction pro5les (only half the column is plotted, as the simulation is axisymmetric). Note that the column head is occupied by gas, as it can be seen in the volume fraction contour plot, and therefore the tracer concentration in that region is of no importance. Fig. 5 shows radial gas holdup pro5les in the riser at di?erent axial points. By inspecting these results, we draw the following conclusions regarding the structure of the compartmental model: • The liquid phase 6ow 5eld is dominated by a circulation loop, characteristic of airlift reactors, due to the di?erential gas holdup between the riser and downcomer. This loop is the main macromixing mechanism and must be accounted for by the structure of the network. • The downcomer resembles a plug 6ow, with no observable gradient in the lateral direction. Therefore there is no need for lateral compartments in that region. • The volume fraction plot (Fig. 5) shows the gas to be concentrated in the middle of the riser and reduce towards the wall. Consequently the reactant is fed mostly throughout the central region of the riser, and disperses towards the walls mainly through turbulent mixing. The change in the direction of the 6ow, from upward in the riser to downward in the downcomer, suggests the use of two vertical layers of cells for the riser and downcomer, respectively. The fact that the gas holdup is concentrated in

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089

3085

Free surface

1.0000E+00 8.3334E-01 6.6667E-01 5.0001E-01 3.3334E-01 1.6668E-01 1.0385E-05

1.0705E+00 8.9211E-01 7.1368E-01 5.3526E-01 3.5684E-01 1.7842E-01 0.0000E+00

5.2040E-04 4.6600E-04 4.1159E-04 3.5719E-04 3.0279E-04 2.4839E-04 1.9398E-04

Gas insertion Fig. 4. CFD results, left to right: volume fraction, mean velocity pro5le (m/s) and dimensionless scalar concentration (resulting from the virtual tracer experiment). Axisymmetric simulation, half of the column shown, axis on the right.

the middle of the riser suggests that the riser be further splitted into 2 vertical layers, one for the gas-rich and one for the gas-poor region. Since the reactant is introduced to the bulk from wherever gas is present, this way we emulate the feeding policy of the reactor. Based on the above, a network comprising 10 horizontal and 3 vertical layers of cells is constructed. The cells are initially linked with 6owrates re6ecting only the convective mean 6ow, averaged along the cell boundaries. This network is depicted in Fig. 6 (without the dashed lines, which signify turbulent exchange 6owrates to be added later). 6.3. Turbulent exchange Figs. 7 and 8 show the axial concentration pro5les in the two riser layers and in the downcomer. It is evident that this network yields a rather poor description of the riser

concentrations, due to the lack of direct material exchange between the two vertical layers of riser compartments. The network is thus re5ned to include lateral turbulent exchange (Fig. 6 including the dashed lines). Agreement between the network and CFD is signi5cantly improved, as shown in Figs. 9 and 10. The discrepancy is now in the order of 1– 4%, and is deemed acceptable for the purpose of this study. 6.4. Superposition of chemical reactions Assume the following system of mass transfer accompanied by a system of parallel–consecutive reactions: CO2(aq) + OH−
HCO− 3 HCO− 3

−

+ OH


CO= 3

[i]; [ii]:

(17)

3086

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089

0.12

volume fraction

0.10

10

11

30

9

12

29

8

13

28

7

14

27

6

15

26

5

16

25

4

17

24

3

18

23

2

19

22

1

20

21

0.08 H = 2.2 m 0.06 0.04 0.02 0.00 0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.12

volume fraction

0.10 0.08 H = 1.5 m 0.06 0.04 0.02 0.00 0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.12

volume fraction

0.10 0.08

H = 0.8 m 0.06 0.04 0.02 0.00 0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.12

volume fraction

0.10 0.08 H = 0.3 m

0.06 0.04

Gas-rich

0.02 0.00 0.00

0.02

0.04

0.06

0.08

0.10

Gas-poor

0.12

distance from axis (m)

Fig. 5. Radial gas holdup pro5les in the column at four axial locations.

Kinetics and equilibrium parameters are taken from Astarita (1967) and are shown in Table 1. At pH=12–13 the equilibrium of this system is shifted towards CO= 3 ; below pH = 11, HCO− 3 exists in signi5cant amounts, and around pH = 10 it coexists with CO= 3 at equal concentrations, while at pH = 9 the molar fraction of CO= 3 is less than 0.1. As pH approaches the neutral range, CO2(aq) dominates. Figs. 11 and 12 show the time course of the concentration and pH pro5les, averaged over the whole column. In the beginning of the process almost all of the absorbed gas is = converted to CO= 3 . The increase in CO3 concentration soon reaches a peak however, and is subsequently being reduced, while HCO− 3 starts to increase. This event is re6ected in the pH curve by a change in slope, and can be explained by the shift of equilibrium of reaction (ii) in favour of HCO− 3 occurring around pH 11. At lower pH values the equilibrium

Riser

Downcomer

Fig. 6. Compartmental network with two vertical layers of riser cells, to emulate the feeding policy of a bubble column. The network is subsequently modi5ed to account for lateral turbulent exchange (the extra turbulent exchange 6owrates being represented by the dashed lines).

of reaction (i) shifts towards CO2 , and another change in the pH slope is observed. After this point, chemical reactions become negligible and the e?ect is that of physical absorption. Later absorption stops and pH becomes constant, as there is no driving force due to the high bulk concentration of dissolved CO2 . The shape of the pH curve is in qualitative agreement with the experiments of Fleischer, Becker, and Eigenberger (1996). In Fig. 13 we see a comparison of the total amount of CO2 absorbed in the end of the chemical process to the amount that would have been absorbed during physical absorption. It is much higher in the chemical process, partly but not only due to the enhancement of mass transfer by the reaction at the interface. In the chemical process, the bulk conversion

1.40

1.40

1.30

1.30

dimensionless concentration

dimensionless concentration

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089

1.20 1.10 1.00 0.90 0.80

Outer region, CFD

0.70

Central region, CFD Central region, network

0.60

Outer region, network

0.50 0.40 0

50

100

150

200

1.20 1.10 1.00 0.90 0.80 0.70 CFD Network

0.60 0.50 0.40 0

50

100

distance from bottom, cm

150

200

distance from bottom, cm

Fig. 7. Concentration pro5les in the riser, comparison of CFD with the network of Fig. 6.

Fig. 10. Concentration pro5les in the downcomer, comparison of CFD with the network of Fig. 6 with turbulent exchange.

Table 1 Parameter values

1.40 dimensionless concentration

3087

1.30 1.20 1.10 1.00 0.90 0.80 0.70 0.60

CFD Network

0.50 0.40 0

50

100 150 distance from bottom, cm

200

Fig. 8. Concentration pro5les in the downcomer, comparison of CFD with the network of Fig. 6.

Parameter (units)

Value

Ki (m3 mol−1 ) Kii (m3 mol−1 ) ki (mol−1 m3 s−1 ) kii (mol−1 m3 s−1 ) H (atm m3 mol−1 ) D (m2 s−1 ) db (m) usl (m s−1 ) (m) (s) P (atm)

6:1 × 104 5.88 5.0 1 × 104 0.0346 2:2 × 10−9 0.005 0.25 0:1 × 10−6 0.02 1

13 1.40

12

1.20 1.10

11

1.00 0.90 0.80

Outer region, CFD

0.70

Central region, CFD

pH

dimensionless concentration

1.30

Central region, network

0.60

10 9

Outer region, network

0.50 0.40 0

50

100

150

200

8

distance from bottom, cm

7 Fig. 9. Concentration pro5les in the riser, comparison of CFD with the network of Fig. 6 with turbulent exchange.

− of CO2 to CO= 3 and HCO3 keeps the concentration of dissolved CO2 in the bulk close to zero, thus maintaining the mass transfer driving force at high levels, until the time of the pH drop. On the other hand, the physical absorption process is sooner halted due to the accumulation of CO2 in the bulk. This example illustrates the need for rigorous dynamic modelling in the analysis of batch absorption processes.

0

10

20

30

40

50

60

70

time (min) Fig. 11. Evolution of pH.

7. Concluding remarks The main concept of this methodology is that mixing, both due to mean 6ow and turbulent exchange at the meso- and macro-scale, can be accounted for via averaging the 6ow

3088

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089

20

Notation C D FD ki

15

kii

30

concentration (mol/m3)

25

CO2 gmol/m3

10

Ki

CO3= gmol/m3 HCO3- gmol/m3

Kii

5 0 0

10

20

30

40

50

60

70

time (min)

Greek letters

Fig. 12. Evolution of chemical species’ concentrations.

 G L 

35 30 total CO2 absorbed (mol)

H P usl

concentration, mol m−3 di?usion coe:cient, m2 s−1 drag force, kg m s−2 kinetic constant of sub-reaction (i), mol−1 m3 s−1 kinetic constant of sub-reaction (ii), mol−1 m3 s−1 equilibrium constant of sub-reaction (i), m3 mol−1 equilibrium constant of sub-reaction (ii), m3 mol−1 Henry’s constant, atm m3 mol−1 pressure, atm slip velocity of bubbles, m s−1 penetration depth, m volume fraction of gas, volume fraction of liquid, density, kg m−3 contact time, s

25 20

References

15 Chemical absorption 10 Physical absorption 5 0 0

10

20

30

40

50

60

70

time (min) Fig. 13. Physical and chemical absorption.

and concentration 5elds resulting from a CFD 6ow simulation and a computational tracer experiment, and using these results to parametrise a generalised compartmental model. The approach introduces no new phenomenological concepts, and its accuracy depends on the accuracy of the CFD (whose main phenomenological aspects are the turbulence and multiphase 6ow models) and the 5ne-tuning between the CFD and compartmental model tracer simulations. The former will bene5t from any future advances; the later depends on the time and e?ort the process designer can afford to invest. Application to bubble column reactors was accomplished here by superimposing a phenomenological model of interfacial mass transfer and chemical reaction on the generalised compartmental model. The main concept, however, is more broad and may be applied to several other 5elds such as stirred tank reactors (in conjunction with phenomenological micromixing models), biochemical and precipitation reactors (in conjunction with population balance models).

Astarita, G. (1967). Mass transfer with chemical reaction. Amsterdam: Elsevier. Baldyga, J., Bourne, J. R., & Hearn, S. J. (1997). Interaction between chemical reactions and mixing on various scales. Chemical Engineering Science, 52, 457–466. Bauer, M., & Eigenberger, G. (1999). A concept for multi-scale modelling of bubble columns and loop reactors. Chemical Engineering Science, 54, 5109–5117. Bezzo, F., Macchietto, S., & Pantelides, C. C. (2000). A general framework for the integration of computational 6uid dynamics and process simulation. Computers and Chemical Engineering, 24, 653–658. Danckwerts, P. V. (1953). Continuous systems: Distribution of residence times. Chemical Engineering Science, 2, 1–13. Danckwerts, P. V. (1958). The e?ect of incomplete mixing on homogeneous reactions. Chemical Engineering Science, 8, 93–100. Delnoij, E., Kuipers, J. A. M., & van Swaaij, W. P. M. (1997). Computational 6uid dynamics applied to gas–liquid contactors. Chemical Engineering Science, 52, 3623–3638. Drew, D. A. (1983). Mathematical modelling of two-phase 6ow. Annual Review of Fluid Mechanics, 15, 261–291. Dudukovic, M. P., Larachi, F., & Mills, P. L. (1999). Multiphase reactors—revisited. Chemical Engineering Science, 54, 1975–1995. Fleischer, C., Becker, S., & Eigenberger, G. (1996). Detailed modelling of the chemisorption of CO2 into NaOH in a bubble column. Chemical Engineering Science, 51(10), 1715–1724. Glasscock, D. A., & Rochelle, G. T. (1989). Numerical simulation of theories for gas absorption with chemical reaction. A.I.Ch.E. Journal, 35(8), 1271–1281. Higbie, R. (1935). The rate of absorption of a pure gas into a still liquid during a short time of exposure. Transactions of A.I.Ch.E., 31, 365–389. Hindmarsh, A. C. (1983). Odepack, a systematized collection of ode solvers. In R. S. Stepleman, et al. (Eds.), Scienti@c Computing. Amsterdam: North-Holland, pp. 55 – 64.

S. Rigopoulos, A. Jones / Chemical Engineering Science 58 (2003) 3077 – 3089 HostomskWy, J., & Jones, A. G. (1995). A penetration model of the gas– liquid reactive precipitation of calcium carbonate crystals. Transactions of the Institution of Chemical Engineers, 73, 241–245. Ishii, M. (1975). Thermo-?uid dynamic theory of multi-phase ?ow. Eyrolles, Paris. Jakobsen, H. A. (2000). Phase distribution phenomena in two-phase bubble column reactors. Chemical Engineering Science, 56, 1049–1056. Jakobsen, H. A., Sannaes, B. H., Grevskott, S., & Svndsen, H. F. (1997). Modelling of vertical bubble-driven 6ows. Industrial and Engineering Chemistry Research, 36, 4052–4074. Joshi, J. B. (2001). Computational 6ow modelling and design of bubble column reactors. Chemical Engineering Science, 56, 5893–5933. Kuipers, J. A. M., & van Swaaij, W. P. M. (1997). Application of computational 6uid dynamics to chemical reaction engineering. Reviews in Chemical Engineering, 13, 1–118. Launder, B. E., & Spalding, D. B. (1972). Lectures in mathematical modelling of turbulence. New York: Academic Press. Mann, R., & Hackett, L. A. (1988). Fundamentals of gas–liquid mixing in a stirred vessel: An analysis using networks of backmixed zones. Proceedings of the sixth European conference on mixing, BHRA, Pavia (p. 321). Mecklenburgh, J. C., & Hartland, S. (1975). Theory of backmixing. New York: Wiley-Interscience. Montante, G., Lee, K. C., Brucato, A., & Yianneskis, M. (2001). Experiments and predictions of the transition of the 6ow pattern with impeller clearance in stirred tanks. Computers and Chemical Engineering, 25, 729–735. Nauman, E. B., & Bu?ham, B. A. (1983). Mixing in continuous ?ow systems. New York: Wiley. Patankar, S. V. (1980). Numerical heat transfer and ?uid ?ow. Washington, DC: Hemisphere Publishing. Pope, S. (2000). Turbulent ?ows. Cambridge: Cambridge University Press.

3089

Ranade, V. V. (1995). Computational 6uid dynamics for reactor engineering. Reviews in Chemical Engineering, 11(3), 229–289. Rigopoulos, S., & Jones, A. G. (2001). Dynamic modelling of a bubble column for particle formation via a gas–liquid reaction. Chemical Engineering Science, 56(21–22), 6177–6183. Romanainen, J. J. (1997). Numerical approach to modelling of dynamic bubble columns. Chemical Engineering and Processing, 36, 1–15. Schwarz, M. P., & Turner, W. J. (1988). Applicability of the standard k– turbulence model to gas-stirred baths. Applied Mathematical Modelling, 12, 273–279. Sokolichin, A., & Eigenberger, G. (1995). Gas–liquid 6ow in bubble columns and loop reactors: Part I. Detailed modelling and numerical simulation. Chemical Engineering Science, 49, 5735–5746. Toor, H. L., & Marchello, J. M. (1958). Film-penetration model for mass transfer and heat transfer. A.I.Ch.E. Journal, 4, 97–101. Torvik, R., & Svendsen, H. F. (1990). Modelling of slurry reactors. A fundamental approach. Chemical Engineering Science, 45, 2325–2332. Van den Akker, H. E. A. (1997). On the status and merits of computational 6uid dynamics. In W. A. Nienow (Ed.), Proceedings of the fourth international conference on bioreactor and bioprocess ?uid dynamics (pp. 407– 431). Van Elk, E. P., Borman, P. C., Kuipers, J. A. M., & Versteeg, G. F. (1999). Modelling of gas–liquid reactors: stability and dynamic behaviour of gas–liquid mass transfer accompanied by irreversible reaction. Chemical Engineering Science, 76, 223–237. Van Swaaij, W. P. M., & Versteeg, G. F. (1992). Mass transfer accompanied with complex chemical reactions in gas–liquid systems: An overview. Chemical Engineering Science, 47, 3181–3195. Wachi, S., Jones, A. G., & Elson, T. P. (1991). Flow dynamics in a draft-tube bubble column using various liquids. Chemical Engineering Science, 46, 657–663. Zauner, R., & Jones, A. G. (2000). Scale-up of continuous and batch precipitation processes. Industrial and Engineering Chemistry Research, 39(7), 2392–2403.