Modeling of Semibatch Agglomerative Gas−Liquid Precipitation of

Ind. Eng. Chem. Res. 2003, 42, 6567-6575

6567

Modeling of Semibatch Agglomerative Gas-Liquid Precipitation of CaCO3 in a Bubble Column Reactor Stelios Rigopoulos† and Alan Jones* Department of Chemical Engineering, University College London, Torrington Place, London WC1E 7JE, U.K.

The outcome of gas-liquid precipitation in industrial reactors, such as bubble columns, is determined by the interplay between multiphase fluid dynamics; gas-liquid reaction engineering; and crystallization mechanisms such as nucleation, growth, and agglomeration. In this work, a modeling approach that takes the above phenomena into account is proposed and applied to investigate the semibatch precipitation of CaCO3. The main elements of the approach are a dynamic population balance equation including nucleation, growth, and agglomeration, discretized with a finite element method; a phenomenological model of interfacial mass transfer and reaction based on penetration theory; and rigorous prediction of gas holdup via an EulerianEulerian multiphase CFD code. Experiments on CaCO3 precipitation via reaction of CO2 and Ca(OH)2 in a 21-L bubble column are conducted and simulated with the aid of the model. The evolution of the process is adequately reproduced, and qualitative comparisons with the product CSD are possible, but more fundamental work on precipitation kinetics in the gas-liquid hydrodynamic environment is required to obtain quantitative agreement. Introduction Gas-liquid precipitation is the formation of a solid crystalline product via chemical reactions involving a gas phase and a liquid phase. It is a process of considerable industrial importance, as it is involved in the production of a variety of speciality chemicals.1,2 The ability to quantify the phenomena that effectively determine the product crystal size distribution (CSD) and establish a predictive model based on design and operating parameters would be of great value for process design and scale-up purposes, as it would lead to improved product quality and enable tailoring of the product properties to new applications. Gas-liquid reactive precipitation was first studied in laboratory-scale reactors with a flat gas-liquid interface.3-5 Other laboratory-scale devices that have been employed are an agitated tank in which the gas was dispersed in the form of bubbles6,7 and a Couette-Taylor reactor.8 The latter was employed to study the morphology of precipitated CaCO3 by careful control of solution composition and shear. A proper analysis of industrial-scale equipment, however, requires consideration of hydrodynamics, and accurate prediction of the CSD can be accomplished by solution of the population balance equation. The coupling of the above approaches requires excessive computational resources, and it has been realized1 that a stagewise model should be employed rather than a differential one. The complexity is further increased by the fact that such processes are usually operated in batch or semibatch mode with time-varying supersaturation levels and precipitation rates, so that a timedependent model is required and numerical methods must be employed for its solution. Earlier work by the present authors9 employed a stagewise dynamic model for the simulation of a tall * To whom correspondence should be addressed. E-mail: [email protected]. † Present address: Department of Mechanical Engineering, Imperial College, Exhibition Road, South Kensington, London SW7 2BX, U.K.

bubble column carrying out precipitation of CaCO3, based on an implementation of penetration theory and CFD predictions. Mean crystal size was predicted with the method of moments, resulting in qualitative reproduction of the main trends, but several issues such as the rapid increase in mean size after the pH drop, possibly attributed to agglomeration, were left unaddressed. This was because both the experimental results and the model predictions reported only mean values of the CSD. Agglomeration seems to be pronounced during gas-liquid precipitation, as experiments7 have shown. It has a profound effect on the CSD and might determine the physical properties of the end product as well as its potential for separation; depending on the application, either its promotion or prevention might be desirable. The mathematical analysis of agglomeration is complex, however: it requires solution of the population balance, and no definitive kinetic model is available. In this work, a recently developed finite element method for solving the time-dependent volume-based population balance11 is coupled with a reaction engineering model, to simulate the evolution of the CSD while accounting for agglomeration in the kinetic model. In addition, new experiments on CaCO3 precipitation in a bubble column, during which the entire CSD was retrieved by laser light scattering, are reported. Thus, deeper insight into the gas-liquid precipitation process and the underlying phenomena can be obtained, while it is also shown that a model combining CFD, reaction engineering, and population balance can provide design and scale-up data in reasonable CPU time. Experimental Part The objective of this section is to describe the equipment and analytical methods employed in the experimental work and to obtain several clues regarding the precipitation mechanism from a preliminary inspection of the experimental results. These clues will prove to be of vital importance for configuring the kinetic model of CaCO3 precipitation. It must be stressed that the present work does not constitute an exhaustive study of CaCO3 precipitation kinetics. Such a study would

10.1021/ie020851a CCC: $25.00 © 2003 American Chemical Society Published on Web 11/05/2003

6568 Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003

Figure 2. Evolution of mean crystal size during a typical experiment. Figure 1. Experimental rig used for gas-liquid precipitation.

require the isolation of individual mechanisms,8,12 which is not possible in a bubble column, an apparatus whose simplicity of design gives rise to extraordinary complexity in the physical and chemical phenomena. Instead, our purpose is to explore whether the complex hydrodynamic, reaction engineering, and crystallization phenomena taking place in an industrial type of reactor, such as the bubble column, can be investigated with the aid of a mathematical model. Gas-Liquid Precipitation Rig. The experimental work was conducted using a bubble column reactor, equipped with an internal concentric draft tube to enhance circulation and mixing. The outer cylinder was 0.77 m high and had an internal diameter of 0.215 µm; the corresponding dimensions of the draft tube were 0.45 and 0.135 m, respectively. Gas (a mixture of CO2 and N2) was injected through a sparger ring situated at the bottom of the inner annulus (Figure 1). In each experiment, the column was filled with approximately 22 L of Ca(OH)2 solution [AnalaR grade Ca(OH)2deionized water]. Before initiation of the experiment, pure N2 gas was injected for some time to enhance the dissolution of the Ca(OH)2 by intense circulation and mixing. Different gas flow rates and initial solution compositions were explored. Because the aim was to observe the dynamic behavior of the process, during the course of each experiment, regular samples were taken and analyzed immediately for pH, calcium ion concentration ([Ca2+]), and CSD. The CSD analysis was carried out with a laser scattering instrument (Malvern Mastersizer S), equipped with lenses valid for measuring particles in the range 0.05900 µm. To gain more insight into the particle form, selected samples were dried and further analyzed by scanning electron microscopy (JEOL JSM-820). Time Course of Precipitation. Typical plots of the temporal evolution of pH, calcium ion concentration, mean crystal size, and CSD are shown in Figures 2, 5, and 7. The main features can be summarized as follows: Most of the precipitation occurs during a period of high pH (11-12), with calcium ion being consumed at a slightly decaying rate. This period is followed by a rapid fall in the pH, while the calcium ion concentration has reached a miminum. Subsequently, the calcium ion concentration increases again. After an induction period, the mean crystal size slowly increases. As soon as the pH drops, a rapid increase in crystal size occurs, eventually stabilizing around 5-6 µm.

The first CSD measurement (immediately after the induction period) is dominated by a single bell-shaped peak around 1 µm. Later, a new peak appears at the size range 0.1-1 µm. The bimodal distribution remains until the end. Product Morphology. Figure 3 a,b shows a largescale and a closer view, respectively, of the precipitate at an early stage of the process (soon after the end of the induction time), evidently comprising mostly single crystals. Occasional agglomerates are also present, and in Figure 3c is presented a close view at one of them. The initial agglomerates have sharp, clearly defined edges, and the single crystals that comprise them can be clearly identified. At a later stage, after the pH drop, the product is dominated by agglomerates (Figure 3d), confirming the increase in mean size. Meanwhile, the edges of the agglomerates become increasingly smoother, until reaching a stage where only agglomerates exist, exhibiting an almost amorphous shape (Figure 3e). At the same time, the solids concentration is apparently lower. Interpretation of Results. The absorption of CO2 into alkali solutions is accompanied by the formation of several ionic species, of which CO32- and HCO3predominate in the alkaline region. Figure 4 shows the equilibrium concentrations as functions of pH, calculated using the data of Astarita.13 Note that CO32-, which contributes to the precipitation process, is present at a significant proportion only at pH > 10. This implies that, after the occurrence of the pH drop in the process, CO32- is no longer present, and precipitation has virtually terminated. The increase in mean size after the pH drop is in accordance with previous findings,1,10 which were explained in the following way: During the high pH stage, particles experience slow growth, while agglomeration is being inhibited by the electrical charge on the particles’ surface. Once the pH has fallen, the mean size increases rapidly through agglomeration.1 This interpretation, however, leaves an important question unanswered: How does agglomeration occur after the depletion of supersaturation, which is necessary to create solid bridges and unify the loose aggregates? Furthermore, the experiments of Kotaki and Tsuge7 on the same system clearly evidenced agglomeration at high pH; our own experiments also detected agglomerates during early stages (Figure 3c). The evolution of crystal shape provides vital clues for the understanding of the process. The loss of the sharp edges in the agglomerates, coinciding with the reappearance of free calcium ion after the pH drop, suggests

Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 6569

Figure 3. (a) Large-scale view of the precipitate, shortly after the induction time. (b) Closer view of the single product crystals. (c) Close view of an agglomerate during an early stage of the process. (d) View of the precipitate during a later stage of the process. (e) Closer view of an agglomerate during a very late stage.

Figure 4. Equilibrium concentrations of the CO32--HCO3--CO2 system

that the crystals are redissolving. Edges exhibit more rapid dissolution, thus leading to the almost amorphous material of Figure 3e. The redissolution can be explained by the shift of the equilibrium from CO32- ion toward HCO3- ion. At high pH, when CO32- ion is favored, a balance is maintained between the migration of ions from the crystal lattice to the surrounding solution and deposition of ions in the lattice. After the pH drop, however, the CO32- ion that migrates is converted to HCO3- ion, and as a result, the balance is disrupted. The overall dissolution process can be described by the following reaction

CaCO3 + H2O T Ca

2+

-

-

+ HCO3 + OH

The rate of dissolution of a crystal is proportional to its surface area. Smaller crystals with a higher surface area-to-volume ratio disappear faster, leaving only the larger agglomerates in the system. This explains the fact that agglomerates dominate the precipitate after the pH drop, despite the fact that supersaturation has been depleted. The increase in the mean crystal size after the pH drop is therefore only an apparent growth. Bimodality in the CSD is usually explained by agglomeration or breakage. Although agglomerates were indeed detected in the precipitate through SEM, the fact that the second peak appears in the small size range indicates that agglomeration is not the cause. Disruption of the weak aggregates is another candidate, but

it must be ruled out by the fact that the new peak appears at sizes smaller than the size of the single crystals that comprise the agglomerates (∼1 µm). The crystals comprising the new peak must, therefore, be produced by a new wave of nucleation-growth. This could be the effect of secondary nucleation, which is proportional to the existing crystal mass.14,15 Size-independent growth exhibits a third-order polynomial increase in crystal mass with respect to time. It is reasonable, then, that, at some point, a sharp increase in crystal mass will rapidly generate a large amount of secondary nuclei, thus inducing a new wave of crystal growth. To confirm this speculation, an experiment was conducted during which the introduction of CO2 was halted just before the new peak was expected to appear. N2 continued to be introduced, maintaining the total gas flow rate constant, to ensure that the turbulent conditions would remain unchanged. After 30 min, the distribution remained unchanged (i.e., single peak), indicating that disruption or breakage of the existing crystals alone does not account for the new peak. As soon as CO2 was introduced again, the new peak appeared, thus demonstrating that it consists of freshly grown particles and does not occur without the presence of supersaturation. Process Simulation Population Balance. Most of the previous studies of gas-liquid precipitation have employed a nucleationgrowth population balance model, solved using the method of moments.5,9 Experiments, however, in both gas-liquid4 and liquid-liquid16,17 precipitation of CaCO3 have evidenced the presence of agglomeration and demonstrated its importance in determining the product CSD. A new method for solving the time-dependent population balance equation with combined nucleation, growth, agglomeration, and breakage was recently developed by the authors11 and will be employed in this study. It belongs to the class of finite element methods and emphasizes conservation of moments and stability. A brief description follows. The time-dependent population balance equation that describes the evolution of the CSD in a finite, spatially

6570 Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003

uniform domain, with particle volume as the “internal” coordinate, and including nucleation, growth, and agglomeration, is formulated as

∂n(v,t) nin(v,t) - n(v,t) ∂ ) - [G(v) n(v,t)] + ∂t τ ∂v 1 v B0δ(v - v0) + 0 βa(v - v′,v′) n(v - v′,t) n(v′,t) dv′ 2

∫

n(v,t)

∫0∞βa(v,v′) n(v′,t) dv′

(1)

where n(v,t) and nin(v,t) stand for the population density at the reactor and at the inlet, respectively; G(v) represents the volumetric growth rate; and B0, βa, and v0 stand for the nucleation rate, agglomeration kernel, and volume of the nuclei, respectively. The following initial and boundary conditions apply

n(v,0) ) n0(v) n(0,t) ) 0

(initial distribution)

(2)

(no crystals of zero size)

(3)

The solution is obtained by dividing the domain of the independent variable (particle volume) into a number of elements, over which the solution is approximated with first-order trial functions that are subsequently collocated at the element boundaries. The use of firstorder functions contributes to the stability of the method and allows it to capture the nearly discontinuous solutions that can arise from the combination of nucleation and growth in precipitation. Furthermore, the method emphasizes conservation of the moments during the evaluation of the aggregation source integral. The time-dependent distributed population balance equation is transformed into a system of ordinary differential equations describing the temporal evolution of the population density at the element boundaries, and the solution is propagated along time. Reaction Engineering. The concept of penetration theory13 is adopted here to treat simultaneous mass transfer and chemical reactions at the interface. Any reaction system can thus be incorporated without the need to resort to simplified assumptions with respect to the reaction regime (e.g., slow, instantaneous). The implementation is similar to our previous work:9 the interface and the bulk are considered as two separate dynamic reactors that operate independent of each other and interact at discrete time intervals. The diffusion and reaction of chemical components is described by the following system of PDEs (note that variables at the interface and bulk are denoted by superscripts I and B, respectively)

∂CIi ∂t

)D

∂2CIi ∂x

2

K

+

rk(CI1, CI2, ..., CIn) ∑ k)1

(4)

with initial and boundary conditions

dCBi dt

K

)

∑ rk(CB1 ,CB2 ,...,CBn ,n1,n2,...,nm) k)1

dnj ) fj(CB1 ,CB2 ,...,CBn ,n1,n2,...,nm) dt

t ) 0, x > 0 w x ) 0, t > 0 w

CIi

x ) δ, t > 0 w

dCIi )0 dx (nonvolatile species) (7)

(5)

C/i

) (volatile species) (6)

(8)

(9)

Here, fj(CBi ,nj) stands for the right-hand side of the ODEs into which the original population balance is transformed via the finite element discretization. Initial conditions are calculated from the mixing of bulk and interface at the end of the previous contact time. Note that the mass balances and the ODEs resulting from the population balance are coupled via the source terms. The coupling results because (a) precipitation kinetics are strongly dependent on supersaturation and (b) the rate of consumption of the precipitating species is the rate of volumetric crystal growth and must be calculated from the population balance. The solution of the interface PDEs is obtained numerically with an implicit iterative scheme, while the bulk ODEs are propagated with the LSODE set of subroutines18 implementing an Adams method. Reaction Kinetics. The reaction of CO2 into alkali solutions has been extensively studied.13,19 The first step to consider in a CO2 absorption system is the gas-liquid equilibrium, which determines the boundary condition for the volatile species

CO2(g) f CO2(aq)

CIi

)

CBi

Previous studies of the flat-interface reactor4,5 have accounted for precipitation occurring at the interface, because the low mass-transfer coefficient in that apparatus results in very long values of the contact time5 (2-16 s). In most industrial equipment, however, such as the bubble column, the contact time is several orders of magnitude lower (0.02 s for 5-mm bubbles moving at a slip velocity of 0.25 m/s, values typical for such equipment and assumed in this work), and the extent of precipitation phenomena at the interface is negligible compared to the bulk. If significant spatial variation in the concentration profiles in the reactor were observed (as in ref 9), a network of backmixed compartments would be required to describe the bulk. However, the bubble column used in the present work was found to be almost fully backmixed because of its short height and intense recirculation; moreover, the uncertainty in precipitation kinetics means that a too-detailed model would be wasted. However, nonideal mixing should certainly be considered in future work, possibly in conjunction with techniques such as CFD or network-of-zones, provided that satisfactory models of precipitation kinetics become available. In the present model, the bulk is considered to be fully mixed and described by the species’ mass balances and the ODEs for the nodal values of the population density

(i)

This step can be approached by an equilibrium equation such as Henry’s law. To calculate Henry’s constant in an ionic system, the effect of ionic strength must be taken into account19

log

H Hw

)-

∑i Iihi

(10)

Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 6571

hi ) h+ + h- + hg

(11)

Table 1. Parameter Values Used in the Simulations parameter (units)

The values for the contributions of cations (h+), anions (h-), gas (hg), and compounds (hi) are taken from ref 21. The kinetics of CO2 absorption into alkali solutions are determined by the conversion of CO2(aq) into HCO3-, a reaction that proceeds at a great, but finite rate. This is followed by an instantaneous ionic reaction

CO2(aq) + OH- h HCO3-

(ii)

HCO3- + OH- h CO32-

(iii)

The values of the kinetic and equilibrium constants are from Astarita.13 The kinetic constant of instantaneous reaction iii is approximated by a very high value. The precipitation of CaCO3 can be written concisely as

Ca

2-

+ CO3

h CaCO3(s)

d[CaCO3(s)] ) dt

FCaCO3

∫vv G(v) n(v) dv‚MWCaCO ∞

0

(12)

3

Note how the CSD enters this expression, in the form of the population density n(v). Thus, the rate of crystal mass production is inherently coupled with the population balance, and to estimate it, we need a complete kinetic model of precipitation that accounts for all crystallization mechanisms, namely, nucleation, crystal growth, and interparticle phenomena such as agglomeration and breakage. The kinetics of the precipitation mechanisms are, in turn, functions of the chemical species’ concentrations, as they depend on supersaturation. It is generally agreed that the driving force of supersaturation is best related to R+βxλs - 1, where R and β are the kinetic orders of the cation and anion in the ionic reaction, respectively, and λs is the saturation ratio. In an ionic solution, this is best expressed in terms of activities17,22

xλs - 1 )

x

RCa2+RCO32Ksp

-1

(13)

Activities are calculated according to the Debye-Hu¨ckel

6.1 × 104 5.88 5.0 1 × 104 0.0047 6.2 × 10-11 7.5 × 106 6 1.6 × 1014 0.1 1 1.2a 1.8b

a At 0.6 × 10-4 m3/s gas flow rate. b At 1.1 × 10-4 m3/s gas flow rate.

equation, with a modification proposed by Davies23

Ri ) γiCi

(

(iv)

This expression, however, obscures the underlying crystallization mechanisms that determine the rate. The rate of CaCO3(s) production is essentially that of volumetric crystal growth, but it is size-dependent even when the linear growth is size-independent (McCabe’s ∆L law). To obtain the rate of change for the whole crystal mass, we need to integrate the volumetric growth function over the whole range of crystal volumes

mol-1)

equilibrium constant of reaction i, Ki equilibrium constant of reaction ii, Kii (m3 mol-1) kinetic constant. of reaction i, ki (mol-1 m3 s-1) kinetic constant of reaction ii, kii (mol-1 m3 s-1) solubility product, Ksp (mol2 m-6) growth constant, kg (m/s) kn1 (eq 20) (s-1) kn2 (eq 20) kn3 (eq 20) (kg-1) kn4 (eq 20) kn5 (eq 20) agglomeration constant, ka (m3/s)

log(γi) ) -ADHzi2 2+

value (m3

(14)

)

xI - 0.3I 1 + xI

(15)

where ADH is the Debye-Huckel constant and I is the ionic strength, given in terms of the concentration and charge on each ionic species

I)

1

Cizi2 ∑ 2 i

(16)

Precipitation Kinetics. The following expression is usually employed to describe growth rate kinetics

Gl ) kg(λs - 1)2

(17)

Here, Gl is the linear growth rate, i.e., the increase in particle diameter or radius, and kg is a kinetic constant. The experiments of several authors17,24 seem to agree on the order of magnitude of the parameter kg; the data of the former17 were used in this work. See Table 1 for a summary of the parameter values used. Nucleation, on the other hand, is much more difficult to study because of the variety of mechanisms (homogeneous, heterogeneous, secondary) that might be responsible for it. Primary nucleation depends mainly on supersaturation, and a power law of the following form is commonly used to describe it

B0 ) kn1(λs - 1)kn2

(18)

Secondary nucleation is induced by the existing crystals14,15 and is a function of the crystal mass

B0 ) kn1(λs - 1)kn2Mckn3

(19)

Although several studies on calcium carbonate have appeared,25,26 their results cannot be extrapolated to equipment of different type and scale, where nucleation might occur via a different mechanism. Moreover, relatively few data have been reported on gas-liquid systems. Most important is the fact that, during a batch process, different nucleation mechanisms might prevail during different periods. In our case, the experiments indicate that, in the beginning, high supersaturation levels induce primary nucleation, but later, secondary nucleation gives rise to a new wave of crystal growth.

6572 Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003

The two mechanisms are fundamentally different and can operate independent of each other. To account for these facts, our overall nucleation model consists of the sum of the two models, primary and secondary

B0 ) kn1(λs - 1)kn2 + kns(λs - 1)kn3Mckn4

(20)

Finally, although the prevailing polymorph appears to be calcite, it must be acknowledged that all three polymorphs of CaCO3 are likely to appear during the precipitation process. The derivation of a kinetic model accounting for their simultaneous presence is an exceedingly difficult task, however, and will not be attempted here; one is referred to ref 22 for this purpose. Agglomeration of crystals is a two-step process. The first stage of agglomeration, i.e., the formation of flocculates through collisions and interparticle attraction, is akin to similar phenomena occurring in colloids and aerosols. Thus, the starting point for the formulation of a mathematical description of agglomeration are the theories developed in those fields, most importantly the pioneering work of Smoluchowski.27 The second step is the growth of crystalline material between the clusters at spots called cementing sites. This step differentiates crystal agglomeration from other clustering processes, and few mechanistic models of it have appeared so far. Certain studies28,29 have approached the problem by considering crystal growth at a cementing cite between a pair of spherical particles. Although these approaches have yielded valuable insights, the formulation of a definitive model of the agglomeration process as a whole is still faced with severe obstacles: (1) Particles of various states (flocculates, aggregates of various ages, firm agglomerates) are likely to be present in a crystal agglomerating system. CSD measurements, however, provide no information on these states. (2) The interparticle phenomena that modify the CSD are also various. Unsuccessful collisions, aggregate disruption, and the breakage of firm agglomerates are all possible, depending on the conditions, and it is very difficult to distingush between these events. (3) Current mechanistic models tend to describe aggregation as a two-particle event; in reality, however, it might occur through multiparticle clusters with bridges of various ages and strengths. (4) The size-independent kernel seems to be in good agreement with most agglomeration experiments in stirred tanks and is used extensively. This agreement, however, is likely to be due to mutual cancellation of competing mechanisms, and we do not know whether it would be maintained in an arbitrary environment. As a result, agglomeration kinetics are heavily systemdependent and difficult to scale-up. The dependency of the agglomeration kernel on hydrodynamics has been investigated mainly in stirred tanks,29 whose flow fields are very different from those of bubble columns. For this reason, no explicit dependency on hydrodynamics will be employed here; instead, the agglomeration constant will be adjusted to fit experiments with different gas flow rates. Among the most common theoretical agglomeration kernels, we found the fluid shear kernel (orthokinetic aggregation) to provide the best fit of our experiments. This is in good aggreement with the theory, but in contrast with most of the experimental data on stirred tanks, which showed the size-independent kernel to be more suitable.29 It must be borne in mind, however,

Figure 5. Experimental and simulated evolution of pH and calcium ion concentration: (a) [Ca(OH)2] ) 3 mol/m3, CO2/N2 ) 0.000 01:0.0001 m3/s; (b) [Ca(OH)2] ) 4 mol/m3, CO2/N2 ) 0.000 01: 0.0001 m3/s; (c) [Ca(OH)2] ) 4 mol/m3, CO2/N2 ) 0.000 01:0.000 05 m3/s.

that, when agglomeration is the only interparticle event in the population balance, the results refer to the net effect of an agglomeration-disruption equilibrium. Stirred tanks exhibit higher shear forces, and therefore, disruption of the weak, not fully cemented aggregates (especially the larger ones) is likely to occur more often, yielding an apparent size-independent agglomeration. On the other hand, the more gentle agitation in the bubble column might be the reason the effect of the sizedependent kernel becomes apparent. MSMPR experiments29 have shown agglomeration to be roughly proportional to growth. The underlying cause of both phenomena is supersaturation, and during a batch experiment, supersaturation varies. It seems appropriate, then, to correlate agglomeration to supersaturation directly. Given the second-order dependency of growth on supersaturation, we can infer a similar dependency to the agglomeration model. Our formulation is thus

βa(v′,v - v′) ) ka(λs - 1)2[(v′)1/3 + (v - v′)1/3]3 (21) Apart from the growth kinetics, the remaining parameters must be adjusted to fit the experiments. This

Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 6573

Figure 6. Simulation of the evolution of supersaturation in a typical experiment.

is especially difficult in a batch experiment, where different mechanisms can prevail at different time

periods and each mechanism is a nonlinear function of supersaturation. In our study, the determination of the parameters was guided by (a) the mean size of the crystals immediately after the induction period, which is determined by the balance between primary nucleation and growth; (b) the evolution of the Ca2+ concentration; and (c) the appearance of the second peak at the low size range, which indicates the effect of secondary nucleation. Hydrodynamics. The most important hydrodynamic parameter in gas-liquid precipitation is the gas holdup, as it determines the rates of the chemical phenomena. As in our previous work,9 the Eulerian-Eulerian multiphase CFD model is used to obtain the gas holdup. This approach is the best choice for modeling large-scale equipment, as it does not raise extreme computational demands and yields sufficiently accurate results with respect to overall (averaged) properties, although it is less successful in reproducing more fine details such as

Figure 7. Experimental and simulated CSDs, [Ca(OH)2] ) 3 mol/m3, CO2/N2 ) 0.000 01:0.0001 m3/s.

Figure 8. Experimental and simulated CSDs, [Ca(OH)2] ) 4 mol/m3, CO2/N2 ) 0.000 01:0.0001 m3/s.

6574 Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003

Figure 9. Experimental and simulated CSDs, [Ca(OH)2] ) 4 mol/m3, CO2/N2 ) 0.000 01:0.000 05 m3/s.

the radial phase distribution. Details on its theoretical derivation can be found in refs 30 and 31, and its application to vertical dispersed flows is discussed in ref 32. Drag force was specified according to the approach of Scwarz and Turner.33 Turbulence in the liquid phase was implemented with the aid of the k- model.33 The implementation was carried out with the CFX 4.3 code (AEA Technology, Harwell, UK). Simulation Results. Figure 5 shows a comparison of the simulated and experimental results for the evolution of pH and calcium concentration during the precipitation experiments. Overall, the time course of the process is well simulated. The measured pH curve seems to stabilize at a higher pH (7-8) than predicted, which falls to pH 5-6. This can be attributed to the kinetic model of the CO2-OH- process that included only those subreactions that are important in the high pH regime. At near-neutral pH, other ionic species apart from the ones considered (CO32-, HCO3-) might be present. Moreover, the pH meter was calibrated for the alkaline region. The Ca2+ plot is indicative of the precipitation phenomena. During the induction time, i.e., the time during which the nuclei are formed and grow to observable size, consumption of substrate is minimal, and the Ca2+ concentration is constant. Subsequently, a rapid drop in Ca2+ concentration indicates a high growth rate. This is due to very high supersaturation levels, generated by the accumulation of CO32- during the induction period. Rapid growth implies rapid consumption, and supersaturation soon falls to low levels, resulting in a lower growth rate. This rate is almost constant until the pH drop, indicating a balance between consumption and generation of supersaturation. The redissolution of the crystals, resulting in the reappearance of free calcium ion, is not included in the model. The above phenomena are reflected in the supersaturation plot, shown in Figure 6: a rapid increase during the induction time is followed by a rapid drop. Then, supersaturation is constant, until dropping to zero upon the depletion of calcium ion. Figures 7-9 show comparisons of the experimental CSD measurements and the simulation results. In most

cases, the main trends are reproduced, although agreement is not quantitative. The first CSD measurement usually exhibits the highest deviation from the modeling results, partly because of experimental error. This measurement was usually taken immediately after the induction time, when the solution was still too dilute and the obscuration of the Mastersizer S was below the limit required for a reliable measurement; large foreign particles would easily shift the CSD toward greater sizes. This happened in Figures 7 and 8, where the first measurement showed a much higher proportion of large particles, which is clearly unrealistic. The cumulative undersize distribution reveals that the deviation occurs mostly at the largest sizes, sometimes larger even than 20 µm, but such crystals were seldom observed by the SEM, which supports the possibility of the CSD being biased by foreign particles. The intermediate and final measurements, however, are in much better agreement with the predictions. Most important, the combination of secondary nucleation and agglomeration succeeds in reproducing the bimodality in the CSD. The lack of quantitative agreement can be attributed to the uncertainty in the precipitation kinetics, as the reaction engineering model was validated by the concentration predictions. This was to be expected, as the precipitation mechanisms have not been adequately studied in the gas-liquid environment. The modeling of agglomeration seems to pose the most difficult obstacles, because of the lack of a definitive mechanistic model of the cementation step. This problem is especially pronounced in this study for two reasons: (a) in a batch process, the agglomeration kinetics are timevarying, depending on supersaturation, and (b) the hydrodynamic environment that determines the aggregation-disruption equilibrium (or aggregation efficiency) in gas-liquid flows is fundamentally different from that in stirred tanks. The agglomeration model used in this study takes supersaturation into account and thus achieves reasonable predictions of the extent and time scale of the phenomena occurring during the semibatch process. The sustaining or disruption of the weak aggregates is not mechanistically incorporated,

Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 6575

however, and this might be the reason the exact shape of the CSD is not accurately reproduced. Conclusions In this work, population balance, reaction kinetics, and hydrodynamic principles were integrated into a modeling approach for the gas-liquid precipitation process. By accounting for all of the above considerations, the model aims to provide tools for the analysis and scale-up of industrial-class equipment. The approach was subsequently used to simulate the precipitation of CaCO3 by CO2 absorption into lime, and experiments on a laboratory-scale rig were carried out for comparison. The main findings can be summarized as follows: (1) The increase in mean size after the pH drop is due to the disappearance of the smaller crystals by dissolution. (2) Agglomeration starts to take place at relatively high pH and proceeds to a considerable extent, possibly because the aggregates are less frequently disrupted than in stirred tanks. For the same reason, the shear (orthokinetic) kernel succeeded in reproducing the net result of the aggregation-disruption processes better than the size-independent kernel. (3) Bimodality in the CSD was shown, by means of both experiment and model, to be the result of a new wave of nucleation-growth induced by the existing crystals (secondary nucleation). The conjunction of penetration theory and CFD predictions of the gas holdup seems to yield an adequate description of the reactor performance, as the prediction of the evolution of the species concentrations shows. Therefore, the discrepancies in the CSD predictions can be attributed mainly to the uncertainty in precipitation kinetics; the main trends, however, could be explained and reproduced. Systematic studies of the kinetics, as well as improved mechanistic models of agglomeration (especially in relation to hydrodynamics) and nonideal mixing, should be considered in future work to approach quantitative agreement. Literature Cited (1) Wachi, S.; Jones, A. G. Aspects of gas-liquid reaction systems with precipitate solid formation. Rev. Chem. Eng. 1995, 11 (1), 1-51. (2) Jones, A. G. Crystallization Process Systems; ButterworthHeinemann: Woburn, MA, 2001. (3) Yagi, H. Kinetics of solid production accompanying gasliquid reaction. In Proceedings of World Congress III of Chemical Engineers; Institution of Chemical Engineers: London, 1986; Vol. 4, pp 20-23. (4) Wachi, S.; Jones, A. G. Mass transfer with chemical reaction and precipitation. Chem. Eng. Sci. 1991, 46 (4), 1027-1033. (5) Hostomsky, J.; Jones, A. G. A penetration model of the gasliquid reactive precipitation of calcium carbonate crystals. Trans Inst. Chem. Eng. 1995, 241-245. (6) Yagi, H.; Iwazawa, A.; Sonobe, R.; Matsubara, T.; Hikita, H. Crystallization of calcium carbonate accompanying chemical absorption. Ind. Eng. Chem. Fundam. 1984, 23, 153-156. (7) Kotaki, Y.; Tsuge, H. Reactive crystallization of calcium carbonate by gas-liquid and liquid-liquid reactions. Can. J. Chem. Eng. 1990, 68, 435-442. (8) Jung, W. M.; Kang, S. H.; Kim, W.; Chang, K. C. Particle morphology of calcium carbonate precipitated by gas-liquid reaction in a Couette-Taylor reactor. Chem. Eng. Sci. 2000, 55, 733-747. (9) Rigopoulos, S.; Jones, A. G. Dynamic modelling of a bubble column for particle formation via a gas-liquid reaction. Chem. Eng. Sci. 2001, 56 (21-22), 6177-6183. (10) Jones, A. G.; Wachi, S.; Delannoy, C. C. Precipitation of calcium carbonate in a fluidised bed reactor. In Proceedings of the

Seventh Engineering Foundation Conference on Fluidization; Potter, O. E., Nicklin, D. J., Eds.; Engineering Foundation: New York, 1992; pp 407-414. (11) Rigopoulos, S.; Jones, A. G. Finite element scheme for solution of the dynamic population balance. AIChE J. 2003, 49 (5), 1127-1139. (12) Mumtaz, H. S.; Hounslow, M. J Aggregation during precipitation from solution: An experimental investigation using Poiseuille flow. Chem. Eng. Sci. 2000, 55, 5671-5681. (13) Astarita, G. Mass Transfer with Chemical Reaction; Elsevier Publishing: Amsterdam, 1967. (14) Botsaris, G. D. Secondary nucleation: A review. In 6th Symposium of Industrial Catallisation; Mullin, J. W., Ed.; Plenum Press: New York, 1976; pp 3-22. (15) Garside, J.; Davey, R. J. Secondary contact nucleation: Kinetics, growth and scale-up. Chem. Eng. Commun. 1980, 4, 393424. (16) Tai, C. Y.; Chen, P.-C. Nucleation, agglomeration and crystal morphology of calcium carbonate. AIChE J. 1995, 41 (1), 68-77. (17) Collier, A. P.; Hounslow, M. J. Growth and aggregation rates for calcite and calcium oxalate monohydrate. AIChE J. 1999, 45 (11), 2298-2305. (18) Hindmarsh, A. C. Odepack, a systematized collection of ode solvers. In Scientific Computing; Stepleman, R. S., Carver, M., Peskin, R., Ames, W. F., Vichnevetsky, R., Eds.; NorthHolland: Amsterdam, 1983; pp 55-64. (19) Danckwerts, P. V.; Sharma, M. M. The absorption of carbon dioxide into solutions of alkalis and amines (with some notes on hydrogen sulphide and carbonyl sulphide). Chem. Eng. 1966, 44 (3), 244-280. (20) Juvekar, V. A.; Sharma, M. M. Absorption of CO2 in a suspension of lime. Chem. Eng. Sci. 1973, 28, 825-837. (21) Pohorecki, R.; Moniuk, W. Kinetics of reaction between carbon dioxide and hydroxyl ions in aqueous electrolyte solutions. Chem. Eng. Sci. 1988, 43, 1677-1684. (22) Chakraborty, D.; Bhatia, S. K. Formation and aggregation of polymorphs in continuous precipitation. 2. Kinetics of CaCO3 precipitation. Ind. Eng. Chem. Res. 1996, 35, 1995-2006. (23) Davies, C. W. Ion Association; Butterworths: London, 1962. (24) Nancollas, G. H.; Reddy, M. M. The crystallisation of calcium carbonate: II. Calcute growth mechanism. J. Colloid Interface Sci. 1971, 37, 824. (25) Pakter, A. The precipitation of sparingly soluble alkalineearth metal and lead salts: Nucleation and growth order during the induction period. Inorg. Phys. Theor. 1968, 859, 862. (26) Swinney, L. D.; Stevens, J. D.; Peters, R. W. Calcium carbonate crystallisation kinetics. Ind. Eng. Chem. Fundam. 1982, 21, 31-36. (27) Smoluchowski, M. V. Mathematical theory of the kinetics of coagulation of colloidal systems. Z. Phys. Chem. 1917, 92, 129168. (28) David, R.; Marchal, P.; Klein, J.-P.; Villermauxx, J. Crystallization and precipitation engineering III. A discrete formulation of the agglomeration rate of crystals in a crystallization process. Chem. Eng. Sci. 1991, 205-213. (29) Hounslow, M. J.; Mumtaz, H. S.; Collier, A. P.; Barrick, J. P.; Bramley, A. S. A micro-mechanical model for the rate of aggregation during precipitation from solution. Chem. Eng. Sci. 2000, 2543-2552. (30) Ishii, M. Thermo-Fluid Dynamic Theory of multiphase Flow; Eyrolles: Paris, 1975. (31) Drew, D. A. Mathematical modelling of two-phase flow. Annu. Rev. Fluid Mech. 1983, 15, 261-291. (32) Jakobsen, H. A.; Sannaes, B. H.; Grevskott, S.; Svndsen, H. F. Modelling of vertical bubble-driven flows. Ind. Eng. Chem. Res. 1997, 36, 4052-4074. (33) Schwarz, M. P.; Turner, W. J. Applicability of the standard k- turbulence model to gas-stirred baths. Appl. Math. modell. 1988, 12, 273-279. (34) Launder, B. E.; Spalding, D. B. Lectures in Mathematical Modelling of Turbulence; Academic Press: New York, 1972.

Received for review October 28, 2002 Revised manuscript received September 16, 2003 Accepted September 17, 2003 IE020851A