Comparison of Various Micromixing Approaches for Computational

Ind. Eng. Chem. Res. 2009, 48, 999–1007

999

Comparison of Various Micromixing Approaches for Computational Fluid Dynamics Simulation of Barium Sulfate Precipitation in Tubular Reactors ¨ ncu¨l,* Ga´bor Janiga, and Dominique The´venin Alper A. O Institut fu¨r Stro¨mungstechnik und Thermodynamik, Lehrstuhl fu¨r Stro¨mungsmechanik und Stro¨mungstechnik (ISUT/LSS), Otto-Von-Guericke-UniVersita¨t-Magdeburg, UniVersita¨tsplatz 2, 39106 Magdeburg, Germany

Precipitation (also called reactive crystallization) processes are widely employed in the chemical industry in a variety of reactors and for different purposes. The properties and quality of the resulting particulate product are highly dependent on the mixing conditions in the reactor, since precipitation reactions are mostly fast: the reaction time is usually of the order of or shorter than the mixing time. Therefore, experimental and numerical investigations of precipitation processes under different hydrodynamic conditions constitute an interesting topic. In this work, precipitation of barium sulfate (BaSO4) is numerically investigated in two different tubular reactors, taking micromixing effects into account. Different models are employed with an increasing level of complexity. These models and the reaction kinetics are first explained. RANS-based axisymmetric computations are performed using the industrial CFD code FLUENT 6.3 at steady-state. The results of the simulations are then presented and discussed in detail. Moreover, quantitative comparisons with experimental data obtained from literature are also employed for validation. Finally, concluding remarks and possible model improvements are listed. 1. Introduction The term “micromixing” refers to the molecular mixing process which takes place at a length scale below Kolmogorov and Batchelor scales.1-3 Micromixing is an important mechanism that directly influences fast chemical reactions since the reactants locally come in contact with each other through micromixing. In this respect, micromixing must be distinguished from “macromixing”, which refers to the flow mixing processes on the scale of the whole reactor, controlling the mean concentration and the global residence time distribution.3 In that case, the processes between these two scales, which refer to the inertial-convective part of the concentration spectrum and to the extent of turbulent diffusion, should be called “mesomixing”. Notice that especially in the combustion literature mixing in the inertial-convective part of the spectrum is called micromixing as well, although the size of the corresponding eddies can be as large as the integral scale of turbulence. Turbulent mixing processes are nowadays simulated by computational fluid dynamics (CFD) for many engineering applications (e.g., precipitation, combustion, etc.). Although the macromixing process is easily solved using appropriate transport equations, micromixing must be modeled additionally in order to predict its effect on the reactive mixing processes when applying the standard Reynolds-averaged Navier-Stokes (RANS) equations, relying on a relatively coarse grid. For this purpose, a number of micromixing models have been formulated and presented in the literature.2,4-7 The most common models and methods for implementing them can be designated as ME (multienvironment), Eng (engulfment), EDD (engulfment-deformation-diffusion), IEM (interaction by exchange with the mean), and DQMOM-IEM (direct quadrature method of momentssinteraction by exchange with the mean). They can be formulated either in a Lagrangian or an Eulerian framework. The aim of this work is to implement and test the Eng, ME, and DQMOM-IEM micromixing models in order to describe the precipitation process of BaSO4 in two coaxial pipe mixers. * To whom correspondence should be addressed. Tel.: +49-3916718195. Fax: +49-391-6712840. E-mail address: [email protected].

Although the Eng and ME models have been applied to the prediction of the precipitation processes in tubular reactors for many years, the DQMOM-IEM micromixing model has still not been applied widely, since it is quite recent. Therefore, the predictions of this recent model in comparison to those of other approaches represent an interesting aspect. The employed numerical model is based on complementary elements. The population balance equation (PBE) is solved using the standard method of moments (SMOM), in which the crystal size distribution (CSD) is expressed in terms of its lower-order moments. Moreover, two presumed probability density function (PDF) formulations have been applied. One of these is the mixture fraction approach, in which the reactant concentrations are defined as a function of local mixture fraction values.8 The second one is the finite mode (FM)-PDF method, which is capable of representing the composition PDF by a finite set of delta functions corresponding to different environments within each volume element.5 Both methods make it possible to perform subgrid level computations with environment-based micromixing approaches. All these methods are explained in details in the next section. Similar numerical studies in tubular reactors have already been made in our research group for different configurations and hydrodynamic conditions, but mostly without considering ¨ ncu¨l et micromixing effects (see, for example, the work of O al.9). More recent computations include micromixing effects, employing various models. For validation, these numerical results are compared with experimental data obtained from literature. In the present work, corresponding comparisons are presented and discussed. 2. Numerical Models 2.1. Micromixing Models. In order to understand the idea of defining the micromixing model in terms of “fluid environments”, it is useful to consider the ME model in more details. The ME model dates back to the 1970s10 and is still being applied for many studies in the form of two-,2,11 three-,12,13 and four-environment systems.14,15 The basic principle of this

10.1021/ie800364k CCC: $40.75  2009 American Chemical Society Published on Web 11/26/2008

1000 Ind. Eng. Chem. Res., Vol. 48, No. 2, 2009 Table 1. Physical Properties of the Environments and Representation of Their Interaction (Arrows)a

a

Abbreviations: iuf, initially unmixed fluid; mf, mixed fluid; pf, pure fluid; r, reaction.

approach is to divide the composition space into a number of environments which interact through micromixing.6 These environments are characterized by their volume fractions and species concentrations. In this way, subgrid level fluctuations in each computational cell can be considered. Although the ME model might be quite practical for many numerical investigations it may also show some limitations due to its simplified structure. For instance, it is known that this idealized model does not always describe properly the physics of mixing2 and may be inadequate for complex chemistry.6 In the present work, micromixing approaches based on two (denoted 2E) and three environments (3E) are employed. The properties of the environments in 2E and 3E models are explained next in order to clarify the procedure and the notations. The first environment represents the unmixed inlet feed whereas the second one represents the pure reactor medium in the 2E model6 (see Table 1). Note that an initially unmixed fluid (iuf) is similar to a pure fluid (pf) in principle but can only exist as a boundary condition. The whole reactor consists initially of Environment 2. Due to injection Environment 1 increases with time and mixes with the second environment. Consequently, the volume fraction of Environment 2 decreases. Chemical reactions are only allowed in Environment 1. In fact, the Eng model is also formally a 2E model from this point of view, but with another formulation of the micromixing rate, which is based on an engulfment parameter.2 Note that in principle this standard Eng model can only be applied for slow feeding, when the inertial-convective mixing does not affect too much the process. Otherwise, one should apply instead the extended version of the Eng model,16 including two time constants, where the controlling mechanism shifts from the viscous-convective (small Re) to the inertial-convective (large Re) range. In the 3E model, Environments 1 and 2 contain pure (unmixed) fluids. For an easier understanding, these two environments can be considered as two separate feed streams entering the reactor, in a similar manner to the real configuration in this work. The concentrations in these environments remain constant and equal to the inlet concentration. On the other hand, their volume fractions diminish due to mixing leading to Environment 3, where the reaction occurs (Table 1). At this point, it is useful to give further information about the formulation and solution methods of the ME model in the flow computation. Each environment can be considered as existing at a particular spatial location (i.e., in physical space) with a certain probability.6 In order to define the probability of

each environment, a PDF is employed. The corresponding socalled joint composition PDF, denoted by fφ, is expressed as fφ(ψ;x, t)dψ ) P[ψR < φR(x, t) < ψR + dψR, R ) A, B, ... ] (1) where ψ is the composition vector, x is the spatial position vector, t is time, and φR(x, t) is the composition of scalar R. The PDF methods can be distinguished as full and presumed PDF methods. In the full PDF methods, the PDF is described with a set of notional particles that obey stochastic differential equations and that move within the computational domain according to the PDF transport equation written in physical and composition space. On the other hand, the functional form of the PDF is assumed a priori in presumed PDF methods, employing well-known functions (β-function, γ-function, Dirac or δ-function, etc.).17,18 The focus is on the FM (finite mode) PDF method in this work. The presumed ME model can thus be defined by the FMPDF method, where the joint composition PDF is approximated by a finite set of delta functions:6 Ne

fφ(ψ;x, t) )

Ns

∑ p (x, t)∏ δ(ψ n

R)1

n)1

R - 〈φR 〉n(x, t))

(2)

where pn(x, t) is the probability of environment (or mode) n at position x and time t, 〈φR〉n(x, t) denotes the value (or mean composition) of scalar R within environment n, while Ne and Ns are the total number of environments and of scalars, respectively. As can be seen from this formula, every cell of the computational domain contains Ne different environments with a certain probability. In this manner, micromixing processes can be taken into account at subgrid level during the CFD simulations. Equation 2 describing the FM-PDF can be rewritten as below when the scalar R is chosen to be the mixture fraction, ξ: Ne

f(ξ;x, t) )

∑ p (x, t)δ(ξ - 〈ξ〉 (x, t)) n

n

(3)

n)1

Implemented in this manner, the mixture fraction approach is indeed also a presumed PDF method, in which ξ ) 0 for one feed stream whereas ξ ) 1 for the other one. Thus, the fluid is “segregated” (or unmixed) when ξ ) 0 or ξ ) 1. On the other hand, the fluid is “micromixed” when 0 < ξ < 1. The degree of this micromixing can be quantified by eq 3 for each volume element in the CFD calculations. It is therefore possible to obtain local PDFs, i.e. local micromixing properties, everywhere in

Ind. Eng. Chem. Res., Vol. 48, No. 2, 2009 1001

the domain. Moreover, it is possible to reconstruct the local PDF by using the lower-order moments of the mixture fraction PDF.18 The zeroth-moment of this PDF (which represents the sum of the probabilities) is automatically unity: Ne

∑p

n)1

(4)

n)1

The local mean mixture fraction (or the first-moment of the mixture fraction PDF), 〈ξ〉, and the further moments of this PDF, denoted as 〈ξj〉, can be calculated by Ne

〈ξ〉 )

Ne

∑ p 〈ξ〉 ) ∑ s n

n

(5)

jg2

(6)

n

n)1

n)1

Ne

〈ξ 〉 ) j

∑ p 〈ξ〉 n

j n,

n)1

where sn is the weighted mixture fraction (sn ) pn〈ξ〉n). This last relation is used during the computation to determine the mixture fraction knowing the weighted mixture fraction. The mixture fraction variance (i.e., the square of the standard deviation of the mixture fraction PDF) is estimated by (7) 〈ξ′ 2 〉 ) 〈ξ2 〉 - 〈ξ〉2 The mixture fraction moment values are expected to become more accurate as the number of environments increases since the corresponding predictions converge toward the results of computations based on full PDF methods.7 The following relation between the mixture fraction (which is a conserved scalar) and local reactant concentrations can be obtained for a simple, instantaneous and irreversible reaction (A + B f P) resulting from the mixing of two nonpremixed feeds:19 ξ)

CA - CB + CB0

(8)

CA0 + CB0

where the reactant concentrations with subscript 0 denote the inlet values. In the case of a finite-rate, single chemical reaction, individual reactant concentrations can be expressed in terms of mixture fraction and reaction progress variable, Y, as follows CB CA ) ξ - ξsY, ) (1 - ξ) - (1 - ξs)Y CA0 CB0

(9)

ξs )

(10)

CA0 + CB0

For a nonreacting system Y ) 0. For an instantaneous reaction, since the two reactants cannot coexist at the same location, Y can be evaluated as follows17

(

Y ) Y∞(ξ) ) min

ξ 1-ξ , ξs 1 - ξs

)

where the micromixing rate, γ, and the spurious dissipation rate (which is a correction term necessary to eliminate the spurious scalar dissipation resulting from the finite-mode representation of turbulent diffusion5), γs, are ε γ ) Cφ k (13) 2Γt | ∇ 〈ξ〉1|2 γs ) 2 〈ξ〉1 where k and ε denote turbulent kinetic energy and dissipation rate, respectively. Differently from the 2E model, γ is expressed as follows in the Eng model


Vε

γ ) 0.058

where the stoichiometric mixture fraction, ξs, is defined as CB0

More precisely, the DQMOM-IEM model relies on the DQMOM model,23 that solves the Eulerian composition PDF transport equation, in which the molecular mixing term is closed by the IEM model. This model is also defined in terms of environments like the ME model. However, contrary to the ME model, reaction may take place in all environments of the DQMOM-IEM model due to the fact that there is fundamentally exchange of fluid between all environments (Table 1). Theoretically, as the number of environments increases the predictions of this model would approach those of an IEM-based full PDF method since the micromixing mechanisms are exactly the same.7 In this respect, the DQMOM-IEM model represents a promising numerical tool which would yield reliable results at a considerably lower computational cost than the full PDF method. As such, DQMOM-IEM should solve some, if not all, of the known issues and deficiencies associated with the simple IEM approach.7,24 The two- and three-environment (also called two- and three-node) formulations of this approach, denoted DQMOM-IEM(2E) and DQMOM-IEM(3E), are employed in the present work (Table 1). Both are again coupled with presumed PDFs in a standard manner. The source terms of the RANS-based transport equations of the above-mentioned mixing parameters are given now for the 2E2,6,7 and Eng2,6 models (note that the transport equation for s2 is not needed here explicitly since 〈ξ〉2 is constant due to the nonreacting condition of Environment 2; hence, the relation sn ) pn〈ξ〉n is employed): Sp1 ) γp1p2 - γsp1 Ss1 ) γp1s2 - γsp1〈ξ〉2 (12) Ss2 ) -γp1s2 + γsp1〈ξ〉2

(11)

Consequently, Y is zero in the feed streams where no reaction occurs. The source term for Y can be reformulated for a finiterate chemistry as in the case of the present work and is given in the next section. The second type of micromixing model used in this work is the DQMOM-IEM model introduced recently by Fox.6 This model builds on top of the classical IEM model, which is one of the simplest micromixing models, as described by Harada et al.,20 Costa and Trevissoi,21,22 and Villermaux and Devillon.4

(14)

where V is the kinematic viscosity. Similarly, for the 3E model2,6,7,25 (note that the transport equation for s1 (respectively s2) is not needed here explicitly since 〈ξ〉1 (respectively 〈ξ〉2) is constant due to the nonreacting condition of Environment 1 (respectively, 2); hence, the relation sn ) pn〈ξ〉n is employed): Sp1 ) -γp1(1 - p1) + γsp3 Sp2 ) -γp2(1 - p2) + γsp3 Ss1 ) -γ(1 - p1)s1 + γsp3〈ξ〉1 (15) Ss2 ) -γ(1 - p2)s2 + γsp3〈ξ〉2 Ss3 ) γ[(1 - p1)s1 + (1 - p2)s2] - γsp3(〈ξ〉1 + 〈ξ〉2) where 〈ξ′ 2 〉 ε k [p (1 - p )(1 - 〈ξ〉 )2 + p (1 - p )〈ξ〉 2] 1 1 3 2 2 3 2Γt 2 γs ) | ∇ 〈ξ〉3| 1 - 2〈ξ〉3(1 - 〈ξ〉3)

γ ) Cφ

(16)

1002 Ind. Eng. Chem. Res., Vol. 48, No. 2, 2009

Finally, for the DQMOM-IEM(2E) model6,7,26

The nucleation rate, that has been estimated by Bałdyga and Orciuch,29 using the experimental data of Nielsen,30 is as follows

Sp1 ) 0 Γt (p | ∇ 〈ξ〉1|2 + p2| ∇ 〈ξ〉2|2) ξ1 - ξ2 1 Γt Ss2 ) -γ(p1s2 - p2s1) (p | ∇ 〈ξ〉1|2 + p2| ∇ 〈ξ〉2|2) ξ1 - ξ2 1 (17)

Ss1 ) γ(p1s2 - p2s1) +

and for the DQMOM-IEM(3E) model6,7 Sp1 ) 0 Sp2 ) 0 Ss1 ) γ(p1(s2 + s3) - (1 - p1)s1) Γt [p | ∇ 〈ξ〉1|2(2ξ1 - ξ2 - ξ3) + + (ξ1 - ξ2)(ξ1 - ξ3) 1 (p2| ∇ 〈ξ〉2|2 - p3| ∇ 〈ξ〉3|2)(ξ2 - ξ3)] Ss2 ) γ(p2(s1 + s3) - (1 - p2)s2) (18) Γt 2 [p | ∇ 〈ξ〉2| (2ξ2 - ξ1 - ξ3) + + (ξ2 - ξ1)(ξ2 - ξ3) 2 (p3| ∇ 〈ξ〉3|2 - p1| ∇ 〈ξ〉1|2)(ξ3 - ξ1)] Ss3 ) γ(p3(s1 + s2) - (1 - p3)s3) Γt [p | ∇ 〈ξ〉3|2(2ξ3 - ξ1 - ξ2) + + (ξ3 - ξ1)(ξ3 - ξ2) 3 (p1| ∇ 〈ξ〉1|2 - p2| ∇ 〈ξ〉2|2)(ξ1 - ξ2)] where γ ) Cφ

ε 2k

(19)

In all of these equations, the micromixing constant, Cφ, is assumed to be unity, following the work of Marchisio et al.,27 who investigated the same chemical process in their tubular reactor. Moreover, the turbulent diffusivity constant, Γt, is defined as Veff Γt ) Sct

(20)

where Veff is the effective kinematic viscosity and Sct is the turbulent Schmidt number which is set to 0.7, a typical value used in the computations of such turbulent processes.6 The mixing variables, for which a transport equation is necessarily defined in the present study, are listed in Table 1. Note that transport equations for p2 in the Eng, 2E, and DQMOMIEM(2E) models and for p3 in the 3E and DQMOM-IEM(3E) models are not defined additionally in order to reduce the CPU time and effort. These variables can be obtained by considering the summation of the probabilities to unity as given in eq 4. Finally, the reaction progress variable, Y, is solved by considering transport equations identical with that describing the mixture fraction, but adding on the right-hand side the source term defined later in eq 25, since Y is not a conserved scalar. 2.2. Reaction Kinetics and PBE. The precipitation reaction of BaCl2 and Na2SO4 producing BaSO4, as shown in eq 21, is chosen as the chemical process in this work since this reaction has been widely examined in literature and in the previous works of the authors.9,28 The main mechanisms of this precipitation reaction are nucleation and growth. Possible aggregation and breakage processes are assumed to be negligible considering the obtained results in previously published investigations.27,29 BaCl2 + Na2SO4 f 2NaCl + BaSO4 V

(21)

[

]

A (22) ln2 S where A ) 44.60 and B ) 1.06 × 1012 m-3 s-1 for heterogeneous nucleation whereas A ) 3020 and B ) 1.50 × 1045 m-3 s-1 for homogeneous nucleation. The local nucleation rate is defined by the sum of homogeneous and heterogeneous nucleation rates. Here, S indicates the supersaturation ratio: J ) B exp -

(

S ) γ(

CBaCSO4 Ksp

)

0.5

(23)

where γ( is the activity coefficient and Ksp is the solubility product of BaSO4, that is equal to 1.1 ×10-10 kmol2 m-6 at 25 °C. The local activity coefficient values are computed according to Bromley’s method.31 The practical computation of all terms involved in the precipitation kinetics, using concentrations of reacting and nonreacting ions, has been described in a previous publication.9 The crystal growth rate is described by a two-step diffusionadsorption model as given below29

[( ) ]

G ) kg γ(

/ C/BaCSO 4

Ksp

0.5

2

-1

) kd(CBa - C/Ba)

/ ) kd(CSO4 - CSO ) 4

(24)

where the growth rate constant, kg, is equal to32 4.0 × 10-11 m s-1 and the mass transfer coefficient, kd, is fixed at 1.0 × 10-4 (m s-1)(m3 kmol-1). Here, the concentrations with stars denote the values at the crystal surface. This equation is solved by using a simple Newton-Raphson iterative algorithm. The source term for the reaction progress variable used in eq 9 can then be expressed for example as SY )

3kvFBaSO4µ2G ξsCBa0WBaSO4

(25)

where the crystal volume shape factor, kv, is taken equal to 5.0 and 58.0 for Cases 127 and 2,29 respectively, FBaSO4 and WBaSO4 indicate the density and the molecular weight of BaSO4 crystals, respectively, and µ2 is the second moment of the CSD. As stated previously, the PBE is solved employing SMOM in which the source terms of the moments (here, all moments up to the fourth-moment are considered) are defined as Sµj ) 0jJ + jGµj-1, j ) 0 ... 4

(26)

Using the moments, the volume-averaged mean crystal size values, d43, are calculated by the following formula to validate the numerical model, since corresponding reference experimental data for this parameter can be found in the literature d43 )

µ4 1 µ3 Ψ

(27)

The sphericity of the particles, Ψ, has been estimated via experimental observations27,29 and is equal to 0.333 and 0.208 for Cases 1 and 2, respectively. At this point it is worth reminding that the source terms given in eqs 25 and 26 are just the general expressions of those terms. Since the related scalars are defined with their probabilities and local fractions (similarly to the mixing parameters), their source terms must additionally be coupled with the corresponding

Ind. Eng. Chem. Res., Vol. 48, No. 2, 2009 1003

Figure 1. Tubular reactor used in this work. Dimensions (in millimeters) are as follows: l ) 168 (240), L ) 2100 (2040), Ri ) 0.5 (0.9), Ro ) 5 (16) for Case 127 (respectively, for Case 229).

micromixing and correction terms in each micromixing model. This can simply be managed by substituting the related scalar existing in environment n instead of 〈ξ〉n in the source term equations listed in the previous section. More detailed information on this procedure can be obtained in the literature.6,7,25,26 3. Experimental Work The configuration and dimensions of the tubular reactor (see Figure 1) are based on the experimental studies of Marchisio et al.27 and Bałdyga and Orciuch.29 The study of Marchisio et al. is denoted as Case 1 and that of Bałdyga and Orciuch as Case 2 in the whole paper. Two independent experiments are chosen in the present work in order to refine the validation of the applied numerical models. In Case 1, two different inlet configurations have been analyzed, both yielding a mean pipe Reynolds number Re of 10 000 (keeping the mean velocities of the inner and outer feed streams equal to each other). BaCl2 is fed through the inner tube when Na2SO4 flows in the annular region (conf A) and vice versa (conf B). In both configurations, the concentration of the reactant in the inner pipe is kept constant at 0.034 kmol m-3 whereas that of the other reactant is equal to 0.034R kmol m-3. The concentration ratio, R, takes values between 0.1 and 3.0. On the other hand, BaCl2 and Na2SO4 solutions are always fed in Case 2 through the inner and outer pipes, respectively, with various velocity values corresponding to a mean pipe Reynolds number Re ranging from 20 000 to 60 000 (keeping the mean velocity of the inner pipe equal to the mean velocity in the reactor) and with concentrations equal to 1.50 and 0.015 kmol m-3, respectively. 4. CFD Computations The number of cells in the employed structured grid has been set around 5000 and 50 000 for Cases 1 and 2, respectively. These numbers have been obtained after checking in detail the grid-independency of the obtained results for one set of parameters. Note that Case 2 has been found to be much more sensitive to the grid size compared to Case 1, probably due to a higher supersaturation ratio. The resulting two-dimensional axisymmetric simulations are carried out using the industrial CFD package FLUENT 6.3 at steady state. FLUENT is complemented by appropriate userdefined functions (UDFs) and scalars (UDSs), which are externally implemented in order to solve the reaction kinetics equations and the source terms of micromixing as well as of the PBE. The standard k-ε approach is employed as the turbulence model, since it should deliver reasonable accuracy in short computational times for such a simple flow pattern. The simulations are performed in three steps for an easier and faster convergence. The flow field is solved first, then further computations with the micromixing model are carried out, and finally the precipitation process is calculated. Desired convergence has been obtained after around 2000-5000 iterations for

Figure 2. Local mean mixture fraction profiles as a function of radial position for 3 different axial locations obtained without micromixing (lines) and with the 3E micromixing model (diamonds) at Re ) 40 000 in Case 2. The smaller plot at the upper right corner corresponds to a magnification of the profiles near the r-axis.

the first step, after around 1000-16 000 supplementary iterations for the second step, and after around 1500-25 000 final iterations for the third step depending on the hydrodynamic and reaction parameters as well as on the employed micromixing models. The under-relaxation factors for the scalars are set between 0.75-0.90 for the micromixing parameters and between 0.10-0.90 for the reaction and PBE parameters, depending on the specific conditions. Due to the numerical stiffness induced by fast chemical reactions, under-relaxation parameters as low as 0.10 are indeed needed in some cases to maintain stability. 5. Results and Discussion 5.1. Mixing Process. Before analyzing the results of the precipitation process, it is useful to compare the local mean mixture fraction values, 〈ξ〉, as a verification of the predictions of the micromixing models, based on the fact that the micromixing model leaves the mean unchanged for such a nonreacting, conserved scalar. Therefore, the 〈ξ〉 profiles obtained, for instance, by the 3E model at Re ) 40 000 in Case 2 are shown for various axial locations in Figure 2 in comparison to those obtained without micromixing. As can be seen easily, the profiles are identical. This confirms the correct functionality of the implemented micromixing model. The profiles obtained by the other micromixing approaches and at various Re values have also been checked in a similar manner for both cases. The spatial distribution of all parameters controlling micromixing in the reactor are presented in Figure 3 (again for Case 2 employing the 3E model and Re ) 40 000) as an example. This plot might help the reader to better understand the underlying models. As the two pure-fluid-containing environments (Environments 1 and 2) diminish, the third environment (Environment 3), where the reaction occurs, extends as a result of the mixing process. The resulting mixture fraction and weighted mixture fraction plots of Environment 3 (i.e., ξ3 and s3) are presented as well. 5.2. Precipitation Process. The mean crystal size values (d43) averaged at the reactor outlet in Case 1 are shown in Figure 4 as a comparison between different micromixing models and with reference to the experiments, considering various R values. Additionally, the corresponding averaged relative errors between the results of each micromixing model and experimental data are listed in Table 2. It is worth reminding that the experimental

1004 Ind. Eng. Chem. Res., Vol. 48, No. 2, 2009

Figure 3. Spatial distribution of all micromixing parameters according to the 3E model (Re ) 40 000) in Case 2. (a) Probability of Env 1, p1. (b) Probability of Env 2, p2. (c) Probability of Env 3, p3. (d) Mixture fraction of Env 3, ξ3. (e) Weighted mixture fraction of Env 3, s3. Note that all three probability contours (a-c) are identified with one common colormap. The contour plots are magnified 75 times in the radial direction for a better visualization, due to the high aspect ratio of the reactor.

Figure 4. Mean crystal size profiles averaged at reactor outlet in Case 1 according to our numerical calculations in comparison with the experimental data for (a) conf A and (b) conf B.

Ind. Eng. Chem. Res., Vol. 48, No. 2, 2009 1005 Table 2. Averaged Relative Error Values between Numerical Results and Experimental Data for Each Micromixing Model in Case 1sConf A (ErrC1-A) and Conf B (ErrC1-B)sas Well as in Case 2 (ErrC2)a model\relative error (%) 1E model (no micromixing) Eng model 2E model 3E model DQMOM-IEM(2E) model DQMOM-IEM(3E) model

ErrC1-A ErrC1-B ErrC2 Errmean 19.6 20.9 20.5 22.9 44.4 22.0

52.9 38.1 38.8 63.8 82.7 37.0

15.9 17.9 15.8 11.5 18.8 14.0

29.4 25.6 25.0 32.7 48.7 24.3

CPU time in hours (for reaction) 0.5 1.0 1.0 1.2 28.0 35.0

a Mean relative errors (Errmean) based on these three analyses and the computing (CPU) times for the simulation of the mixing-precipitation reaction are also listed for each approach.

data including the error bars are those presented in the work of Marchisio et al.27 The numerical average at reactor outlet is based on the local mass flowrate and should thus be directly comparable to the experimental results. As can be seen both qualitatively and quantitatively, the global agreement in conf A is more satisfactory than that in conf B. This discrepancy probably comes from the employed precipitation kinetics rather than from the micromixing approach. The retained kinetic equations are probably unable to predict fully the effects of inverse feed conditions, even if they are based on the local ionic activity. Therefore, this problem might possibly be solved by applying alternative and more sophisticated kinetic equations.26,33 At the lowest supersaturation ratio (i.e., R ) 0.1), the influence of the micromixing model is negligible due to the lack of competition between the mixing and reaction mechanisms. This situation can be understood better when the relevant timescales are compared. The reaction time scale, τr, until obtaining 90% chemical conversion under homogeneous conditions (computed as in the previous work of the authors34) for R ) 0.1 and 1.0 is equal to 61.3 and 12.0 s, respectively, whereas the mean micromixing time scale, τm ) 1/γ, is around 79.1 s. Since the slowest mechanism controls the process, it is obvious that micromixing plays a more significant role at R ) 1.0. Contrarily, the reaction and micromixing phenomena affect the whole process almost equivalently at R ) 0.1. Various micromixing approaches may yield quite distinct results at higher supersaturation levels. For instance, the difference between the results of 2E and DQMOM-IEM(2E) models at R ) 1.0 in both configurations is more than 2-fold. However, the Eng, 2E, and DQMOM-IEM(3E) models predict almost the same results in this case. Additionally, these models do not yield a clear peak around R ) 1.0-1.5, unlike the other models. They deliver instead a rather flat evolution, in better agreement with the experimental data. It is in fact not surprising that the Eng and 2E yield similar results, since the principle of both approaches is exactly the same, except for the micromixing rate. The different micromixing rates used in these two models do not have a great effect on the results in this particular case. However, the fact that the DQMOM-IEM(3E) model predictions are quite close to those of two-environment based closures in Case 1 is another interesting point, which could deserve further investigations. Moreover, both in the ME and DQMOM-IEM models, the three-environment based approaches yield results obviously different from those of their two-environment based variants. But there is unfortunately no regular relative behavior of this distinction. Finally, it can be concluded that the computations employing the Eng, 2E, and DQMOM-IEM(3E) models reproduce the experimental data with the best general agreement in both configurations A and B (see Table 2).

Figure 5. Mean crystal size profiles averaged at reactor outlet in Case 2 according to our numerical calculations in comparison with the experimental data.

Figure 5 depicts the results obtained from the simulations performed for Case 2. It might again be useful to remind that the experimental data including the error bars are obtained from the study of Bałdyga and Orciuch.29 As can be seen, all predictions, even without any micromixing approach, are quite close to each other. The influence of the micromixing model on the mean crystal size is always lower than 30%, with the largest influence for the lowest Reynolds number. Moreover, the influence of micromixing for all models decreases as expected when the Reynolds number increases, due to the more efficient macroscale mixing induced by more turbulent hydrodynamic conditions. The agreement is globally quite satisfactory considering the fact that the averaged relative errors are all below 20% (Table 2). The smallest difference between the numerical and experimental results are observed at Re ) 40 000 (Figure 5). The results of all micromixing models are close to the experiments at higher Reynolds numbers. The 3E model predicts the experimental data in a more accurate manner in the lower Re region, but it is difficult to define a globally optimal model. It is nevertheless worth noting that the three-environment models (3E and DQMOM-IEM(3E)) lead to a flatter slope which is closer to that of the experimental data. When summing up the relative differences, the 3E model provides the best general agreement with the experiments (averaged relative error of 11.5%) in Case 2. On a standard PC, the resulting computing time is around 1 h for the flow field computation. The further computing (CPU) times for the reaction simulations (taking into account both mixing and precipitation steps) are listed in Table 2 for all approaches. Although the DQMOM-IEM(3E) model predicts the experimental data with the lowest mean relative error (24.3%) as a whole, the needed computational effort with this approach is the highest (35 h). The mean error level of Eng and 2E are similar (25-26%) but are obtained for a drastically lower CPU time (1 h). Consequently, taking the accuracy of the numerical results in the two cases and the computational costs into account, the 2E model appears to be globally the best compromise for representing the mixing-precipitation process in the investigated configurations. However, it should be kept in mind that this conclusion might be, of course, configurationdependent. Other closures may be more suitable in configurations and/or hydrodynamic conditions different from those

1006 Ind. Eng. Chem. Res., Vol. 48, No. 2, 2009

studied in the present work. It is therefore impossible to decide from such results that a model will be generally more appropriate. 6. Conclusions In this work, three conventional micromixing approaches, namely Eng, ME, and DQMOM-IEM, have been employed in CFD simulations of BaSO4 precipitation reaction in two different single-jet coaxial tubular reactors. The results of these simulations have been compared with those of the calculations performed without micromixing model in order to examine the influence of micromixing on this process. Moreover, a validation of the numerical models is carried out by comparison with reference experimental data obtained from literature. The formulation of the implemented micromixing models is based on two PDFs. One of these is the mixture fraction PDF, defining the reactant concentrations as a function of local mixture fraction values. The second one is the FM-PDF method, which is capable of representing the composition PDF by a finite set of delta functions corresponding to different environments (or modes) within each volume element. The correct local activity coefficient values have been considered in the precipitation kinetics equations (rather than assuming it to be equal to unity), in order to obtain more accurate results. The calculations applying SMOM to solve the PBE have been performed in successive steps in order to reduce the needed computing time. According to the results, it can be concluded that although the influence of the micromixing model on the results is generally clearly visible, it might be negligible for example at higher Reynolds numbers or lower supersaturation ratios. Moreover, the predictions of the employed closures may show similarities with each other and/or may show drastic deviations from the experimental data. In some cases, the computations without any micromixing model lie even closer to the experimental values. Thus, it is not an easy task to come up with a definite optimum approach. Nevertheless, the 2E model represents the best compromise for describing micromixing in the present configurations, when both the overall accuracy of the obtained results and the computational times are taken into account. DQMOM-IEM(3E) is slightly more accurate, but requires 35 times more computing effort. No model leads to a perfect agreement with the experiments for all cases and configurations. It is of course not absolutely clear if the discrepancy results from numerical approximations, inappropriate physical models, or even experimental uncertainties. To check further this issue, the simulations in the tubular reactor will be repeated by taking the influence of aggregation and breakage into account, even if both processes are not expected to be essential in the present case. Even more important might be the implementation of alternative precipitation kinetics33 in the further analysis. Acknowledgment The authors thank D. L. Marchisio and W. Orciuch for helpful discussions during this work. Nomenclature Ci ) reactant concentration, kmol m-3 Ci/ ) reactant concentration at crystal surface, kmol m-3 Ci0 ) feed reactant concentration, kmol m-3 Cφ ) micromixing constant d43 ) volume-averaged mean crystal size, µm fφ ) joint composition PDF G ) growth rate, m s-1

J ) nucleation rate, m-3 s-1 k ) turbulent kinetic energy, m2 s-2 kd ) mass transfer coefficient, (m s-1)(m3 kmol-1) kg ) growth rate constant, m s-1 kv ) crystal volume shape factor Ksp ) solubility product, kmol2 m-6 L ) axial position in reactors, mm n ) environment (or mode) Ne, Ns ) total number of environments and scalars, respectively pn ) probability of environment n r ) radial position in reactors, mm Ri, Ro ) radius of inner and outer pipes, respectively, mm sn ) weighted mixture fraction of environment n S ) supersaturation ratio Spn, Ssn ) source term for probability and weighted mixture fraction of environment n, respectively SY ) source term for chemical reaction progress variable Sµj ) source term for moments, mj m-3 s-1 t ) time, s x ) position vector, m W ) molecular weight, kg kmol-1 Y ) reaction progress variable Greek Symbols R ) scalar in the numerical model and reactant concentration ratio of the feeds ε ) turbulence dissipation rate, m2 s-3 φR ) composition of scalar R γ ) micromixing rate, s-1 γs ) spurious dissipation rate, s-1 γ( ) activity coefficient Γt ) turbulent mass diffusivity coefficient, m2 s-1 µj ) jth moment of the CSD, mj m-3 V ) kinematic viscosity, m2 s-1 F ) density, kg m-3 τm, τr ) time scale for micromixing and reaction, respectively, s 〈ξ〉 ) local mean mixture fraction 〈ξ′2〉 ) mixture fraction variance ξn ) mixture fraction of environment n ξs ) stoichiometric mixture fraction Ψ ) sphericity of particles ψ ) composition vector Dimensionless Numbers Re ) Reynolds number Sct ) turbulent Schmidt number AbbreViations CFD ) computational fluid dynamics conf ) configuration CPU ) central processing unit CSD ) crystal size distribution DQMOM ) direct quadrature method of moments E ) environment EDD ) engulfment-deformation-diffusion Eng ) engulfment Err ) relative error FM ) finite mode IEM ) interaction by exchange with the mean iuf ) initially unmixed fluid ME ) multienvironment mf ) mixed fluid PBE ) population balance equation PDF ) probability density function pf ) pure fluid r ) reaction RANS ) Reynolds-averaged Navier-Stokes

Ind. Eng. Chem. Res., Vol. 48, No. 2, 2009 1007 SMOM ) standard method of moments UDF ) user-defined function UDS ) user-defined scalar

Literature Cited (1) Klein, J.-P.; David, R.; Villermaux, J. Interpretation of Experimental Liquid Phase Micromixing Phenomena in a Continuous Stirred Reactor with Short Residence Times. Ind. Eng. Chem. Fundam. 1980, 19, 373–379. (2) Bałdyga, J.; Bourne, J. R. Turbulent Mixing and Chemical Reactions; John Wiley and Sons: New York, 1999. (3) Vicum, L.; Ottiger, S.; Mazzotti, M.; Makowski, L.; Bałdyga, J. Multi-scale Modeling of a Reactive Mixing Process in a Semibatch Stirred Tank. Chem. Eng. Sci. 2004, 59, 1767–1781. (4) Villermaux, J.; Devillon, J. C. Repre´sentation de la coalescence et de la redispersion des domaines de se´gre´gation dans un fluide par un mode`Ie d’interaction phe´nome´nologique In Proceedings of the Fifth European (Second International) Symposium on Chemical Reaction Engineering, Amsterdam, Netherlands, 1972; Paper B1-13. (5) Fox, R. O. On the Relationship between Lagrangian Micromixing Models and Computational Fluid Dynamics. Chem. Eng. Process 1998, 37, 521–535. (6) Fox, R. O. Computational models for turbulent reacting flows; Cambridge University Press: Cambridge, 2003. (7) Wang, L.; Fox, R. O. Comparison of Micromixing Models for CFD Simulation of Nanoparticles Formulation. AIChE J. 2004, 50, 2217–2232. (8) Bałdyga, J.; Orciuch, W. Closure Problem for Precipitation. Chem. Eng. Res. Des. 1997, 75, 160–170. ¨ ncu¨l, A. A.; Sundmacher, K.; The´venin, D. Numerical Investigation (9) O of the Influence of the Activity Coefficient on Barium Sulphate Crystallization. Chem. Eng. Sci. 2005, 60, 5395–5405. (10) Nishimura, Y.; Matsubara, M. Micromixing Theory via the TwoEnvironmental Model. Chem. Eng. Sci. 1970, 25, 1785–1797. (11) Tsai, K.; Gillis, P. A.; Sen, S.; Fox, R. O. A Finite-Mode PDF Model for Turbulent Reacting Flows. J. Fluids Eng.: Trans. ASME 2002, 124, 102–107. (12) Taguchi, K.; Garside, J.; Tavare, N. Mixing, Reaction and Precipitation: Semibatch Barium Sulphate Precipitation. Inst. Chem. Eng. Symp. Ser. 1999, 146, 395–419. (13) Marchisio, D. L.; Barresi, A. A. CFD Simulation of Mixing and Reaction: The Relevance of the Micro-Mixing Model. Chem. Eng. Sci. 2003, 58, 3579–3587. (14) Villermaux, J.; Falk, L. A Generalized Mixing Model for Initial Contacting of Reactive Fluids. Chem. Eng. Sci. 1994, 49, 5127–5140. (15) Piton, D.; Fox, R. O.; Marcant, B. Simulation of Fine Particle Formation by Precipitation Using Computational Fluid Dynamics. Can. J. Chem. Eng. 2000, 78, 983–993. (16) Bałdyga, J.; Bourne, J. R.; Hearn, S. J. Interaction between Chemical Reactions and Mixing on Various Scales. Chem. Eng. Sci. 1997, 52, 457– 466. (17) Marchisio, D. L. Precipitation in Turbulent Fluids. Ph.D. Dissertation, Politecnico di Torino, Italy, 2002.

¨ ncu¨l, A. A.; The´venin, D. Techniques for (18) John, V.; Angelov, I.; O the Reconstruction of a Distribution from its Moments. Chem. Eng. Sci. 2007, 62, 2890–2904. (19) Fox, R. O. Computational Methods for Turbulent Reacting Flows in the Chemical Process Industry. ReV. Inst. Fr. Pet. 1996, 51, 215–243. (20) Harada, M.; Arima, K.; Eguchi, W.; Nagata, S. Micro-mixing in a Continuous Flow Reactor (Coalescence and Redispersion Model). Memoirs Facul. Eng.-Kyoto UniV. 1962, 24, 431–446. (21) Costa, P.; Trevissoi, C. Some Kinetic and Thermodynamic Features of Reactions Between Partially Segregated Fluids. Chem. Eng. Sci. 1972, 27, 653–668. (22) Costa, P.; Trevissoi, C. Reactions with Non-linear Kinetics in Partially Segregated Fluids. Chem. Eng. Sci. 1972, 27, 2041–2054. (23) Marchisio, D. L.; Fox, R. O. Solution of Population Balance Equations Using the Direct Quadrature Method of Moments. J. Aerosol Sci. 2005, 36, 43–73. (24) Pope, S. B. PDF Methods for Turbulent Reactive Flows. Prog. Energy Combust. Sci. 1985, 11, 119–192. (25) Marchisio, D. L.; Barresi, A. A.; Fox, R. O. Simulation of Turbulent Precipitation in a Semi-Batch Taylor-Couette Reactor Using CFD. AIChE J. 2001, 47, 664–676. (26) Gavi, E.; Rivautella, L.; Marchisio, D. L.; Vanni, M.; Barresi, A. A.; Baldi, G. CFD Modelling of Nano-Particle Precipitation in Confined Impinging Jet Reactors. Chem. Eng. Res. Des. 2007, 85, 735–744. (27) Marchisio, D. L.; Barresi, A. A.; Garbero, M. Nucleation, Growth, and Agglomeration in Barium Sulfate Turbulent Precipitation. AIChE J. 2002, 48, 2039–2050. ¨ ncu¨l, A. A.; The´venin, D.; Sundmacher, K. Numerical Simulation (28) O of BaSO4 Precipitation in a Coaxial Pipe Mixer with Micromixing Effects. Presented at the AIChE Annual Meeting, San Francisco, CA, November 2006; Paper 70c. (29) Bałdyga, J.; Orciuch, W. Barium Sulphate Precipitation in a PipeAn Experimental Study and CFD Modelling. Chem. Eng. Sci. 2001, 56, 2435–2444. (30) Nielsen, A. E. Nucleation and Growth of Crystals at High Supersaturation. Krist. Tech. 1969, 4, 17–38. (31) Bromley, L. A. Thermodynamic Properties of Strong Electrolytes in Aqueous Solutions. AIChE J. 1973, 19, 313–320. (32) Wei, H.; Garside, J. Application of CFD Modelling to Precipitation Systems. Trans. Inst. Chem. Eng. 1997, 75A, 219–227. (33) Mersmann, A. Crystallization Technology Handbook, 2nd Ed.; Marcel Dekker: New York, 2000. ¨ ncu¨l, A. A.; Sundmacher, K.; Seidel-Morgenstern, A.; The´venin, (34) O D. Numerical and Analytical Investigation of Barium Sulphate Crystallization. Chem. Eng. Sci. 2006, 61, 652–664.

ReceiVed for reView March 5, 2008 ReVised manuscript receiVed October 16, 2008 Accepted October 22, 2008 IE800364K