Reactive mixing in a tubular jet reactor: a comparison of PDF simulations with experimental data

Chemical

Engineering

Science, Vol. 49, No. 248, pp. S229mS241, 1994 Copyright fQ 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved OW-2509194 57.C0+O.CQ

OCKW-2509(94)00299-l

REACTIVE MIXING IN A TUBULAR JET REACTOR: A COMPARISON OF PDF SIMULATIONS WITH EXPERIMENTAL DATA MASSIMO Dipartimento

PIPINO

di Scienza de1 Materiali e Ingegneria Chimica, Politecnico di Torino, C. so Duca degli Abruzzi 24, 10124 Torino, Italy

and RODNEY

0.

FOX*

College of Engineering, Kansas State University, Manhattan, (Accepted for publication

26 September

KS 66506, U.S.A.

1994)

Abstract-The simulation of turbulent reactive flows in chemical reactors can be extremely useful for practical applications. In particular, both industrial and academic chemical reaction engineers are keenly interested in the exact description of the composition field in a reactor volume. Results from two different probability density function (PDF) descriptions of a laboratory tubular jet reactor are compared with experimental data for acid-base neutralization in turbulent liquid media. The first description is based on the Lagrangian joint PDF of velocity and composition, while the second is based on the Eulerian composition PDF. In both descriptions chemical reactions are treated exactly, but molecular mixing must be modeled. A computational fluid dynamics code (FLUENT) provides the mean velocity field and turbulence quantities. Results are found to be sensitive both to turbulent diffusivity and to the model for molecular mixing. A new model for the effect of scalar integral-length scale relaxation on the molecular mixing rate is proposed. By comparison with experimental data, the superior performance of a new model over existing molecular mixing models is demonstrated.

1. INTRODUCTION

The simulation of turbulent reactive flows in real reactors can be extremely useful for practical applications. In particular, both industrial and academic chemical engineering researchers are keenly interested in predicting the composition fields in chemical reactors. The conventional deterministic approach to this problem based on solving moment equations for the composition variables is of limited use due to the well-known difficulties in closing the nonlinear reaction rate terms in turbulent flows. The probability density function (PDF) approach overcomes these difficulties by describing each point in the turbulent flow by an n-dimensional joint PDF of random variables such as the velocity and/or the chemical compositions. Depending on the level of closure, the joint PDF may include the fluctuating velocity components, thereby avoiding the closure problem (i.e. the gradient-diffusion assumption) associated with the composition-velocity fluctuation transport term in the composition moment balance equation. However, the most important advantage of the PDF approach for reactive flow applications is the exact treatment of the reaction term. PDF methods are thus a promising tool for chemical reaction engineering applications. *To whom correspondence should be addressed.

This work investigates two different PDF descriptions of a laboratory tubular reactor, comparing the predictions to experimental data for acid-base neutralization: A+BdP

(1)

with k, = 1.0 x 10’ ma/km01 s. The first description is based on the joint velocity-composition PDF transport equation for Lagrangian fluid particles traveling through the reactor (Pope, 1985). Each notional particle in the resulting Monte-Carlo simulation carries with it random variables specifying its location, velocity, and chemical composition. The turbulent motion is modeled by a stochastic differential equation whose coefficients are functions of the local mean velocity, the local turbulent kinetic energy, and the local turbulent energy dissipation (Pope, 1985; Tsai and Fox, 1994). Unlike in standard moment closures, the transport of chemical species due to velocity fluctuations is treated exactly at this level of description. However, the molecular mixing term must be closed, using a local characteristic mixing time found from the turbulence fields. As shown in the sequel, the proper specification of the molecular mixing rate is crucial for successfully modeling the experimental data. The second description is based on the transport equation of the Eulerian joint composition PDF. For

5229

MASSIMO PIPINO and RODNEY 0. Fox

5230

I

HCI

NaOH

Fig. 1. Experimental setup for tubular jet reactor this case, the velocity does not appear as a random

variable in the joint PDF, and thus a computational fluid dynamics code (FLUENT) is required to determine the mean velocity and the turbulent diffusivity. A novel Monte-Car10 solution algorithm has been developed that interfaces directly with FLUENT. The description, although in theory less rigorous than the previous one due to the need to model transport by turbulent velocity fluctuations, is nevertheless very promising for industrial applications since the lower dimensionality of the PDF allows for full 3D simulations or, in 2D cases, requires less memory and CPU time. In the present code, a gradient-diffusion approximation is employed. Using the PDF codes, simulations have been carried out to predict the composition fields in a laboratory tubular reactor with nonpremixed coaxial feeds. The extent of the reaction zone is determined experimentally by the use of a chemical indicator (phenol red) added to both feed streams. The axial concentration profile for the key reactant (fed in the inner tube) has been determined by a spectrophotometer with fiber-optic sensors able to move along the axis of the reactor in correspondence to the diametrical plane. Data have been collected for various Reynolds numbers, various inlet concentration ratios of the two reactants, and several internal tube geometries, and compared to simulation results from the two PDF codes. The remainder of this work is organized into five sections. Section 2 describes the experimental setup. Section 3 reviews the PDF approach and the models that close the PDF transport equation. A new model is introduced for the rate of molecular mixing which, unlike previous formulations, explicitly accounts for the relaxation of the concentration spectrum. Section 4 briefly describes the numerical procedures employed in the PDF simulations. Results from the PDF codes are compared with the experimental data in Section 5. Section 6 draws conclusions and discusses some possible extensions.

2. EXPERIMENTAL

SETUP

In order to check the results given by the simulations, axial concentration profiles have been measured experimentally in a tubular jet reactor with nonpremixed coaxial feed. This choice is due to the simple hydrodynamic field involved in such a geometry and to the fact that this configuration allows simple measurements of concentration decay. The experimental setup (see Fig. 1) was realized in the micromixing laboratory of Politecnico di Torino with the use of a silica tube (10 mm interval diameter, 1 m long) mounted on an optical bench, in which different coaxial feed tubes (0.7-1.2 mm internal diameter) could be inserted. Note that the outer-to-inner tube diameter ratio of approximately 10 implies that both radial mixing due to turbulent diffusivity and turbulent molecular mixing will strongly affect the measured profiles. The measurement of axial concentration profiles was carried out using a fiber-optic spectrophotometer, able to record absorbance along the reactor axis. Both nonreactive and reactive experiments were carried out. Results from the former were used to check the rate of radial mixing, and from the latter to test the molecular mixing models. In the nonreactive case, methylene blue was used as an inert tracer, and the axial absorbance profile was measured using visible light (658 nm). For the reactive cases, NaOH-HCI neutralization was employed, and the residual base concentration visualized by relating it to the absorbance of phenol red measured at 558 nm. With the present experimental setup, soda (001 N) could be fed into the internal stream at various velocity ratios (mean internal velocity/mean external velocity) ranging from l-10 and various concentration ratios (acid concentration: base concentration) ranging from 0.5-20. In order to allow the comparison of absorbance data with the “absorbance” field produced by the simulation codes, the linearity of the Beer-Lambert

5231

Reactive mixing in a tubular jet reactor law in the experimental concentration range was verified, together with the influence of the local pH on the change of color of the indicator. In the simulation codes, the local absorbance field was computed from the local mixture fraction PDF using the method described by Li and Toor (1986). 3. PDF

METHODS

Excellent introductions to PDF methods can be found in O’Brien (1980) and Pope (1985), and a summary of their application to chemical reaction engineering problems can be found in Tsai and Fox (1994) where the PDF transport equation is derived and models for closing the velocity-composition joint PDF are reviewed. Here, discussion will be limited to a brief review of the relevant PDF transport equations and the closure models employed in the PDF codes. In particular, the molecular mixing model will be discussed in detail due to its strong influence on the predicted axial concentration profiles for the reactive cases.

appears in a similar form in the composition joint PDF transport equation and is discussed in detail below. 3.2. Composition joint PDF The Eulerian composition joint PDF approach is based on the solution of the PDF transport equation for &(JI; x, t):

3.1. Velocity-composition joint PDF The transport equation for the velocity-composition joint PDF f(u, C; x, t) can be written for a constant density flow as

The terms on the left-hand side of eq. (6) represent, respectively, the change in time of the PDF, convective transport in physical space due to the mean velocity field and convective transport in composition space due to chemical reactions. These terms appear in closed form. On the right-hand side of eq. (6), diffusive transport in composition space due to molecular mixing, and transport in physical space due to turbulent velocity fluctuations, can be identified. They appear in terms of conditional expectations and must be modeled. The modeling of the molecular mixing terms is common to both the composition joint PDF and the velocitycomposition joint PDF approach and is discussed below. Compared to the velocity-composition PDF approach, the composition PDF presents the disadvantage of modeling turbulent transport in physical space. In this work, turbulent transport is assumed to follow the gradient-diffusion model:

where the scalar flux JF and the viscous force zij have the following forms:

(7) where I-r represents the turbulent diffusivity. When the k-e model is used to describe the turbulent motion field, we can express the turbulent diffusivity as

In eq. (2), the mean pressure ( p) can be obtained from the Poisson equation: V’(P)

=

d2(UiUj> --P

Z’

(5)

A remarkable feature of the PDF approach can be seen by observing that all the terms on the left-hand side of eq. (2) are closed, which include the turbulent convection, the production by the mean pressure gradient, gravity and chemical reactions. However, the three conditional expectation terms on the right-hand side, which represent transport in velocity space by viscous forces, and the fluctuating pressure gradient and transport in composition space by diffusive molecular mixing, must be closed. The first two of these terms involve models that are independent of the composition fields and are discussed by Pope (1985) and Tsai and Fox (1994). The final term

OTC

where C,, is an empirical constant (set here to 0.09) and bT is the turbulent Schmidt number. The standard value of bT is 0.7; however, comparisons of the nonreactive data with simulation results suggest that this value is too small (radial mixing is too fast). The gradient-diffusion assumption can fail in variable-density reactive flows, but it is accurate in representing many constant density flows. The effect of varying CJ=has been found to be significant for the tubular jet reactor used in this work. Good agreement with experimental results for the nonreactive case has been found using CT== 1.0. This value has thus been employed in all isokinetic jet simulations using the composition PDF code reported in this work. For the non-isokinetic jet simulations, a low Reynolds number correction (M. Pipino, 1994, personal communication) shown in Fig. 2, has been used in some cases. In the

5232

MASSIMO PIPINO and RODNEY 0.

Fox

constant density with no reaction gives:

we see that the molecular mixing term is modeled as

For nonhomogeneous conditions, the same model is used but with Q and (4.) replaced by their local values. Fig. 2. Low Reynolds number correction factor for turbulent diffusivity employed in nonisokineticjet simulations.

joint velocity-composition code br = 1.0 corresponds to a value of C,, in the range of 8-9 [standard value is 2.1 (S. B. Pope, 1993, personal communication)]. All simulations using the joint velocity-composition code reported in this work have been found with ca = 8. 3.3. Molecular mixing A successful model for closing the molecular mixing term must address two problems: first, a passive scalar PDF should relax to a Gaussian form regardless of its initial condition; second, the effect of turbulent stretching on the rate of scalar mixing should be properly accounted for. Several models have been proposed to address the first problem (Curl, 1963; Dopazo, 1975; Janicka et al., 1979; Pope, 1982, Gao, 1991; Valifio and Dopazo, 1991; Fox, 1992, 1994; Norris, 1993; Fox and Tsai, 1994). In this work, the LMSE (linear mean square error) mode1 will be employed. This model is exact when the initial PDF is Gaussian (O’Brien, 1980) and has the simplest mathematical form. For the simple acid-base reaction under consideration, the results are insensitive to the PDF shape; hence, the problem of turbulent stretching is more significant. In this section, a new model for the molecular mixing rate is proposed to explicitly account for relaxation of the concentration spectrum to its isotropic form through the actions of turbulent mixing. 3.3.1. LMSE mixing model. The LMSE model in isotropic, homogeneous turbulence with an isotropic scalar field reads

!&_at

-;

c&I - <&>I,

where r,+is the molecular mixing time scale determined by turbulent stretching. The corresponding PDF evolution equation can be written as

By comparing eq. (10) with the evolution equation of &, [Eq. (6)] under conditions of homogeneity and

3.3.2. Molecular mixing rate. For modeling the effect of molecular mixing on nonpremixed chemical reactions, the choice of the molecular mixing rate is paramount. Several models have been suggested in the literature to relate z+ to local turbulence quantities such as k and E. Baldyga and Bourne (1989) proposed a simplified micromixing model that incorporates the “hydrodynamic lifetime of a vortex” that they compute to be proportional to (v/s)‘/‘, e.g. the Kolmogorov time scale (Baldyga and Bourne, 1984). Using their model, the molecular mixing time can be written as 54.1 = 17.31

0 1

E

l/2

.

(13)

Due to the choice of turbulence time scale, this model is applicable when the rate of mixing is controlled by vortex engulfment at fine scales. The standard model for high Reynolds number isotropic turbulent flows relates the molecular mixing time to the eddy turnover time which is proportional to k/E (Pope, 1985): 2k

(14)

Q.,2 = ,-,

where C, has a limiting value of 2 for high Reynolds number flows. This choice of r+ (at least in the asymptotic decay regime*) is supported by direct numerical simulations of a nonpremixed scalar field in isotropic, homogeneous turbulence (Eswaran and Pope, 1988). Use of r+,s is appropriate when scalar gradients are created by integral-scale eddies and the rate of transport from integral scales to small scales is the rate-controlling step. Use of r’o,2 also assumes that the scalar spectrum is fully developed at all scales. Equations (13) and (14) differ by a factor proportional to the turbulence Reynolds number defined by (Hinze, 1975): 60 Re,=-- 43&

k

Thus, as the turbulence Reynolds number increases, T+.i will decrease much faster than T&.~. If scalar

*That below.

is, at the end of the relaxation

period modeled

Reactive mixing in a tubular jet reactor

gradients are initially created at scales near the integral scale of the turbulence field as is the case for turbulent jets, at high Reynolds number, r9,2 will be the appropriate limiting molecular mixing time. On the other hand, if the scalar gradients are created on smaller scales, as performed by Baldyga et al. (1994) by using a separate feed tube, r+,t may become the limiting molecular mixing time. The spectral relaxation model presented below explicitly accounts for different initial scalar integral-length scales by introducing a cascade of scalar dissipation from integral scales, through the Kolmogorov scale, down to the Batchelor scale. Corrsin’s model (Corrsin, 1964), which combines r+. I and rb. z was also tested:* 70.3

for SC >> 1. Note number, eq. (16)

=

i

i

+

5

ln(Sc)

Y 0&

5233

3.3.3. Spectral relaxation model. A scalar field advected by an incompressible turbulent flow is governed by

W - = n-Q,

(17)

Dt

and the magnitude of its gradient by (see Appendix)

where xi = @/axi and eij = &,/axi is the strain rate tensor (repeated indices imply summation). Reynolds averaging of eq. (17) multiplied by 4’ yields (see Appendix)

r/2

(19)

(16)

that at high turbulence Reynolds approaches eq. (14) if C, = 4/3.’

Note also that Corrsin assumed that the scalar spectrum was fully developed when deriving t+, 3. Experience, in this work, with the above models suggests that none are entirely adequate for a tubular jet reactor. For example, C, in eq. (14) must be set smaller than the limiting value in order to achieve even qualitatively satisfactory agreement with experimental data. In confined flows, where reaction occurs near the bounding wall, E attains high values in the boundary layer (i.e. Re, decreases) resulting in an overprediction of the mixing rate (Tsai and Fox, 1994). Furthermore, even for free jets, all three models may overpredict mixing rates near the injection point due to the lack of initial integral-length scale information and the assumption that the scalar spectrum is fully deve1oped.t Indeed, even in stationary isotropic turbulence, if the initial integral-length scale of the scalar field is larger than the integral-length scale of the velocity field, rg approaches k/c only asymptotically from above (Eswaran and Pope, 1988; Kosaly, 1989; Jiang and O’Brien, 1991).6 These observations suggest that an improved molecular mixing model is needed that includes the effect of Re, and incorporates integral-length scale information in the form of a spectral relaxation period where r+ moves to its limiting value.

*The integral-length scale of the velocity field has been taken to be proportional to k3/‘/c (Hinze, 1975) and the proportionality constant fixed by forcing the original model [e?. (69) in C orrsin (1964)] to agree with ro.s at SC = 1. This differs from C, = 2 (valid for SC = 1.0) due to a dependence on SC (Corrsin, 1964). :In pure mixing studies, this problem has been avoided in the past by adjusting the data to eliminate the initial decay region at the injector (e.g. Brodkey, 1966). However, for reacting flows this region cannot be ignored since a considerable percentage of fast reactions can occur there (Baldyga et al., 1994). “These studies have also demonstrated that the initial scalar-to-velocity micro-length scale ratio is unimportant in determining scalar decay rates.

where the source term can be expressed as $2 Likewise, Reynolds Appendix):

=

a<4> --. axi averaging

of eq. (18)

(20) yields (see

The scalar dissipation can be defined as E, = f-xf2, and the molecular mixing rate as

such that r+ = l,l. Note that both E+ and r4 are random fields since they involve the fluctuating variable xJ2. side of eq. (21) The second term on the right-hand corresponds to turbulent stretching of the scalar gradient, and the first term to its molecular dissipation. The last two terms are production terms: the first due to the mean scalar gradient and the second to the mean velocity gradient. Fox (1994), based on studies of diffusion in random scalar fields and direct numerical simulation of scalar mixing in isotropic stationary turbulence, proposed the following models for the molecular dissipation and turbulent stretching terms:

and

r

= -G

0E “2 ;

(251

where C, = 0.54 (Fox et al., 1992; Leonard, 1989) is near Batchelor’s constant (Batchelor, 1952), and C, = 3 for one-dimensional diffusion in random

MAWMO PIPINOand RODNEY 0.

5234

Fox

2 E l/Z 3(azt3 Re, - 1) 0 Y
lamellar systems (Fox, 1994). Finally, if (xix>> is assumed to be nearly isotropic, the source term due to the mean velocity gradient can be neglected for incompressible flows.* Following Newman et al. (1981), the source term due to the mean scalar gradient can be modeled by

(29) 2

E 112

3(a 2’3 Rr, - 1) 0 ; s 112 2 <%) ln(Sc) (-1v


or (27)

0 E

where C, in the spectral relaxation model is uniquely determined in terms of C, and C, as discussed below. Note that the source term for (S,,/T) is proportional to .+ divided by the square of the Batchelor length scale, and thus represents a source for the maximum attainable scalar dissipation (i.e. with all scalar variance located at dissipative scales.) The source term differs from that employed by Mantel and Borghi (1991) where the Taylor length scale is used in place of the Batchelor length scale in order to obtain the correct limiting behaviour for . The spectral relaxation model describes the cascade of scalar dissipation from large to small scales by introducing additional variables that correspond to “potential” scalar dissipation residing at length scales larger than the Batchelor length scale. These variables are referred to as potential scalar dissipation because they are formed by dividing the scalar variance in a given range of length scales by the square of the Batchelor length scale. Thus they represent the maximum attainable scalar dissipation for a given scalar spectrum. The cascade model we propose is similar in spirit to the three-stage model proposed by Baldyga (1989) for the scalar variance. We denote the potential scalar dissipation ith stage of the cascade by E$. The cascade

in the

will be composed of up to three stages representing transport from the integral scale of the turbulence field, L,, to the initial7 integral scale of the scalar field, L,; transport from L, to the Kolmogorov scale of the turbulence field; and transport from the Kolmogorov scale to the Batchelor scale (Corrsin, 1964). Based on the form of eq. (21), we define the three-stage cascade model for SC r

1 and ReF3”

< c( = 2Lg/3L,

I

1 bys

(28)

*In Mantel and Borghi (1991), the source term is modeled by assuming it to be proportional to (uiu;> which is null in isotropic turbulence. +The “usual” initial condition is L+ = 2L,/3 (Corrsin, 1964). However, in some cases the scalar may be generated by a different mechanism than the turbulence field. For example, in direct numerical simulations the fields may be initialized independently of one another. Experimentally, L, << L, can be realized by using a separate feed tube for key reactants (Baldyga et al., 1994).

x +Cw where
Re,

=




=k. fi

-

V

l/Z


(33)

Note that eqs (19) and (23) are used to relate the scalar variance to the (E+). Comparison of the above model with direct numerical simulation data at Re, = 49, SC = 0.7 (two-stage model) and experimental data at Re, = 35, SC = 450, yield good agreement as discussed below. The source term is placed at the velocity integral scale since it describes gradient production due to the product of the scalar flux and the mean scalar gradient, both of which are integral-length scale terms. For the case where all scalar dissipation is produced by the source term (all dissipation variables are initially zero), a = 1 and stage 1 can be dropped. A further refinement of the model could include a dynamical equation for a in order to handle nonsteady flows. The dependence of the decay rates on Re, and SC follow from Corrsin (1964) and essentially represent the time needed for a peak in the scalar spectrum at one stage to arrive at the next stage. For small values of Re, and SC, these decay rates are large compared to the scalar dissipation rate, and stages 2 and 3 have little influence on the behavior of <.Q >. For this case, a one-stage model wherein eq. (29) decays directly into eq. (31) can be employed. For our reactor (Re, z 35, SC = 450), a one-stage model found by dropping eqns (28) and (30), taking a = 1 and Re,, SC >> 1 in eq. (29), and C, = 0.18 in eq. (31) gives good agreement with the full model.

*For 0 < I, the integral scale of the scalar field is above the integral scale of the turbulence field and stage 1 is dropped. For c( 5 ReF3” 9the integral scale of the scalar field is below the Kolmogorov scale.

Reactive mixing in a tubular jet reactor

Three first three terms on the right-hand side of the balance equation for [eq. (31)] model, respectively, the cascade (source) from larger scales, molecular dissipation, and turbulent stretching in the viscousdiffusive subrange of the scalar spectrum. The latter is taken directly from the model proposed by Fox (1994) [eq. (25)]. The characteristic dissipation rate in the molecular dissipation term has been modified [cf. eq. (24)] to include the sum of both the actual and potential scalar dissipation rates. This modeling assumption is justified for two reasons: (1) it yields the correct limiting behavior as Re, + 00, and (2) C, can be taken as independent of the initial scalar-tovelocity integral-length scale ratio and of Re,. The dependence on SC has been added to agree with Corrsin (1964) in the constant mean gradient case. Limiting cases are discussed below. Direct numerical simulation results (Fox et al., 1992) indicate that C, is independent of the initial scalar-to-velocity integrallength scale ratio, but that C, in eq. (24) exhibits a strong dependence on it. When this ratio is greater than unity, C, is initially very large and relaxes slowly to its final value. By including the source cascade (source) term and by modifying the molecular dissipation term, C, = 3, independent of SC, Re, and the initial scalar-to-velocity integral-length scale ratio.

5235

Fig. 3. The mechanical-to-scalar time-scale ratio, R, for isotropic, stationary turbulence (Re, = 49, SC = 0.7) found from the spectral relaxation model with t* = et/k. Solid line: (a = 0.989) = 0. Dashed line: (u = 0.456) = 8.25E/k, (~(0)) = 0.25&/k. Dash-dot line: (a = 0.300) = 0.75&/k. Asterisks: (a = 0.168) = 0.67&/k, = 11.25&/k, = 1.756/k. Dotted line: (a = 0.121) = 0.67&/k, = 1.0&/k. Crosses: (K = 0.121) = 0.67&/k, = 2.75&/k.

actual scalar dissipation quickly becomes constant: 3.3.4. Limiting cases. The limiting case of homogeneous turbulence and scalar fields with S,, constant (constant mean scalar gradient) can be used to relate C, to C, and C,. When Re, is very large, the scalar dissipation rate is controlled by the rate of transport from large to small scales so that (rd} + E/k (SC = 1). For SC >> 1, the limiting behavior has been determined by Corrsin (1964). Setting the material derivatives equal to zero in eqs (19) and (28)~(31), letting G(= 1, and taking the limit Re, + m, we find that c,=-.

CL0 c, - 1

(34)

For finite Re,, t,+ = k/e + 1.23(v/e)“’ when SC = 1. The limiting case of isotropic, stationary turbulence with no mean scalar gradient (S,l = 0) is also of interest. For this case, the material derivatives in eqs (29)-(31) become total derivatives and the behavior of _* In particular, at moderate to large Reynolds numbers, if > 0 the limiting value of (r+) is E/k (see Fig. 3) as seen in direct numerical simulations (DNS) (Eswaran and Pope, 1988). In the other extreme, if = 0 and {r&O)} > 0, corresponding to an initial scalar spectrum with all its energy near the Batchelor length scale, then has a limiting value proportional to (s/v) I/*. Note, however, that regardless of the inuial conditions the sum of potential plus *The expressionsfor and can be found from their definition and eqns (19). (28)-(3 1) after replacing the material derivatives by total derivatives.

+

C,(2 + SC-‘) 2C,-2-Sc-’

s l/Z 0;

C,(2 + SC- ‘1 I<& = 2c, - 2 -SC-’ n;

(35)

where 1, is the Batchelor length scale. This result is consistent with our definition of potential scalar dissipation as being proportional to the scalar variance in a range of length scales divided by the square of the Batchelor length scale. Moreover, it implies that the characteristic length scale for scalargradient dissipation is proportional to A,,. A qualitative comparison of the model predictions with the DNS results of Eswaran and Pope (1988) at ReA = 49 (SC = 0.7) can be made. Since SC -=z1 for this case, we will drop eq. (30) and let stage 2 cascade directly into eq. (31). For this Reynolds number, the maximum mechanical-to-scalar time-scale ratio [R = 2k for other combinations of initial conditions will depend on the initial scalar-to-velocity integral-length scale ratio, i.e. the value of u. We will fix the initial conditions so that eq. (35) is satisfied for all t values. Results are presented in Fig. 3, found using the same values for 0: as Eswaran and Pope (1988), and in general show good agreement with their resuIts. The insensitivity of R to the initial scalar micro-length scale,* i.e. to , (Eswaran and Pope, 1968; Jiang

*Recall that scalar micro-length scale is proportional JCG&

to

MASSIMO PIPINO and

5236

and O’Brien, 1991) seen in Fig. 3.

is also predicted

by the model

as

3.35 Stochastic spectral relaxation model. Equations (19), (28)-(31) and the equation for the scalar mean field could be solved using a standard finitevolume code after invoking the gradient-diffusion assumption to close the turbulent flux terms. However, we will instead propose a stochastic version formulated in terms of a set of linear equations that model the values of Q--E+ following a fluid particle (cf. Pope, 1985). These equations will be solved along with a modified version of the LMSE model. Let E: denote a random variable following a notional particle with mean value = where x* is the location of the notional particle at time t. The stochastic spectral relaxation model can then be expressed as la&: 1 a.52 2 at

E*

(36)

E,,

at

2

k E 112 Et,

2

E

SC*+ i CT - 3(&3 J&, _ l) 0 ;

1 as: 2 dt

ST

(38)

and

1 a&x 2

(37)

at

2 In(Sc)

l/2

4) E

Y

ef -

2c x 2+sc-’

x .$ + C,

RODNEY

0.

Fox

4 (e.g. mixture fraction), eq. (9) must be modified as fo1lows:

-

at

= -lJrf(A

where

p =

(41)

‘<9”> ’= -
(42)

(I

and

that in the derivation of the mode1 for scalar dissipation we have assumed that the scalar is nonreactive, e.g. the mixture fraction in nonpremixed cases. Furthermore, we have assumed that the scalar dissipation rate for reactive scalars is the same as that for the nonreactive scalar so that yr,$ is the same for all scalars. In our PDF codes, eqs (36)-(42) are solved* for each notional particle to find the PDF of the mixture fraction 4. Since the reaction rate is large, this information suffices to compute the absorbance field (Li and Toor, 1986). The inlet value of the mixture fraction is 0 in the outer stream and 1 in the inner jet. The inlet values for E:~E$ are equal and null and a = 1, reflecting the nonpremixed nature of the inlet streams and the fact that the gradients in the scalar field are generated by the same mechanism that generates the turbulence field. Note

4. NUMERICAL

E l/2 et. 0Y (39)

wherein the source terms, mean values, and turbulence fields are evaluated at the spatial location of the notional particle. Note that at this level of closure the model contains no random noise terms and is thus similar in nature to the LMSE model. With the exception of the source terms that couple the equations to the mean scalar gradient (fluxes), the right-hand sides of eqs (36)-(39) are similar in nature to the reaction rate functions appearing in the balance equation for reacting scalars and can be treated analogously. Also note that since r$ = $/<#‘> is a random variable, the molecular mixing rates in notional particles at the same spatial location will differ, depending on the past history of each individual notional particle. This would not be the case if the Reynolds averaged version of the model is employed. In this manner, the stochastic formulation explicitly accounts for the “age” distribution of fluid elements and its effect on molecular mixing. 3.3.6. Modified LMSE model. Since r$ is a random variable, in order to insure that the mean scalar concentration remains constant for an arbitrary scalar 4, and that eq. (19) is satisfied during the molecular mixing step (i.e. with S,, = 0) for a nonreactive scalar

- j?.).

SIMULATION

4.1. Simulation of the jlow field The geometry of the experimental apparatus consists of a tube with a 0.01 m inner diameter and 1 m long, into which a smaller coaxial tube (1.17 mm inner diameter, 0.11 mm wall thickness) is inserted. The experimental measurements show that in the reactive cases the reaction goes to completion in the first 0.1 m. The flow field has been simulated for both isokinetic and nonisokinetic conditions using a total length of 0.4 m. The k--E mode1 in FLUENT has been employed for the isokinetic case and the Reynolds stress model for the nonisokinetic case. At the end of the computational domain all the radial concentration gradients for an inert tracer are essentially gone. The inlet conditions are obtained from a pre-simulation of fully developed (annular) pipe flow. It should be noted that the inner stream is laminar when low vetocity ratios are considered, while inside the reactor both of the streams are turbulent. The axisymmetric, 2D simulation domain is represented by 20 x 200 cells with nonuniform spacing, with a particularly fine grid close to the walls and in the first section of the reactor. In the whole domain the yt values are set greater than 25. The highest value of Re, for this inhomogeneous *In the one-stage cascade model employed in most simulations, $7 and E: are dropped and E: decays directly into of, and C, is set to 0.18 as noted earlier.

Reactive mixing in a tubular jet reactor flow is of the order of 35; thus, spectral relaxation effects are significant as seen below. 4.2. PDF simulations 4.2.1. Lagrangian velocity-composition PDF code. The velocity-composition PDF code (Pope, 1991) is initialized by reading the steady-state solution of the flow field at each grid node from the output of FLUENT, including mean velocity , k, E, and mean static pressure p. The averaged values of those quantities are linearly interpolated onto each notional particle according to its location. In order to have uniformly distributed notional particles in each cell, a specific weight, which is proportional to the radial distance from the symmetry line, is assigned to each notional particle. This is only necessary when a non-Cartesian coordinate system is adopted (a cylindrical coordinate system is used in this study). At each time step, new notional particles are initialized from the inlet with proper weights, flow-field properties and scalar values. The number of notional particles in each cell is controlled by separating heavy particles into lighter particles, or clustering lighter particles into heavier particles according to the average particle weight in each cell. Further details concerning the Monte- Carlo algorithm can be found in Tsai and Fox (1994) and Pope (1985). 4.2.2. Eulerian composition PDF code. A translator code reads output files from FLUENT and prepares the input files for the composition PDF code. The latter can adopt (partially or entirely) generic 2D FLUENT simulation files by adopting the same grid, be it regular or body-fitted. The translator code gives the PDF code the grid information and the cell types, together with the turbulence quantities k and E. A standard option in FLUENT outputs the values of the cell-face flow rates instead of the mean velocities: this is employed with body-fitted coordinates when the computational domain is not aligned with the real geometry. In all the runs of the PDF codes, a reduced domain of 20 x 61 cells has been used, corresponding to the first 0.1 m of the reactor. The Monte-Carlo solver expresses the composition PDF in each cell as an ensemble of notional particles (normally 300 for each cell) that move from cell to cell, mix and react. Unlike the velocity-composition PDF code, the composition PDF code preserves almost constant the number of notional particles in each cell. The convective fluxes from cell to cell are discretized in order to simulate the mean motion with a finite integer number of notional particles moving from each cell to its neighbors. At each time step, the truncation error between the continuous and the discrete fluxes is updated and new corrected discrete fluxes are computed in order to keep the instantaneous truncation error below the value (20 + 1)/N where D is the geometric dimension of the domain (2 or 3) and N is the average number of particles per cell. Both cylindrical and rectangular grids can be used by the code, without the introduction of weighted notional

5237

particles. After steady state is reached, the code takes samples of all the notional particles at regular time intervals, generating an averaged composition field, thereby reducing the statistical noise due to the use of a finite number of notional particles to represent the PDF. 4.3. Integrated profiles The absorbance fields obtained by the PDF codes have been post-processed in order to compare them with the experimental data. To do this, it has been assumed that the light beam coming from the spectrophotometer crosses the reactor centered on the axis, and that the light intensity in the whole cross-section of the beam (0.28 mm2) is constant. A numerical routine computes the integration of the indicator absorbance, computed from the local mixture fraction PDF using the method described by Li and Toor (1986), in the 3D domain generated by the intersection of the light beam and the tubular reactor for each axial position. A 3D Gaussian integration algorithm with 24 nodes in each of the three main directions (axial and radial direction of the light beam, axial direction of the tube) gives the total absorbance encountered by the light. A suitable interpolation routine has been realized in order to determine the local absorbance in a point of the 3D domain from the 2D absorbance field output by the PDF simulations. All of the integrated values have been normalized with respect to the total absorbance computed by positioning the light beam at the exit of the inner tube. 5. RESULTS

AND

DISCUSSION

5.1. Zsokinetic jet In this section, experimental data for an isokinetic jet with a jet Reynolds number of 10,000 (mean jet velocity of 1 m/s) are compared with PDF simulation results for three different inlet concentration ratios (acid/base): 0.5, 1.0 and 2.0. Simulations have been performed to compare the results from the velocitycomposition and the composition PDF codes. The predictions were found to be identical within the margins of the statistical error generated by the Monte-Carlo algorithms. 5.1.1. Comparison of PDF simulation results and experimental data. In Fig. 4 results found using Corrsin’s mixing time, Q~, are presented. These results are also representative of PDF simulations using Q. i or T+,Z (C, = 2). Indeed, all standard mixing time models underpredict the length of the base plume and the curves fall off too rapidly as a function of axial position. This effect is due to the large values of E predicted by FLUENT behind the inner-tube walls at the injection point. This zone of high turbulent energy dissipation results in the peak molecular mixing rate being located just at the inlet and a concomitant peak in the local mean reaction rate that appears not to be present in the experimental data.

5238

MASSIMO PIPINO

and

Fig. 4. Comparison of simulation results with experimental data for an isokinetic jet with jet Reynolds number of 10,000

(mean jet velocity of 1 m/s). Simulation results found using the velocity-composition PDF code and Q~. Solid line/ squares: concentration ratio = 0.5. Dash-dot line/triangles: concentration ratio = 1.0. Dashed line/diamonds: concentration ratio = 2.0.

Fig. 5. Effect of C, in T+.~ for the isokinetic jet in Fig. 4 with a concentration ratio of 1.0. Simulation results (squares) found using the velocity-composition PDF code. Solid line: C, = 0.5. Dash-dot line: C, = 1.0. Dashed line: C, = 2.0.

We have verified that the PDF codes predict the correct axial profiles for a nonreacting tracer so that turbulent diffusivity and convection are not contributing factors in the overprediction of the molecular mixing rate. 5.1.2. Effect of C, in z, 2_ In an attempt to improve the model predictions, we have varied C+ in t+,,. Results from three runs for a concentration ratio of 1.0 are presented in Fig. 5. All other conditions are the same as in Fig. 4. It can be seen from Fig. 5 that lowering C, does improve the model predictions. However, even with C, = 1.0 where initially the fit is good, the curve eventually falls below the data. Similar results have also been found using T+,~ and r+ 3 modified by an adjustable constant. In all cases, ‘it was found that the mixing rate is initially overpredieted in the injection zone. 5.1.3. Spectral relaxation model. Figures 4 and 5 illustrate the shortcomings of mixing models that neglect the time needed for the large integral-length scale structures formed at the inlet to reduce in size and for the scalar spectrum to become isotropic. The relaxation model introduced in Section 3 accounts for this behavior by letting r; relax to its final value by a cascade process from large to small scales. PDF simulations have been run using the full model and a one-stage stochastic spectral relaxation model with

RODNEY 0. Fox

Fig. 6. Effect of full spectral relaxation model for the isokinetic jet in Fig. 4. Simulation results found using the

composition PDF code. Solid line/squares: concentration ratio = 0.5. Dash-dot line/triangles: concentration ratio = 1.0. Dashed line/diamonds: concentration ratio = 2.0.

Fig. 7. Effect of one-stage spectral relaxation model with C, = 0.18 for the isokinetic jet in Fig. 4. Simulation results found using the composition PDF code. Solid line/squares:

concentration ratio = 0.5. Dash-dot line/triangles: concentration

ratio = 1.0. Dashed line/diamonds: ratio = 2.0.

concentration

C, = 0.18; the results are shown in Fig. 6 and Fig. 7, respectively. It can be seen that the agreement with the data is excellent. We can conclude that the incorporation of spectral relaxation effects in the specification of r+ is essential in order to correctly account for the experimentally observed behavior. 5.2. Nonisokinetic jet In addition to the isokinetic experiments, runs have also been made with an inner-to-outer jet velocity ratio of 3.0 (mean velocity in the coflow of 1 m/s). Comparisons between experimental data and PDF simulations employing the spectral relaxation model are presented in Figs 8 and 9. In Fig. 8 it can be seen that the model slightly underpredicts the mixing rate for concentration ratios of 1.0 and 2.0 in the near field. These predictions can be improved by employing the low Reynolds number correction for turbulent diffusivity (Fig. 2). Simulations run using the spectral relaxation model and the low Reynolds number correction are shown in Fig. 9. 6. CONCLUSIONS

Comparison of PDF simulations for reactive mixing in a tubular jet reactor with experimental data has demonstrated that existing molecular mixing models for the scalar dissipation rate overestimate the mixing rate in the injection zone, resulting in an overestimation of the mean reaction rate and thus an

Reactive mixing in a tubular jet reactor

5239

the model. Nonetheless, PDF simulations appear to be a powerful tool for simulating mixing-sensitive reactions in turbulent flows.

I-

O.OQ

Fig. 8. Comparison of simulation results with experimental data for a nonisokinetic jet with an inner-to-outer velocity ratio of 3.0 and a coflow average velocity of 1 m/s. Simulation results found using the composition PDF code with one-stage spectral relaxation model and C, = 0.18. Solid

line/squares: concentration ratio = 0.5. Dash-dot line/ triangles: concentration ratio = 1.0. Dashed line/diamonds: concentration ratio = 2.0.

Acknowledgements-The authors would like to thank Prof. S. B. Pope for allowing access to his velocity-composition PDF code PDF2DS. In addition, they would like to thank Mauro Fenoglio and Pa010 Baile for collecting the experimental data. This work was supported by the National Science Foundation under grant CTS-9158124 (PYI Award), Dow Chemical Company, and the Italian Government (40% MURST, Multiphase Hydrodynamics). R.O.F. would also like to thank Profs J. Baldyga and J. R. Bourne for the stimulating discussions on scalar mixing during his stay at the Technisch-Chemisches Laboratorium, ETH Zentrum, Ziirich, Switzerland.

NOTATION

constant in the simple Langevin model constants in the spectral relaxation model empirical constant in the k--E model constant in the molecular mixing time constants in the spectral relaxation model strain rate tensor fluctuating part of the strain rate tensor expected value of the strain rate tensor

C=a/axil

Fig. 9. Comparison of simulation results with experimental data for the nonisokinetic jet in Fig. 7. Simulation results found using the composition PDF code with the one-stage spectral relaxation model (C, = 0.18) and low Reynolds number correction for r, (Fig. 2). Solid line/squares: concentration ratio = 0.5. Dash-dot line/triangles: concentration ratio = 1.O. Dashed line/diamonds: concentration ratio = 2.0.

underestimation of the absorbance across the jet plume. This effect has not been accounted for in previous jet reactor models (e.g. Baldyga et al., 1994) where the radially averaged turbulent energy dissipation (E) along the tube axis is modeled using a function that peaks several tube diameters downstream from the inlet. In this study, the distribution of turbulent energy dissipation, computed using the k--E model in FLUENT, shows a strong off-axis peak just behind the inlet-tube walls that results in a locally intense molecular mixing zone. The radially averaged turbulent energy dissipation (or molecular mixing rate) thus has a peak near the inlet before dropping monotonically to its final value. In order to account for the mismatch between the simulation results and the experimental profiles, a novel spectral relaxation model for the scalar dissipation rate has been introduced, which describes the cascade of scalar dissipation from large to small length scales. greatly

PDF simulations have shown that its improves the numerical predictions

use

for the absorbance field. More detailed experimental measurements of the local scalar field and experimental selectivity data will be required to fully validate

R Re,

Re1 S,

se

S S”,” t u U:: <“i>

u u

xi Y+

probability density function gravity scalar flux turbulent kinetic energy reaction rate constant initial turbulence integral length scale initial scalar integral length scale pressure fluctuations mean static pressure potential scalar dissipation rate at the ith stage of cascade scalar dissipation rate scalar dissipation rate following a notional particle ratio time-scale mechanical-to-scalar C= 2klEl turbulence Reynolds number [ = 2.58k/(v&)“2] Reynolds number in the spectral relaxation model [ = k/(vE)‘/‘] chemical source term of species a source term of scalar variance source term of scalar dissipation Schmidt number (= v/T) time velocity fluctuating part of ui expected value of ui mean velocity phase variable of u spatial variable dimensionless distance in the boundary layer

MASSIMO PIPINO

5240 Greek

and

letters

initial integral-length scale wavenumber ratio constant in the modified LMSE model constant in the modifed LMSE model molecular diffusivity turbulent diffusivity turbulent energy dissipation potential scalar dissipation in the ith stage of cascade scalar dissipation potential scalar dissipation in the ith stage of cascade following a notional particle scalar dissipation following a notional particle Batchelor length scale [ = (v/.~)““r’/*] characteristic length scale for ith stage of cascade viscosity kinematic viscosity density turbulent Schmidt number molecular mixing time [ = l/laxil

square magnitude of the scalar gradient variance of the scalar gradient phase variable of C#I

REFERENCES

Baldyga, J., 1989, Turbulent mixer model with application to homogeneous, instantaneous chemical reactions. Chem. Engng Sci. 44, 1175-1182. Baldyga, J. and Bourne, J. R., 1984, A fluid mechanical approach to turbulent mixing and chemical reaction. Chem. Engng Commun. 28, 243-258. Baldyga, J. and Bourne, J. R., 1989, Simplification of micromixing calculations. I. Derivation and application of new model. Chem. Engng J. 42, 83-92. Baldyga, J., Bourne, J. R. and Zimmerman, B., 1994, Investigation of mixing in jet reactors using fast, competitive-consecutive reactions. Chem. Engng Sci. 49, 1937-1964. Batchelor, G. K., 1952, The effect of homogeneous turbulence on material lines and surfaces. Proc. Roy. Sot. A213, 349-366. Brodkey, R. S., 1966, Turbulent motion and mixing in a pipe. A.I.Ch.E. J. 12, 403404. Curl, R. L., 1963, Dispersed phase mixing: I. Theory and effects in simple reactors. A.1.Ch.E. J. 9, 175-181. Corrsin, S., 1964, The isotropic turbulent mixer: Part II. Arbitrary Schmidt number. A.I.CI1.E. J. 10, 870-877. Dopazo, C., 1975, Probability density function approach for a turbulent axisymmetric heated jet. Centerline evolution, Phys. Fluids 18, 397-404. Eswaran, V. and Pope, S. B., 1988. Direct numerical simulations of the turbulent mixing of a passive scalar. Phys. Fluids 31, 506-520. FOX, R. O., 1992, The Fokker-Planck closure for turbulent molecular mixing: Passive scalars. Phys. Fluids A 4, 1230-1244.

RODNEY 0. Fox

Fox, R. O., 1994, Improved Fokker-Planck model for the joint scalar, scalar gradient PDF. Phys. Fluids 6.334-348. Fox, R. 0. and Tsai, K., 1994, The generalized IEM model for reactive scalar mixing. (submitted). Fox, R. O., Hill, J. C., Gao, F., Moser, R. D. and Rogers, M. M., 1992, Stochastic modeling of turbulent reacting flows, in Proceedings of the I 992 Summer Program, Center for Turbulence Research, Stanford. Gao, F., 1991, An analytical solution for the scalar probability density function in homogeneous turbulence. Phys. Fluids A 3, 511-513. Hinze, J. O., 1975, Turbulence, 2nd edn, McGraw-Hill, New York. Janicka, J., Kolbe, W. and Kollmann, W., 1979, Closure of the transport equation for the probability density function of scalar fields. J. Non-equil. Thermodyn. 4, 47-66. Jiang, T. L. and O’Brien, E. E., 1991, Simulation of scalar mixing by stationary isotropic turbulence. Phys. Fluids A3, 1612-1624.

KosBly, G., 1989, Scalar mixing in isotropic turbulence. Phys. Fluids A 1, 758-760. Leonard, A. D., 1989, Direct numeric1 simulations of chemically reacting turbulent flows. Ph.D. thesis, Iowa State University. Li, K. T. and Toor, H. L., 1986, Chemical indicators as mixing probes. A possible way to measure micromixing simply. Ind. Engng Chem. Fundum. 25, 719-723. Mantel, T. and Borghi, R., 1991, A new model of premixed wrinkled flame propagation based on a scalar dissipation equation. Preprint. Newman, G. R., Launder, B. E. and Lumley, J. L., 1981, Modelling the behaviour of homogeneous scalar fields. J. Fluid Mech. 111,217-232. Norris, A. T., 1993, The application of PDF methods to piloted diffusion flames. Ph.D. thesis, Cornell University. O’Brien, E. E., 1980, in Turbulent Reactive Flows (Edited by P. A. Libby and F. A. Williams), Ch. 5. Springer-Verlag, Berlin. Pope, S. B., 1982, An improved turbulent mixing model Comb. Sci. Technol. 28, 131-145. Pope, S. B., 1985, Pdf methods for turbulent reactive flows. Prog. Energy Cornbust. Sci. 11, 119-192. Pope, S. B., 1991, PDFZDS: A code to solve the velocitycomposition joint pdf equation, unpublished FORTRAN code listing. Tsai, K. and Fox, R. O., 1994, Modeling the effect of turbulent mixing on a series-parallel reaction in a tubular reactor. ICRES Report 9403, Kansas State University. Valifio, L. and C. Dopazo, 1991, A binomial Langevin model for turbulent mixing. Phys. Ffuids A 3, 3034-3037. APPENDIX

The Reynolds-averaged equations for a scalar, its variance, and the magnitude of its gradient can be derived starting from the governing equation for a scalar field in an incompressible turbulent flow:

Dd 84 +,.?_=~a29 _=_ Dt

dt

axiaxi

'axi

(Al)

wherein repeated indices imply summation. Defining xj = differentiation of eq. (Al) with respect to xjr @lax,. yields DXj

ax.

ot---‘tuizI:r at

axj

@X,

axiaxi

ejiXi

C.42)

where we have introduced the strain rate tensor e,, = &/ax,. The equation for the magnitude of the scalar gradient x2 = xjxj follows from eq. (A2): 1 Dx’ 1 ax= i ax2 --=2at+*2Uiax=I-,yj

2 Dt

Reynolds averaging

of eq. (Al)

a% - xjejixi. ax,ari

-

(A3)

yields the well-known

Reactive mixing in a tubular jet reactor expression for the mean scalar field:

04 iG

x4) =dt+h>

( >

aa+ a 7-.+a-d4’>=y3-&. x*

xc

aZc4)

The equation for can be derived from eq. (A2) by first multiplying both sides by ,r; and then averaging:

XL , (A4)

The term involving is the scalar flux term that must be closed, usually by invoking the gradient-diffusion assumption. The final term can be neglected in high Reynolds number flows. Multiplying eq. (Al) by 4’ and averaging yields the equation for the scalar variance:

5241

-2

1 Kx”> _at ah,)

- - .

(A7)

The terms involving second-order spatial derivatives of mean quantities can be neglected in high Reynolds number flows. Rearranging the terms and using = a<4)lax,, we find The next-to-last term is this expression can be neglected in high Reynolds number flows. Rearranging the terms, one finds the final expression for the scalar variance: