5342
Ind. Eng. Chem. Res. 2005, 44, 5342-5352
Interaction between Mixing, Chemical Reactions, and Precipitation Jerzy Bałdyga* and Łukasz Makowski Faculty of Chemical and Process Engineering, Warsaw University of Technology, Warsaw, Poland
Wojciech Orciuch† Department of Chemical Reaction Engineering, Chalmers University of Technology, Gothenburg, Sweden
The way in which reagents are mixed can have a large influence on the product distribution of chemical reactions and the size distribution of particles of the solid product. To model effects of mixing on various scales on the course of chemical reactions and precipitation processes, a nonequilibrium multiple-time-scale mixing model and a beta distribution of the mixture fraction are applied in combination with a simple conditional moment closure based on linear interpolation of local instantaneous reactant concentration values. The mixing model is linked to CFD (standard k- model) and the model predictions are compared with experimental data for fast parallel chemical reactions and barium sulfate precipitation both carried out in the single-feed semibatch stirred-tank reactor. 1. Introduction The way in which reagents are contacted and mixed can have a large influence on the product distribution of fast complex chemical reactions and on the size distribution and morphology of particles of solid products. These effects are of high practical importance because many chemical reactions leading to desirable intermediate and end products (such as pharmaceutical intermediates, agrochemicals, and many other fine chemicals) are accompanied by side reactions producing undesired byproducts, which decreases reaction yield and complicates product separation. The simplest and straightforward way to improve selectivity is by enhancing competition between reactions, for example, by addition of a homogeneous catalyst that can increase the rate of the desired reaction. However, there is a natural limit for such procedure. Namely, when the desired reaction becomes very fast, its rate becomes controlled by mixing, rather than the reaction kinetics. The final selectivity results then from competition between the desired, mixing controlled chemical reaction and other, slower reactions that are often also affected by mixing. This shows that controlling of selectivity of complex chemical reactions requires knowledge of chemical reaction kinetics and thermodynamics as well as detailed information about bringing reagents together on the molecular scale, which in the case of homogeneous chemical reactions means mode of operation, method of feeding, and last but not least mixing of reagents on the molecular scale. The feature of sensitivity of chemical reactions to mixing allows using them as test systems to characterize mixing and to evaluate the quantitative predictions of models claiming to describe the coupling between mixing and reaction. Namely, to predict, control, and optimize mixing effects on chemical reactions, one needs to apply so called micromixing models often in combination with CFD. To * To whom correspondence should be addressed. Tel.: +48 22 660 6376. Fax: +48 22 825 6037. E-mail: baldyga@ ichip.pw.edu.pl. † Present address: Faculty of Chemical and Process Engineering, Warsaw University of Technology, Warsaw, Poland.
evaluate predictions of models describing reactive mixing, one can employ specially designed mixing sensitive test reactions.1 Mixing is a complex process that can be characterized by several time scales characterizing various stages of distributing one material in the other followed by homogenization of the mixture as discussed in Section 2 of this paper. The test reactions should be fast relative to mixing so that the time constants of chemical reactions are either smaller or of comparable magnitude as the time constants for mixing. The models describing the stages of mixing faster than the chosen test chemical reactions cannot be tested by using such reactions. Fast multiple reactions are the best to be employed as test reactions because in their case the product distribution stores the mixing history in addition to observed effects of mixing on spatial distribution of reactants (extend of the reaction zone) or variation of concentration in time. When the test reactions are designed, one needs to know exact kinetics, accurate quantitative analysis of reactants and products should be possible, and the reagents, products, etc. should be screened for toxicity, disposability, and fire and explosion hazards.1 Good examples of such reactions are as follows: (1) diazo coupling between 1-naphthol and diazotized sulfanilic acid,2 (2) simultaneous diazo coupling between 1- and 2-naphthols and diazotized sulfanilic acid,3 (3) competitive neutralization of hydrochloric acid and alkaline hydrolysis of monochloroacetate esters,4-7 (4) iodate/iodine reaction with neutralization,8 and (5) acetal hydrolysis with simultaneous neutralization.9 In this work we use the parallel reaction system that includes competitive neutralization of hydrochloric acid and alkaline hydrolysis of the ethyl chloroacetate: k1
NaOH + HCl 98 NaCl + H2O (A) (B) (P) k2
(1)
NaOH + CH2ClCOOC2H5 98 (A) (C) CH2ClCOONa + C2H5OH (2) (S)
10.1021/ie049165x CCC: $30.25 © 2005 American Chemical Society Published on Web 01/25/2005
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5343
This system of test reactions was used before to characterize mixing in the tubular reactor,10 the singlefeed5,7 and double-feed11 semibatch stirred-tank reactors, and the continuous-flow stirred-tank reactor.12 In this work the model predictions will be compared with new experimental data for the semibatch single-feed stirred-tank reactor. However, the main problem considered in this work is precipitation of sparingly soluble material C (barium sulfate) from two aqueous ionic solutions A and B (sodium sulfate and barium chloride) in the semibatch stirred-tank reactor.
An+ + Bn- f CV
(3)
Precipitation of sparingly soluble materials is an important unit operation in chemical engineering practice that is employed in the production of bulk and fine chemicals. The process of precipitation involves a mixing controlled chemical reaction and subsequent crystallization of the product, which includes nucleation and growth of particles and their possible agglomeration. The solid product has often a wide crystal size distribution (CSD), which determines its quality. Another factor affecting the product quality is the morphology of precipitating particles. Both the CSD and morphology can strongly depend on the method of contacting reactants and mixing conditions during the process. This results from the fact that the elementary subprocesses forming the overall precipitation process including chemical reaction, nucleation, and growth of crystals are usually very fast, so that mixing can affect their course. We will come back to this problem in Section 2 after presentation of the turbulent mixer model and defining of the time constants for precipitation subprocesses. Here we can just state that chemical reaction, nucleation, and crystal growth are essentially molecular level processes, and mixing on the molecular scale (micromixing) affects directly their course. However, other mixing subprocesses can indirectly affect or even control the course of micromixing. Mixing affects the generation of supersaturation of the product to be crystallized; in our case the instantaneous reaction is simply controlled by mixing, which obviously increases supersaturation. Mixing, however, also influences the redistribution of supersaturation in a system, dilutes the supersaturated solution with less supersaturated or unsaturated fluids, and allows contact of particles formed earlier with regions of high supersaturation where the particles grow, which obviously decreases supersaturation. Increasing of mixing rate also decreases the time period when the supersaturation in fluid elements is high. The nucleation and growth steps of crystallization usually differ significantly in sensitivity to supersaturation, which explains effects of supersaturation distribution on the particle size. As mixing and precipitation both affect this distribution, one observes both decrease and increase of the particle size with increasing mixing rate. The interpretation of mixing effects on precipitation through the modeling of the supersaturation distribution was studied in the subject literature by modeling of the supersaturation distribution in the limits of micromixing intensity,13 with the use of mechanistic micromixing models,14-18 and by applying CFD based methods usually in combination with the closure schemes.19-21 In this work the process of precipitation in the single-feed semibatch stirred-tank reactor is modeled by using CFD in combination with the simple conditional moment closure scheme based on linear
interpolation of the values of the local instantaneous reactant concentration between the values characterizing instantaneous and infinitely slow precipitation. The same closure method was used previously in our group to interpret the precipitation process in the tubular reactor,21 whereas precipitation in the singlefeed, semibatch, stirred-tank reactor was modeled using the mechanistic model of micromixing.15 Mixing effects on the precipitation carried out in the tubular reactor22 and in the Taylor-Couette-flow reactor23 were modeled recently using CFD and the multi-zone finite-mode PDF models. The main problem in modeling the course of the precipitation process carried out in the semibatch stirred tank using CFD is the long feed time. Difficulties result from the fact that it is necessary to follow the process using CFD during long feed times, many orders of magnitude longer than the time period of revolution of the impeller. This problem was solved in our previous work for the case of parallel reactions carried out in the semibatch stirred-tank reactors7,11 by discretizing the feed into the finite number of portions and solving for each portion the differential balance equations for mass, momentum, and species until the quasi-steady-state solution was obtained. This means that the semibatch feeding was simulated using a sequence of the quasisteady-state solutions. In the case of precipitation the problem is, however, more difficult than in the case of parallel chemical reactions. In the case of chemical reactions with one of them instantaneous it was possible to distinguish in the reactor the chemically active zone close to the feeding point and chemically passive bulk, whereas during most of the process one needs to follow precipitation progress in the bulk as well. 2. Multiple-Time-Scale Turbulent Mixer Model and Time-Scale Analysis The local instantaneous value of the concentration variance
b,t) ) 〈[c(x b,t) - 〈c(x b,t)〉]2〉 σS2(x
(4)
is a measure of the local instantaneous departure of concentration from its local instantaneous mean value that can be interpreted as a departure from locally perfect mixing, whereas dissipation of concentration variance can be interpreted as mixing on the molecular scale. More generally, the mixing process can be interpreted as convection, turbulent diffusion, as well as creation and dissipation of the concentration variance of a passive scalar. The equation governing the transport of the concentration variance reads
The term of generation of scalar fluctuations by gradients in 〈c〉 and local fluctuations of concentration and velocity (term I on the r.h.s) and the terms of diffusive transport expressed by term III (turbulent
5344
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
diffusion term IIIa and the molecular diffusion term IIIb) are modeled in what follows using the standard k- model, whereas the dissipation term is modelled employing the concept of the multiple-time-scale turbulent mixer.24 To this end we use the mixture fraction f, which represents a normalized concentration of a nonreacting passive scalar of the local instantaneous concentration c0A and feed concentration cA0. Mixture fraction has a meaning of the mass fraction of the fluid of the feed concentration cA0, and in the case of fluids of constant density it is equivalent to the volume fraction.
c0A f) cA0
(6)
The process of mixing is interpreted now as dissipation of the concentration variance
σS2 ) 〈(f ′)2〉
(7)
that can be modeled using the balance equation
[
]
∂σS2 ∂σS2 ∂σS2 ∂ ) (D + DT) + + 〈uj〉 ∂t ∂xj ∂xj ∂xj 2DT
( ) ∂〈f〉 ∂xj
2
- RD (8)
where RD is the rate of dissipation of σS2 and the turbulent diffusion coefficient DT is calculated from the k- model.
DT )
νT k2 ) 0.1 ScT 〈〉
(9)
In the equilibrium models of scalar dissipation it is usually assumed that the rate of dissipation of the concentration variance RD is proportional to the concentration variance σS2 with τm being the time constant for mixing.
RD )
σS2 τm
(10)
Of course molecular mixing is observed only at small length scales in the viscous-diffusive subrange of the spectrum of concentration fluctuations that can be represented by the Batchelor length microscale ηB and wavenumbers comparable with the Batchelor wavenumber kB ) 1/ηB. -1
kB
( )
D2ν ) ηB ) 〈〉
1/4
-1/2
) ηSc
(11)
where η represents the Kolmogorov length microscale,
η)
() ν3 〈〉
1/4
(12)
Equations 11 and 12 show that dissipation of the concentration variance occurs at scales much smaller than the Kolmogorov length microscale where the local flow of the fluid has a laminar or even creeping flow character. Mixing in the viscous-diffusive subrange represents the final stage of decay of concentration
fluctuations and can be strongly affected by mixing at larger scales in the viscous-convective and inertialconvective subranges. In the viscous-diffusive subrange the concentration fluctuations become so small that molecular diffusion can compete with deformation in affecting the concentration spectrum. Note, however, that ηB is the length scale which yields the same diffusion time as the characteristic deformation time τK that is equivalent to the Kolmogorov time microscale.
τK )
( ) ν 〈〉
1/2
)
ηB2 D
(13)
This observation results in fact from the more general feature of the viscous-convective and viscous-diffusive subranges. Namely, in the viscous-convective subrange at high Schmidt numbers the characteristic time scale of a scalar spot of a tracer is independent of scale or wavenumber, being determined by the time scale associated with the larger scale (of the order of η) straining motion characterized by τK. The statistical equilibrium means the quasi-steady behavior of the spectrum and then it can be assumed that eq 10 with the complete variance σS2 can characterize well the mixing rate in the viscous-diffusive subrange because the time constant τm should contain then relations between all the subranges involved. When, however, can the scalar variance spectrum be assumed quasi-steady? This would be possible when the time scale characterizing dynamics of scalar blobs is much shorter than the time scale of variation of the mixing rate, which is hardly the case in chemical engineering practice. For example in the case of mixing in the viscous-convective subrange the criteria for equilibrium were given by Chasnov25 including high Reynolds number and
y)
( ) ν 〈〉
1 dRD ,1 RD dt
1/2
(14)
In most of the mixers applied in chemical engineering applications the Reynolds number is not very high (in the stirred tank there are usually regions where Re is not high enough), which calls for using the nonequilibrium modeling. In this work the multiple-time-scale model proposed by Bałdyga24 for mixing of liquids of Sc . 1 is used. In the model the local value of the concentration variance σS2 has been divided into three parts according to the scale of segregation and related mechanism of mixing, namely, inertial-convective (σ12), viscous-convective (σ22), and viscous-diffusive (σ32):
σS2 ) σ12 + σ22 + σ32
(15)
Distribution of variances (σi2) is calculated from equations equivalent to eq 8 but formulated separately for σ12, σ22, and σ32.
[
]
∂σi2 ∂σi2 ∂σi2 ∂ ) (D + DT) + RPi - RDi (16) + 〈uj〉 ∂t ∂xj ∂xj ∂xj where RPi and RDi stand for production and dissipation terms. The inertial-convective variance σ12 is produced in the model from the macroscopic inhomogeneity of 〈f〉 as a result of velocity fluctuations
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5345
RP1 ) -2〈u′j f ′〉
( )
∂〈f〉 ∂〈f〉 ) 2DT ∂xj ∂xj
2
2
(17)
with the 〈f〉 distribution computed using the model equation for the average concentration of the passive scalar:
[
]
∂〈f〉 ∂〈f〉 ∂〈f〉 ∂ ) (D + DT) + 〈uj〉 ∂t ∂xj ∂xj ∂xj
(18)
The rate of decay of σ12 (and at the same time the production of σ22 because the variance is conserved) results from
σ12 〈〉 RD1 ) RP2 ) ) σ12R τS k
(19)
3 2/3 L 2 k τS = ) R‚〈〉 R‚〈〉1/3
(20)
where
represents the time constant for the inertial-convective mixing. L is the scale of large energy containing eddies, L ) u3/〈〉, that depends on the scale of the system. Accepting an estimate R ) 2 for the time-scale ratio for decay of velocity and concentration fluctuations which is commonly used in the commercial CFD codes one gets from eq 20
τS ≈
3L2/3 k ) 2〈〉 4〈〉1/3
(21)
The viscous-convective mixing can be described by
RD2 ) RP3 ) Eσ22
τG =
(0.03η) Sc = 17000E D
This results in a closed set of the balance equations for the concentration variance. After eqs 16 for i ) 1, 2, and 3 are added up, the resulting equation describes evolution of the complete variance σS2, σS2 ) σ12 + σ22 + σ32
[
( )
(23)
where E is an engulfment parameter. This yields the time constant for the viscous-convective mixing equal to
( )
1 ν τE ) ) 17.2 E 〈〉
1/2
(24)
in agreement with considerations following eq 13. Decreasing the scale of segregation by viscous deformation increases the wavenumbers and conveys the variance into the viscous-diffusive subrange of the concentration spectrum where mixing on the molecular scale occurs by molecular diffusion in deforming slabs
RD3 ) Gσ32
(25)
( )
∂〈f〉 2DT ∂xj
τC )
C1VR Ndimp3
τD )
Qfeed 〈u〉DT
(
)
and the related time constant reads
(29)
(30)
Before starting modeling, it is very useful to check which of the mixing subprocesses can affect the course of chemical reaction and precipitation. To characterize chemical reaction, one needs to use the adequate characteristic time constant; for example for a single second-order reaction
A + νBB f products
(31)
the time constant should characterize decay of the 〈cA〉‚〈cB〉 term resulting from the effect of chemical reaction, to be compared with the time constants for mixing characterizing decay of 〈c′Ac′B〉. Bałdyga and Bourne1 proposed for such a case
τR )
1 νBk〈cA〉 + k〈cB〉
(32)
In the case of precipitation one considers1,15 the characteristic time for nucleation
τN ) (26)
- Gσ32 (28)
Another time scale to be considered characterizes mesomixing by turbulent diffusion, τD, that characterizes turbulent dispersion of fresh feed stream close to the feed point. Assuming that the inlet from the feed pipe can be represented by a point source, one gets26
where
17050 17000E E= G ) 0.303 + Sc Sc
2
which shows that the only mechanisms of creation of the concentration variance is turbulent diffusion and that the only mechanism of mixing on the molecular scale is molecular diffusion in deforming structures. Equations 20, 24, and 27 represent the time constants for subgrid mixing; to have a complete set of time constants for mixing, we should mention here the time constant for macromixing that characterizes blending of fluids resulting from the macroscopic flow in the system. A simple characteristic of macromixing in a stirred tank is the circulation time τC, which equals the tank volume VR divided by the circulation capacity
(22)
1/2
]
∂σS2 ∂σS2 ∂σS2 ∂ + 〈uj〉 ) (D + DT) + ∂t ∂xj ∂xj ∂xj
and
〈〉 E ) 0.058 ν
(27)
NC RN
(33)
being the time necessary to form a certain characteristic number of nuclei per unit volume. In eq 33 RN is the
5346
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
nucleation rate for fresh premixed feed and NC (average number concentration of crystals in the bulk) can be estimated as
NC )
〈cA〉0M FL h3
(34)
and the estimate of the average crystal size, L h , should be known from experiments. For the crystal growth the characteristic time is expressed in relation to the rate of product concentration decrease resulting from the crystal growth
FG -1 τGcr ) A 〈cC〉 M g
(
)
(35)
where 〈cC〉 is the concentration of precipitating substance. Using the above presented time constants for mixing, one can classify the process (i.e., chemical reaction, precipitation etc.) from the point of view of competition between the process kinetics and mixing for various mechanisms of mixing considered. Mixing subprocesses with time constants comparable or larger than the characteristic time constants for the process can affect the course of considered process and thus should be considered in modeling. Comparison of time constants is also useful when choosing the adequate set of test reactions. 3. Mixing Effects on the Course of Parallel Test Reactions Carried Out in the Semibatch Stirred-Tank Reactor In this section the authors are interested in the influence of mixing on the course of parallel chemical reactions forming a test system based on neutralization of sodium hydroxide (A) and hydrolysis of ethyl chloroacetate (C) with hydrochloric acid (B) as given by eqs 1 and 2 and carried out in the semibatch stirred-tank rector. Semibatch reactors are often applied for production of fine chemicals, and the problems of selectivity (this section) and the product particle size are then of importance, as the product purity and size distribution determine its quality and price. 3.1. Experimental Procedure. The experimental procedure applied in this work was described in detail by Bałdyga and Makowski7 and Bourne and Yu.5 Experiments were carried out at T ) 293 K. At this temperature k1 ) 1.3 × 1011 dm3/(mol s) and from eq 32 one gets the characteristic time for first reaction equal roughly to 10-11 s, which is very short compared to all considered mixing times, and the first reaction can be classified as instantaneous or mixing controlled. The second reaction with k2 ) 23 dm3/(mol s) and characteristic time about 0.05 s is fast compared to inertial-convective and viscous-convective mixing but slow compared to the viscous-diffusive mixing. The stirred tank containing initially VBC,0 ) 0.02079 m3 of the premixture of HCL (B) and CH2ClCOOC2H5 (C), each concentration either 20 or 40 mol m-3, was operated semibatch, while the relative volume of NaOH solution added (cA0 equal either to 1000 mol m-3 or 2000 mol m-3) was 2% of the tank capacity (VA0 ) 4.16 × 10-4 m3). The feed point was located in the impeller discharge stream at the impeller disk level, at a distance of 0.09 m from the axis. The feed tube was of i.d. equal to 0.001 m and the feed time tf equal to 7, 8, 10, 15, and 20 min,
respectively. Measurements were carried out for the rotational stirrer speed in the range 106 rpm e N e 214 rpm, and the concentration of ethyl chloroacetate was measured before and after experiments chromatographically (HPLC). The final selectivity representing a relative amount of the substrate A converted into the byproduct S was expressed by
XS )
cC0 - cC
(36)
cA0
where cC represents the final ester concentration in the tank after the process and cA0 ) VA0cA0/(VA0 + VBC,0) and cC0 ) VBC,0cC0/(VA0 + VBC,0). 3.2. Modeling of Mixing Effects on the Course of Parallel Chemical Reactions. The problem of long feeding time in the semibatch process is solved following the work of Bałdyga and Makowski7 by discretizing the feed into σ portions of equal volume and solving the differential balance equations for mass, momentum, and species until the quasi-steady-state solution is obtained for each portion. This means that once the steady state is reached for the considered portion, then the next one is added and the procedure of reaching steady state is repeated. For each step the process is treated as the steady-state continuous-flow system; however, the flowing out mixture is collected during the time period tfeed/6 and added back to the system before the next portion is introduced. For each addition the k- model is first applied to calculate distributions of the average velocity, u j j, the average kinetic energy, k, and the average rate of the energy dissipation, . Afterward, the mixture fraction distribution, 〈f ′〉, is calculated. In the case of semibatch reactor of long feed time the mixture fraction is interpreted as a mass fraction of the fresh feed that is diluted by the older bulk fluid. It is so because the process is affected by mixing that takes place between the fresh feed and the bulk fluid rather than between fresh feed and fresh fluids present in the vessel initially. Both definitions of the mixture fraction are related for any portion of fresh feed by
f ′ ) (f - fb)/(1 - fb);
fb ) QAt/(VB0 + QAt) (37)
and different methods of calculation of 〈f ′〉are discussed by Bałdyga and Makowski.7 Figures 1a-c illustrate relation between f and f ′. The standard definition yields various distributions of 〈f〉 for different process time values, whereas for the high volume ratio of reagent solutions distribution of 〈f ′〉 is in practice independent of the process time. Using distribution of 〈f ′〉, it is possible now to compute the variances σ12, σ22, σ32, and resulting σS2 from the set of equations given by eq 8. Distribution of σS2 for mixing fresh feed with the bulk is given in Figure 1d. One can see that micromixing effects are observed in the plume of fresh feed close to the feed point. From the spatial distributions of 〈f ′〉 and σS2, one can express the local values of the probability density of f ′ using the Beta probability distribution
Φ(f ′,x b,t) )
(f ′)v-1(1 - f ′)w-1
∫01 yv-1(1 - y)w-1 dy
(38)
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5347
Figure 1. Visualization of the mixing zone in the semibatch stirred-tank reactor; N ) 214 rpm, tf ) 8 min. (a) Distribution of the mixture fractions, f, after 48 s; (b) distribution of the mixture fractions, f, after 432 s; (c) distribution of f ′ both after 48 and 432 s; (d) distribution of σS2 based on f ′ in the mixing zone.
with
[
v ) 〈f ′〉
〈f ′〉(1 - 〈f ′〉) 2
σS
] [
-1 ,w)
]
(1 - 〈f ′〉) 〈f ′〉
(39)
and calculate the local values of the average reaction rate between the ester, C, and the base, A.
∫01 cA(f ′)cC(f ′)Φ(f ′) df ′
〈r2〉 ) -k2〈cAcC〉 ) -k2
(40)
The local concentrations of all reactants of parallel chemical reactions are related to each other and to the mixture fractions f and f ′ by the relations
f)
cA - (cB + cC) + cB0 + cC0 , cA0 + cB0 + cC0 cA - (cB + cC) + cB,bulk + cC,bulk (41) f′) cA0 + cB,bulk + cC,bulk
Application together with eq 41 of the method of interpolation of cC(f) ) c∞C(f) + (cC - c∞C)[c0C(f) - c∞C(f)]/ (cC0 - c∞C) between the concentrations of reagents of the instantaneous c∞C(f) and infinitely slow reaction c0C(f) and using the condition that A and B cannot coexist in the same fluid elements (neutralization reaction is treated as instantaneous) enables one to express all reactant concentrations in terms of mixture fraction and local value of conversion of reactant C. Hence, it is enough to consider the balance of C to obtain concentra-
Figure 2. Influence of the mean energy dissipation rate, 〈〉, on the selectivity, XS, observed in the semibatch stirred-tank reactor; cA0 ) 2 kmol m-3, cB0 ) cC0 ) 0.04 kmol m-3.
tion distributions of all reactants
[
]
∂〈cC〉 ∂〈cC〉 ∂〈cC〉 ∂ + 〈uj〉 ) (D + DT) + 〈r2〉 (42) ∂t ∂xj ∂xj ∂xj For comparison we have performed also computations neglecting micromixing effects, i.e., replacing cA(f ′) and cC(f ′) by 〈cA〉 and 〈cC〉 respectively in eq 40. Comparison of model predictions with experimental data is presented in Figures 2 and 3. The model predicts well the
5348
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
based on the mixture fraction f, the averaged form of the moment transformation of the population balance takes the form
{
}
∂〈mj〉 ∂〈upi〉 ∂〈mj〉 ∂〈mj〉 ∂ + 〈mj〉 ) DpT + + 〈upi〉 ∂t ∂xi ∂xi ∂xi ∂xi
∫01 G(f)mj-1Φ(f) df
0j〈RN〉 + j
Figure 3. Influence of the feed time tf on the selectivity, XS, observed in the semibatch stirred-tank reactor; effect of feed concentrations. N ) 214 rpm.
trends of influence of feeding time, power input, and reactant concentrations. Observed effects agree well with the time-scale analysis: increase of energy input decreases the time scale of mixing decreasing as well the selectivity, increase of the feed time decreases the time constant for turbulent diffusion decreasing again selectivity, and finally decrease of reagent concentrations increases the time constant for chemical reaction and thus decreases selectivity. It is also shown in Figure 3 that neglecting micromixing effects leads to much worse simulation results than the ones obtained including micromixing effects (compare two broken lines in Figure 3) in agreement with earlier observations of Bałdyga and Makowski.7 4. Mixing-Precipitation Modeling of Barium Sulfate Precipitation In this section results of application of the model described in sections 2 and 3 and adapted to the problem of precipitation will be compared with experimental data of Podgo´rska.27 Podgo´rska carried out experiments in a Rushton type single-feed semibatch stirred-tank reactor of diameter T ) 0.242 m and impeller diameter dimp ) 0.075 m, which was operated in a semibatch manner. The feed position considered in what follows was placed near the impeller (r ) 0.060 m, z ) 0.080 m) and the feed pipe i.d. was equal to 0.0015 m. The final volume, being the sum of volumes of barium chloride and sodium sulfate solutions, was kept the same in all experiments while the volume ratio of the solution of barium chloride present in the vessel initially to the volume of sodium sulfate solution fed to the system, R, varied from 20 to 200. Stoichiometrically equivalent amounts of reactants were applied in experiments. The crystal size distributions were measured after the experiment using a Coulter-Counter and particles where photographed with the use of the scanning microscope. No agglomeration of crystals was observed. The closure scheme presented above for fast reactions was adapted to the problem of precipitation by Bałdyga and Orciuch16 and applied to interpret the precipitation process in the tubular reactor. The closure was developed using the moment transformation of the population balance for one internal coordinate L; after Reynolds averaging and application of the conditional closure
(43)
where DpT is the local value of the particle turbulent diffusion coefficient. As shown by Bałdyga and Orciuch28 in the case of small particles of BaS04 the diffusion coefficient for particles DpT does not differ significantly from the diffusion coefficients for the fluids DT. In eq 43 the averaged values of moments are defined by averaging over all fluid elements observed at a given position in the system at steady state and identified using the mixture fraction f:
b,t)〉 ) 〈mj(x
∫0∞ mj(xb,t,f)Φ(xb,t,f) df ) ∫0∞ 〈Ψ(xb,t,L)〉LjdL for j ) 0, 1, 2, ...
(44)
Notice that the regular definition of moments is based on averaging of population of crystals present in fluid elements identified by f
b,t,f) ) mj(x
∫0∞ Ψ(xb,t,L,f)Lj dL
(45)
and the local average density of the particle size distribution is defined by
〈Ψ(x b,t,L)〉 )
∫01 Ψ(xb,t,L,f)Φ(xb,t,f) df
(46)
Equation 43 should be solved together with the species balances
[
]
∂〈cR〉 ∂〈cR〉 ∂〈cR〉 ∂ + 〈ui〉 ) (DR + DT) ∂t ∂xi ∂xi ∂xi ka Fp 1 G(f)m2(f)Φ(f) df (47) 2 MC 0
∫
where R ) Ba2+, SO42-, Na+, Cl-. Of course the rate of precipitation for the sodium and chloride ions is equal to zero (hence G ) 0 for Na+ and Cl-), but concentrations of these species should be calculated to determine the local values of ionic strength and resulting activity coefficients. To express the local values of moments in fluid elements identified by f, the local, instantaneous values of a solid-phase concentration in the suspension, cC(f), are calculated from the local balance, cC(f) ) c0A(f) - cA(f). Assuming that the solid-phase concentration is proportional to such quantities as a number of crystals per suspension unit volume, a cumulative length per unit volume, a crystal surface area per unit volume, and a mass of crystals per unit volume, one gets the final closure relation for moments:
mj/〈mj〉 ) cC(f)/〈cC〉 The problem of long feeding time in the semibatch process is solved in the case of the precipitation process similarly to the case of parallel chemical reactions, i.e., by discretizing the feed into σ portions of equal volume replacing mixture fraction f by f ′ from eq 37 and solving
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5349
the differential balance equations for mass, momentum, species, and population balances. However, in this case the quasi-steady-state solution is obtained for each portion for the meso- and micro-mixing zone, whereas the bulk is treated as the semibatch well-mixed crystallizer where the species concentrations and the population of crystals vary continuously in time. 〈RN(f)〉 appearing in eq 43 represents the averaged sum of homogeneous and heterogeneous nucleation rates expressed by equations of the following form:
RN ) Rmax exp
( ) -A ln2 SA
(48)
where A ) 43.13, Rmax ) 5.36 × 1011 m-3 s-1 for heterogeneous nucleation and A ) 3137, Rmax ) 1.24 × 1049 m-3 s-1 for homogeneous nucleation as given by Vicum et al.29 The supersaturation ratio applied in eq 48 was calculated from eq 49
SA )
(
)
cBa2+cSO42KSP
0.5
γ(
(49)
with log(KSP) ) -10.05 and the activity coefficient calculated from the formulas derived by Bromley30
1 1 + (FBa2+ + FSO42-) log γ( ) -2.044xI 2 1 + xI
(50)
where the Fi terms account for interactions between unlike charged ions in the solution.30 In calculations it has been taken into account that the barium sulfate does not dissociate completely in aqueous solutions forming ion pairs, which decreases the supersaturation ratio compared to the case of complete dissociation observed in the case of very strong electrolytes31 2+ 2Ba(aq) + SO4(aq) a BaSO4(aq)
(51)
In the considered case, dissociation is not complete and thus modeled using the following equilibrium equation
KIP )
cBa2+cSO42-γ(2 cBaSO4(aq)γBaSO4(aq)
(52)
where log(KIP) ) 2.72. The local values of rate of growth were calculated from the two-step model that includes the surface integration step
G ) kr
[(
)
cBa2+,intcSO42-,int KSP
0.5
]
2
γ(,int - 1
) kr(SA,int - 1)2 (53)
and the molecular diffusion step
G ) kD(ci,bulk - ci,int) 2+ 20 , SO4(aq) , BaSO4(aq) (54) for i ) Ba(aq)
with kr ) 9.1 × 10-12 m s-1 (Vicum et al.29) and kD ) 4.0 × 10-5 (m s-1)(m3 kmol-1) (Bałdyga et al.15).
Figure 4. Effects of initial concentration on the particle size: N ) 3 s-1, R ) 50, tfeed ) 1650 s.
The average particle size was recalculated from the moments of the particle size distribution using the following equation
d43 )
m4 1 m φv
(55)
3
where φv represents the particle sphericity; this value was measured by Podgo´rska,27 and for the single-feed semibatch process the value φv ) 0.727 was obtained. Figure 4 compares measured end-predicted effects of initial concentration on the mean particle size, with cA0 ) VA0cA0/(VA0 + VB0) and cB0 ) VB0cB0/(VA0 + VB0). The trend of decreasing the particle size with increasing concentration is quite well-predicted in the region of high concentrations. In the region of small concentrations a clear trend of increase of particle size with increasing concentration is predicted while in experiments this effect is much smaller. This result is very interesting because micromixing effects are expected only at high concentrations (the time constant for nucleation is then much smaller than the time scale for micromixing by engulfment), which supports accuracy of predictions of the micromixing model. In the region of small concentrations the time constants for turbulent diffusion, nucleation, and growth take similar values, and because the rate of heterogeneous nucleation that dominates nucleation kinetics at low supersaturation is less sensitive to the supersaturation than the rate of crystal growth, increasing of concentration increases the size of crystals. It seems thus that the lack of predictability at low concentration results rather from inaccuracy of nucleation kinetics than inaccuracy of the mixing model. Moreover, as stated by Vicum et al.32 application of the thermodynamic Pitzer model33 to predict the activity coefficients can possibly improve model predictions; application of this model will, however, increase drastically computation time. Figure 5a-c explains some of the effects observed in Figure 4. Comparison of Figure 5a and 5b shows that the supersaturation decreases in the system during the process, and as a result most of the particles are nucleated in the initial stages of feeding. Figure 5c shows that in the case of high concentrations the supersaturation ratio is high even at the end of the process (at the beginning of the process it was as high as 107) and as a result one
5350
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
Figure 5. Distribution of supersaturation ratio in the tank: (a) at t ) 500 s for cA0 ) cB0 ) 0.0045 kmol m-3, tfeed ) 1650 s, N ) 3 s-1 and R ) 50; (b) at t ) 1570 s for cA0 ) cB0 ) 0.0045 kmol m-3, tfeed ) 1650 s, N ) 3 s-1 and R ) 50; (c) at t ) 1570 s for cA0 ) cB0 ) 0.0200 kmol m-3, tfeed ) 1650 s, N ) 3 s-1 and R ) 50.
Figure 6. Effect of stirrer speed on the mean particle size; cA0 ) cB0 ) 0.0045 kmol m-3, R ) 50, tfeed ) 1800 s.
observes a larger number of smaller particles. Effects of stirrer speed (Figure 6) and solution volume ratio (Figure 7) were measured only at small concentrations; both figures show rather good agreement between experimental data and model predictions, although a decrease of particle size with increase of the local value of the supersaturation ratio caused by increasing the volume ratio (SA increases by factor 1.5 for increase of R by factor 10 as observed in the distributions similar to that shown in Figure 5) is not predicted by the model. This again shows that the sensitivity of nucleation rate
Figure 7. Effect of volume ratio R on the mean particle size; cA0 ) cB0 ) 0.0045 kmol m-3, N ) 3 s-1, qfeed ) 0.137 × 10-6 m3 s-1.
to the supersaturation is too small at small concentration. Finally, Figure 8 shows comparison of the crystal size distributions measured and predicted by the model, whereas the predicted distribution has been recovered from the first 6 predicted moments using the statistically most likely probability distribution.21 λ
Ψnλ(L) ) B0 exp(
∑ AnLn) n)0
(56)
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5351
Figure 8. Predicted and measured crystal size distributions for cA0 ) cB0 ) 0.0045 kmol m-3, N ) 1.5 s-1, R ) 50, tfeed ) 1800 s.
5. Conclusions It has been shown in this work that the relatively simple CFD based multiple-time-scale mixing model can be used to predict most of the trends observed in experiments carried out for the fast parallel chemical reactions and precipitation of the solid products from the liquid ionic solutions. It has been further shown that in particular in the case of the precipitation process identification of the exact thermodynamics and kinetics of the crystallization is more important for improvement of model predictions than refinement of the mixing model. Interesting that a simple analysis of time constants for mixing and process kinetics enables prediction and understanding of trends observed in experiment and predicted by the model. Nomenclature ci, cR ) concentration of substance i or R, kmol m-3 ci,bulk ) concentration of i at the bulk, kmol m-3 ci,int ) concentration of i at the crystal surface, kmol m-3 c∞i ) concentration of reactant of instantaneous reaction, kmol m-3 0 ci ) concentration of nonreacting tracer, kmol m-3 D ) molecular diffusivity, m2 s-1 DT ) turbulent diffusivity, m2 s-1 DpT ) particle turbulent diffusivity, m2 s-1 d ) particle diameter, m d43 ) mass weighed particle diameter, m, µm E ) engulfment parameter, s-1 f ) dimensionless concentration of nonreacting tracer (mixture fraction) f ′ ) local renormalized mixture fraction G ) molecular diffusion parameter, s-1 G ) linear crystal growth rate, m s-1 I ) ionic strength, kmol m-3 KIP ) ion-pair equilibrium constant, kmol m-3 KSP ) thermodynamic solubility product, kmol2 m-6 k ) kinetic energy of turbulence, m2 s-2 k1, k2 ) second-order rate constant, m3 kmol-1 s-1 ka ) surface shape factor kB ) Batchelor wavenumber, m-1 L ) integral scale of turbulence, m L ) characteristic crystal size, m MC ) molar mass of product C, kg kmol-1 mj ) jth moment of particle size distribution, mj m-3 N ) stirrer speed, rpm, s-1
NC ) average number concentration of crystals in the bulk, m-3 RN ) nucleation rate, m-3 s-1 SA ) thermodynamic saturation ratio Sc ) Schmidt number t ) time, s tf, tfeed ) feeding time, s 〈u〉 ) local average velocity, m s-1 ui ) Eulerian fluid velocity component, m s-1 upi ) Eulerian particle velocity component, m s-1 XS ) product distribution V ) volume, m3 〈〉 ) local average value of the rate of energy dissipation, m2 s-3 〈〉 ) average rate of energy dissipation in the tank, m2 s-3 γ( ) mean activity coefficient for salt ν ) kinematic viscosity, m s-2 νT ) turbulent viscosity, m s-2 Fp ) crystal density, kg m-3 σ ) feed discretization σS2, σi2 ) dimensionless concentration variances τG ) crystal growth characteristic time, s τN ) nucleation time, s τR ) reaction time, s τC ) circulation time, s τK ) Kolmogorov time microscale, s τS ) time scale for the inertial-convective mixing, s Φ(f) ) Beta probability distribution φv ) crystal sphericity Ψ ) population density, m-4
Literature Cited (1) Bałdyga, J.; Bourne, J. R. Turbulent mixing and chemical reactions; Wiley: Chichester, 1999. (2) Bourne, J. R.; Moergeli, U.; Rys, P. Mixing and fast chemical reaction: Influence of viscosity on product distribution. 2 Europ. Conf. Mixing, BHRA, Cranfield, 1977; p 41. (3) Bourne, J. R.; Kut, O. M.; Lenzner, J. An improved system to investigate micromixing in high-intensity mixers. Ind. Eng. Chem. Res. 1992, 31, 949. (4) Bałdyga, J.; Bourne, J. R. The effect of micromixing on parallel reactions. Chem. Eng. Sci. 1990, 45, 907. (5) Bourne, J. R.; Yu, S. Investigation of micromixing in stirred tank reactors using parallel reactions. Ind. Eng. Chem. Res. 1994, 33, 41. (6) Bałdyga J., Roz˘ en´, A.; Mostert, F. A model of laminar micromixing with application to parallel chemical reactions. Chem. Eng. J. 1998, 69, 7. (7) Bałdyga, J.; Makowski, Ł. CFD modelling of mixing effects on the course of parallel chemical reactions carried out in a stirred tank. Chem. Eng. Technol. 2004, 27, 225. (8) Villermaux, J.; Falk, L.; Fournier, M. C.; Detrez, C. Use of parallel competing reactions to characterize micromixing efficiency. AIChE Symp. Ser. 88 1992, No. 286, 6. (9) Bałdyga, J.; Bourne, J. R.; Walker, B. Nonisothermal micromixing in turbulent liquid-theory and experiment. Can. J. Chem. Eng. 1998, 76, 641. (10) Bałdyga, J.; Henczka, M. Turbulent mixing and parallel chemical reactions in a pipe - application of a closure model. Proceedings of The 9th European Conference on Mixing, Paris, France, 1997; p 341. (11) Vicum, L.; Ottiger, S.; Mazzotti, M.; Makowski, Ł.; Bałdyga, J. Multi-scale modelling of a reactive mixing process in a semibatch stirred tank. Chem. Eng. Sci. 2004, 59, 1767. (12) Bałdyga, J.; Henczka, M.; Makowski, Ł. Effects of mixing on parallel chemical reactions in a continuous-flow stirred-tank reactor. Trans. Inst. Chem. Eng., Part A 2001, 79, 895. (13) Garside, J.; Tavare, N. S. Mixing, reaction and precipitation limits of micromixing in an MSMPR crystallizers. Chem. Eng. Sci. 1985, 40, 1485. (14) Pohorecki, R.; Bałdyga, J. The influence of intensity of mixing on the rate of precipitation. In Proceedings of Industrial Crystallization ‘78; de Jong, J., Jancˇic´, S. J., Eds.; North-Holland: Amsterdam, 1979; p 249.
5352
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
(15) Bałdyga, J.; Podgo´rska, W.; Pohorecki, R. Mixing-precipitation model with application to double feed semibatch precipitation. Chem. Eng. Sci. 1995, 50, 1281. (16) Bałdyga, J.; Orciuch, W. Closure problem for precipitation. Trans. Inst. Chem. Eng., Part A 1997, 75, 160. (17) Villermaux, J.; David, R. Effect du micromelange sur la precipitation. J. Chim. Phys. 1988, 85, 273. (18) Barresi, A. A.; Marchisio, D.; Baldi, G. On the role of microand mesomixing in a continuous Couette-type precipitator. Chem. Eng. Sci. 1999, 54, 2339. (19) Pipino, M.; Barresi, A. A.; Fox, R. O. A PDF approach to the description of homogeneous nucleation. In Proceedings of Quarto Convegno Internazionale Fluidodinamica Multifase, Nell’impiantistica Industrialles, Ancona, Italy, 1994. (20) Wei, H.; Garside, J. Application of CFD modelling to precipitation systems. Trans. Inst. Chem. Eng., Part A 1997, 75, 219. (21) Bałdyga, J.; Orciuch, W. Barium sulphate precipitation in a pipe - an experimental study and CFD modelling. Chem. Eng. Sci. 2001, 56, 2435. (22) Bałdyga, J.; Jasin´ska, M.; Orciuch, W. Barium sulphate agglomeration in a pipe -an experimental study and CFD modelling. Chem. Eng. Technol. 2003, 26, 334. (23) Marchisio, D. L.; Barresi, A. A.; Garbero, M. Nucleation growth and agglomeration in barium sulphate turbulent precipitations. AIChE J. 2002, 48, 2039. (24) Bałdyga, J. Turbulent mixer model with application to homogeneous, instantaneous chemical reactions. Chem. Eng. Sci. 1989, 44, 1175.
(25) Chasnov, J. R. The viscous-covective subranges in nonstationary turbulence. Phys. Fluids 1998, 10, 1191. (26) Bałdyga, J.; Bourne, J. R. Interaction between mixing on various scales in stirred tank reactors. Chem. Eng. Sci. 1992, 47, 1839. (27) Podgo´rska, W. Effects of mixing on precipitation. (in Polish) Ph.D. Thesis, Warsaw University of Technology, 1993. (28) Bałdyga, J.; Orciuch, W. Some hydrodynamics aspects of precipitation. Powder Technol. 2001, 121, 9. (29) Vicum, L.; Mazzotti, M.; Bałdyga, J. Applying thermodynamic model to the non-stoichiometric precipitation of barium sulfate. Proc. 15th Int. Symp. Ind. Cryst. Sorento, 2002. (30) Bromley, L. A. Thermodynamic properties of strong electrolytes in aqueous solutions. AIChE J. 1973, 19, 313. (31) Felmy, A. R.; Rai, D.; Amonette, J. E. The solubility of Barite and celesite in sodium sulfate: evolution of thermodynamics data. J. Solution Chem. 1990, 19, 175. (32) Vicum, L.; Mazzotti, M.; Bałdyga, J. Applying a thermodynamics model to the nonstoichiometric precipitation of barium sulfate. Chem. Eng. Technol. 2003, 26, 325. (33) Pitzer, K. Activity coefficients in electrolyte solutions; CRC Press: Boca Raton, FL, 1991.
Received for review September 2, 2004 Revised manuscript received November 8, 2004 Accepted November 12, 2004 IE049165X
