Turbulent Chemical Reactions : Application to Atmospheric Chemistry

Citation

Shu, Winston Rei-Yun (1976) Turbulent Chemical Reactions : Application to Atmospheric Chemistry. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/c507-9y29. https://resolver.caltech.edu/CaltechTHESIS:01202026-223733813

Abstract

The general problem of predicting the rates of chemical reactions in turbulent fluid has been studied. Particular attention has been given to bimolecular (second-order) chemical reactions occurring isothermally. The problem of predicting the rates of turbulent chemical reactions has been formulated in both Eulerian and Lagrangian frameworks. The Eulerian formulation leads to a closure problem in the equations for the moments of concentrations of the reacting species. The Lagrangian formulation, while circumventing the usual closure problem, is based on a transition probability density function of particle displacements which can only be predicted in the simplest of cases.

The principal result of the work is the development of a new closure method for second-order chemical reactions in turbulence. The method is referred to as the diffusion zone model. The essential idea of the model is that one defines the diffusion zone as that region of the fluid in which the reactive species would coexist if they were inert. Then, chemical reaction is presumed to occur only within the diffusion zone; outside the diffusion zone the concentrations of the reactive species are identical to those if the species were inert. The diffusion zone model is tested extensively on available laboratory data and gives good agreement.

The diffusion zone closure method is adapted for predicting the rates of atmospheric chemical reactions occurring in plumes. The effect of inhomogeneous mixing on the reaction rates is explicitly accounted for in the formulation. In assessing this effect, a numerical planetary boundary layer velocity field is employed as a source of the necessary turbulence data. The model is applied to predicting the rates of nitric oxide oxidation in a plume from the Morgantown power plant in Maryland. Predictions are in good qualitative agreement with available measurements on the plume.

Item Type:

Thesis (Dissertation (Ph.D.))

Subject Keywords:

(Chemical Engineering and Environmental Engineering Science)

Degree Grantor:

California Institute of Technology

Division:

Chemistry and Chemical Engineering

Major Option:

Chemical Engineering

Minor Option:

Environmental Science and Engineering

Thesis Availability:

Public (worldwide access)

Research Advisor(s):

Seinfeld, John H.

Thesis Committee:

Unknown, Unknown

Defense Date:

8 December 1975

Funders:

Funding Agency	Grant Number
NSF	UNSPECIFIED
California Institute of Technology	UNSPECIFIED
Systems Applications, Inc.	UNSPECIFIED

Record Number:

CaltechTHESIS:01202026-223733813

Persistent URL:

https://resolver.caltech.edu/CaltechTHESIS:01202026-223733813

DOI:

10.7907/c507-9y29

Default Usage Policy:

No commercial reproduction, distribution, display or performance rights in this work are provided.

ID Code:

17829

Collection:

CaltechTHESIS

Deposited By:

Benjamin Perez

Deposited On:

26 Jan 2026 17:58

Last Modified:

26 Jan 2026 19:15

Thesis Files

PDF - Final Version
See Usage Policy.
46MB

Repository Staff Only: item control page

TURBULENT CHEMICAL REACTIONS : APPLICAT ION TO ATMOSPHERIC CHEMISTRY Thesis by WiMton Ru-Yun Shu In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena , California 1976 (Submitted December 8, 1975) ThMi the/2..L6 ,v., de,cuc.a:te,d to my mothe.Jt -iiiACKNOWLEDGMENT The author wishes to express his sincere appreciation to Professor J. H. Seinfeld for his guidance, interest and encouragement throughout the course of this study. Special thanks go to Dr. R. G. Lamb, without whose contribution of enlightening ideas and patient discussions this work would not have been possible. Financial support from the California Institute of Technology in the form of Graduate Research and Teaching Assistantships and from the National Science Foundation is gratefully acknowledged. Part of of the numerical work was supported by Systems Applications, Inc. The author is also indebted to his wife, Gloria, for her collaboration in preparing this thesis, and to Mrs. Lenore Kerner for her excellent typing. -ivABSTRACT The general problem of predicting the rates of chemical reactions in turbulent fluid has been studied . Particular attention has been given to bimolecular (second-order) chemical reactions occurring isothermally. The problem of predicting the rates of turbulent chemical reactions has been formulated in both Eulerian and Lagrangian frameworks. The Eulerian formulation leads to a closure problem in the equations for the moments of concentrations of the reacting species. The Lagrangian formulation, while circumventing the usual closure problem, is based on a transition probability density function of particle displacement s which can only be predicted in the simplest of cases. The principal result of the work is the development of a new closure method for second-order chemical reactions in turbulence. The method is referred to as the cu.66Mion zone mode.£. The essential idea of the model is that one defines the diffusion zone as that region of the fluid in which the reactive species would coexist if they were inert. Then, chemical reaction is presumed to occur only within the diffusion zone; outside the diffusion zone the concentrations of the reactive species are identical to those if the species were inert. The diffusion zone model is tested extensively on available laboratory data and gives good agreement. The diffusion zone closure method is adapted for predicting the rates of atmospheric chemical reactions occurring in plumes. The effect of inhomogeneous mixing on the reaction rates is explicitly accounted -vfor in the formulation. In assessing this effect, a numerical planetary boundary layer velocity field is employed as a source of the necessary turbulence data. The model is applied to predicting the rates of nitric oxide oxidation in a plume from the Morgantown power plant in Maryland. Predictions are in good qualitative agreement with available measurements on the plume. -viTABLE OF CONTENTS ACKNOLEDGMENT iii ABSTRACT iv I. II. INTRO DU CTI ON 1 1. The Importance of Concentration Fluctuations 3 2. Previous Work 11 3. Scope of this Thesis 18 A DIFFUSION ZONE MODEL FOR TURBULENT CHEMICAL REACTIONS 20 1. Development of the Diffusion Zone Model 2. Estimation of Parameters 20 28 2.1 The Mixing Parameters 30 The Reaction Parameter 37 2.2 3. Application of the Diffusion Zone Model to Laboratory Data 49 4. Generalization of the Diffusion Zone Model to Three Spatial Dimensions 72 III. APPLICATION OF THE DIFFUSION ZONE MODEL TO CHEMICALLY REACTING PLUMES 75 1. Diffusion Zone Model for a Chemically Reacting Plume 78 2. Simulation of a Plume in the Atmospheric Boundary Layer 84 Description of the Velocity Field 84 2.2 Sub-Grid Scale Velocity Simulation 87 2.1 2.3 Diffusion Experiment in Simulated Atmospheric Turbulence 3. Concentration Statistics in a Simulated Plume 96 103 -vii- IV. V. 3.1 Simplification of Equation (3.35) 111 3.2 Sampling Procedure and Statistical Error of the Estimate of Af 115 3.3 121 Numerical Results 4. Study of the N0-03 Reaction in Simulated Atmospheric Turbulence 124 5. Discussion 137 LAGRANGIAN MODELS OF TURBULENT CHEMICAL REACTIONS 140 1. A Laqrangian Description of a Second-Order Chemical Reaction in a Turbu1ent Fluid 142 2. Lagrangian Models of Turbulent Chemical Reactions 152 3. Discussion 163 CONCLUSION AND RECOMMENDATIONS FOR FURTHER WORK 165 APPENDICES A. Multijet Reactor Data 170 B. Numerical Calculation of Reaction Parameter 179 C. Estimation of NO-Emission for a Power Plant 183 D. Considerations of Multiple Kinetics 184 REFERENCES 189 -1I. INTRODUCTION Chemical reactions in turbulent fluids occur in situations of both industrial and natural significance. The main problem in the de- sign of industrial and in the analysis of natural systems is the prediction of the rate of the chemical reactions. Apart from the labora- tory ideals of the perfectly mixed vessel and pure plug flow, chemical reactions in turbulent flow will be influenced by the degree of turbulence and mixing in the system. To illustrate the variety of reactive flows in which turbulence is important; Table 1.1 lists some recent studies. One of the major problems in the analysis of the dynamic behavior of air pollutants is predicting the rates of chemical reactions occurring among atmospheric species. The atmosphere is, of course, in a turbulent state, and since reactive pollutants are always emitted from highly localized sources, there is a period during which the pollutants are inhomogeneously mixed and reactions are taking place. During this period, depending on the speed of the reactions, the speed and extent of mixing may have a significant effect on the reaction rates. tive of this work is twofold. The objec- First, we wish to study the general prob- lem of predicting the rates of chemical reactions occurring in turbulent fluids. Second, we wish to study particularly the effect of atmospheric turbulence on the rates of chemical reactions involving air pollutants. e shall concentrate here on isoH ermd l second- order reactions in turbu- \•1 lent flovJ. Continuous flow sti r red-tank reactor Nitric oxide formation in combustors Turbulent diffusion flame structure Chemical kinetics measurements Turbulent jets and wakes Photochemical smog formation Pollutant disper sion from smokestac ks 2. 3. 4. 5. 6. 7. 8. *Hill (1976) has a more extensive but somewhat different list . Tubular plug -flow reactor l. Davis et al . (1974), Davis and Klauber (1975) Romber g (1974) Donaldson and Hi lst (1972ab) Borghi (1974) , Shea (1976) Glassman and Eberstein (1963) Chung (1971 , 1972) Flagan and Appleton (1974) Evangelista et al . (1969) Toor and Singh (1973) References Examples of Chemical Reactions in Turbulent Flow* Description of Physical Phenomena Table 1.1 I I N -3- I.l The Importance of Concentration Fluctuations To introduce the problem with which we will be dealing, let us consider the isothermal, bimolecular chemical reaction k A+ B - Products occurring between species A and B with rate constant k. The local r ate of change of the concentrations of A and Bin a turbulent system are aA _ aB = - kAB at - ar ( l. l) where k is the kinetic rate constant, and A and Bare the instantaneous values of the concentrations of species A and B, respectively . The in- stantaneous concentrations can be expressed as the sum of ensemble average* and fluctuating concentrations: A= <A> + A1 ( l . 2a) B = + 8 1 ( l . 2b) Substituting (1.2) into (l.l), and ensemble averaging, we find that the mean rate of reaction is given by a<A> = a = _ k( <A> + <A1B1>) at at ( l . 3) *Ensemble average, or statistical average, of a certain random variable is the average value of the ensemble of an infinite number of statistically independent measurements of the variable. If not otherwise ~entioned we shall use 1mean 1 or average 1 to refer to ensemble average 1n the following texts. 1 -4One is generally interested in predicting the ensemble average rate of reaction. The approximation often used in computing this rate is to replace the instantaneous, local expressions (1.1) which are derived from the chemical kinetics alone, by the corresponding expressions involving the ensemble mean concentrations, i.e. a<A> at a at k<A> -- = -- = - (l.4) We see, however, from (1.3) that only when the correlation of the fluctuations <A'B'> is identically zero is (1.4) an exact relation. If <A'B'> > 0, the mean reaction rate is enhanced over that predicted by (1.4); and if <A'B'> < 0, the rate is suppressed. Thus, in order to use (1.4), it is necessary to demonstrate that <A> >> l<A'B'>I the particular system of interest. for Donaldson and Hilst (1972b) have suggested a hypothetical example in which the reaction rate predicted by (1.4) is totally incorrect. -, Consider a situation in which the flow of material by a given point is such that alternate patches of A and B pass the point. Let us suppose that half the time the flow is all A and half the time it is completely Bas depicted below. keeps repeating Concn. ' 1. A B r---, I I I I I I A B r---1 I I I I I I I I I I I I I I I I I I time If this pattern -5and the patches of A and Ball have a concentration of unity, the average concentrations are obviously <A>= = 1/2. When the flow is all A, A' = 1/2 and B' = - 1/2, but when the flow consists of B only, A' = - 1/2 and B' = 1/2. Thus for this flow <A> = 1/4 <A'B > = - 1/4 1 and by (1.3) the correct reaction rate is zero; but according to (1.4) the rate is k/4. The second order correlation <A 1 B1 > arises as a result of inhomogeneous mixing in the system. The next question that naturally follows is - Is it possible to predict <A 1 B1 > for a system of interest and thereby utilize (l .3) instead of (1.4) for predictions? Diffu- sion of one substance through another is usually approached by considering a differential material balance in a fixed volume of space. The resulting differential equation is called the equation of mass continuity and can be solved as a boundary value problem describing the behavior of the concentration of the material of interest. Let us con- sider the equations of continuity for A and Bin a fluid, ~+ at aA _ V a2 A uk axk - A axkaxk - kAB 2 ~+ u ~= VB a B at k axk axkaxk - kAB {l. 5) ( 1 . 6) -6- where uk is the fluid velocity in coordinate direction k (k al ,2,3), and VA and v are the molecular diffusivities of A and Bin the fluid. 8 The summation convention has been employed in (1.5) and (1.6) in which a repeated subscript in a term indicates summation over the three components of that term, e.g. u ~= k axk ! k=l u ~ k axk (l. 7) In a turbulent fluid, the velocities are customarily decomposed into mean and fluctuating components, (1.8) Introducing (l .2), and (1.8) into (1.5) and (1.6) and ensemble-averaging the resulting equations leads to a<A> + a<A> + - a <u'A'> = V ~t k -axk ~ A O oXk k ( l. 9) (1.10) where in obtaining the third term of the left-hand side we have employed the total continuity equation for incompressible flow au'k axk =0 (1.11) -7we see that (1.9) and (1.10) contain the seven new dependent variables, <ukA 1 >,<ukB 1 >(k = 1 ,2,3), and <A 1 B1 >. The last terms on the right hand sides of (1.9) and (1.10) are, of course, the local rates of reaction discussed earlier. Perhaps the first inclination in dealing with the seven new dependent variables is to derive the equations which govern their dynamics. Doing so, we obtain .1... + <u.> _a_ <uk1 A1 > + <uk1 u !> a<A> + a at k , ax. ax. 1 ~ 1 <u'u'A 1 > ki 1 1 a2 <uk.A'> -- VA - - - k[<uk'A'> + <uk'B 1 ><A> + <uk'A'B 1 >] ax.-ax. 1 (1.12) 1 a + <U,> - a + <uku.> , • a + _a_ -at ax. k , k 1 ax. k , -ax. 1 1 1 1 1 a2 <u'B'> = v _ _k__ - k[<A> + + ] k k k 8 ax.ax. , 1 a a<A> 1 - - <u'A > k axk k axk ~ <A'B 1 > + <Uk> _a_ <A'B'> axk "I, (2) ( 1) - (1.13) _a_ + V a <A 8 > 2 axk 1 k axkaxk (3) (4) 1 - k [<A'B'> (<A>+) + <A'B' 2 > + <A' 2 B'> + <B' 2 ><A> + <A' 2 >] (6) (1.14) -8- where in (l .14) we have assumed DA= v8 = V for simplicity. The num- bered terms in (l .14) have the following significance: ( l ) = convection of <A B > by mean velocity 1 1 ( 2) = generation of <A 8 > by mean gradients in <A> and 1 1 (3) = transport of <A1 8 1 > by turbulent velocity fluctuations (4) = molecular diffusion of <A1 8 1 > ( 5) = dissipation of <A'B' > by molecular diffusion (6) = production and decay of <A1 B1 > due to chemical reaction Although we can ascribe a significance to the terms in (l .12-14) these equations cannot be solved in closed form because of the appearance of still more unknowns. the closure problem. This phenomenon in turbulence is known as Some work has been done on propo s ing methods* to close the hierarchy of moment equations for turbulent chemistry, and we shall review this work shortly. First, however, qualitative condi- tions can be developed for when terms of the form <A'B' > can be expected to be significant (Donaldson and Hilst, l972ab). In the absence of mean gradients, two competing processes influence the behavior of <A'B'>: *Closure methods are hypothesized based either on physical insight or strictly analytical approximations. The most common closure methods that have been proposed for the turbulent fluxes <u~c' >, for example, is Prandtl 's mixing length theory (Hinze, 1959) in which <UC > = - K 3<C> k 3X k I I where K is a so-called eddy diffusivity. -9- molecular diffusion and chemical reaction. When two species react to form another species at a certain point, the concentrations of these two reactants instantly become smaller and as a result concentration fluctuations are produced. These fluctuations are dissipated by mo- lecular diffusion (as expressed by term (5) in Eq. (l .14)); but if the rate of dissipation is significantly slower than that of the chemical reaction, then the reactants tend to become segregated and the overall rate of the reaction is perturbed. The ratio of the magnitudes of the reaction and dissipation processes is expressed by the quotients of terms (5) and (6) in Eq. (1.14). Expressing the dissipation in terms of a dissipation length scale ;\(Corrsin, 1958), i.e. 2V<A B > 1 2V < 1 (1.15) ---:;--2--- /\ and assuming that the reaction rate term (6) in (l .14) is of the same order of magnitude as k( <A>+) <A B > ' 1 1 we can write 2 k( <A 2V>+);\ ___ _ (1.16) where TD and TR are the time scales of the dissipation and reaction processes, respectively. Damkohler number. The dimensionless quantity N is called the ~/hen N< < l, the time scale of the reaction is much larger than that of the dissipation ,and concentration fluctuations are eradicated before they had time to affect the chemistry. In this -10- case the mean reaction rate can be adequately predicted by -k<A>. on the other hand, when N>> l, the characteristic time for chemical reaction is short compared ·to that for molecular dissipation. In that case, <A 1 8 1 > will tend to -<A>, and the two species will be poorly mixed. The rate of reaction between A and B will be gover ned not by the kinetics but by the rate of molecular diffusion. A characteristic value of A in the atmosphere is 10 cm. Using a typical value of 0.17 cm 2/sec for V and 0.2 ppm for (<A >+), we calculated values of N for several reactions occurring in photochemical smog and listed them in Table 1.2. Although it is not strictly correct to consider reactions individually which occur in a complicated sequence, we do obtain, nevertheless, an indication that many atmospheric reactions might be 11 diffusion - limited. 11 The extent of the importance of this effect will be examined later in this work. We shall first review previous work relating to the study of this effect. Table 1.2 Estimated Ratios of Characteristic Diffusion Time to Characteristic Reaction Time for Important Atmospheric Reactions Reaction N0 2 + o = N0 + o 3 3 2 N0 3 + NO = 2N0 2 NO+ Ho 2• = NO + ow 2 *Seinfeld (1975) * TD N= TR 2.9xl0 1 4.6xl0- 2 28.4 l. 5xl04 7.0xlO 2 14705.8 k -1 (ppm-min) 0.045 686.3 -11- 1.2 Previous Work Perhaps the first work on the statistical theory of turbulent chemical reactions was that of Corrsin (1958), who considered decay of a single species undergoing an isothermal irreversible reaction in isochoric, isotropic turbulence. He concluded that first-order reactions possess the same general character as those of pure mixing without reaction, but for second-order reactions, the nonlinearity of the kinetics introduces additional unknowns and makes the problem decidedly different from the pure mixing one. Considerable interest has since been focused on the spectral analysis of reactive mixtures in turbulent fluids (Corrsin, 1961, 1962, 1964; O'Brien and Francis, 1962; Pao, 1964). These calcu- lated spectra, however, have not yet received experimental support. (Hill, l976i Later reported work concentrates more on the prediction of statistical moments of concentrations of reactive species as shall be reviewed below. All approaches to turbulent phenomena may be grouped into two broad categories: Eulerian and Lagrangian. The former begins with the mass continuity equation (l .5-6) and at some point introduces a closure hypothesis to achieve a solvable set of equations. The Lagrangian ap- proach arrives at a closed set of model equations by postulating functional forms for certain probability density functions of the displacements of fluid particles*. In the remainder of this section we review *Manin and Yaglom (1971) define a fluid particle as a lump of fluid that is large compared to the mean free path of molecules of the fluid but small enough so that within the framework of continuum mechanics it may be identified as a point moving in the fluid. -12previous attempts of both the Eulerian and the Lagrangian types to analyze the problem of isothermal reactions* in turbulent flows. Eulerian Approaches (Closure Methods) - There have been several studies of the simple Eulerian equation dA(t) = _ kA 2 (t) dt (1.17) with stochastic initial conditions A(O). O'Brien (1966) compared ex- pressions for <A(t)> and <A12 (t)> obtained by averaging the exact solution of Eq. (1.17) with the predictions of these quantities using the third moment discard, quasinormal, and direct interaction approximations. None of the approximations behaved satisfactorily when the relative amplitude of the initial fluctuations, viz <A 12 (0) >/ <A(0) >2 , was large. Later O'Brien (1968) recognized that since A is a nonnegative random variable, the moments of A must sati sfy Liapunov s inequality (Uspen1 sky, 1937): <Ab>a-c ~ <Ac >a-b <Aa >b-c , a > b > c > 0 where a, band care constants. (1.18) As an illustration, the third moment of 1 A must satisfy (1.19) *For treatments of nonisothermal reactions in turbulent flows, see Libby (1972), Dopazo and O'Brien (1973), Flagan and Appleton (1974). -13It is obvious that if we adopt, say, the third moment discard approximation, <A 12 >/<A> 2 must be.:::._ 1 at all times. strictions on the initial fluctuations, To avoid unrealistic re- O'Brien (1968) proposed a so- called inequality preserving closure approximation (IPCA) for <A13 > such that Liapunov's inequality is always satisfied. O'Brien again compared the predictions of his IPCA with the exact solutions of (1.17) and found satisfactory agreement. O'Brien and Eng (1970) generalized the closure for reaction of order one to three, and O'Brien and Lin (1972) used a different IPCA for two-species reaction with spatial variables. Finally, in a recent study Lee (1973) presented a generalized direct interaction approximation to (1 .17). All of these approximations appeared to behave satisfactorily when tested against (1.17). Since each of these closure schemes was developed and tested on a simple system, namely (l .17), there is no assurance that they will hold when applied to the full continuity equations (1.5-6). One of the few closure schemes that have been applied to the full continuity equations is that of McCarthy (1970). By discarding fifth order cumulants, he developed a hierarchy of 78 differential equations for single-point concentration moments and microscales. Later Lin and O'Brien (1972) presented a closure theory which incorporates Lin's (1971) third order IPCA for the reaction terms and Lee's (1966) modification of the quasinormal approximation for the convective terms. Computations of the decay of moments and spectra of A were carried out for various conditions. The decay of <A> and <A12 > was found to -14- depend primarily on the Damkohler number defined as in (1. 16). More , recently, Hilst et al. (1973) combined a third order IPCA (different from that of Lin, 1971) for the reaction term A'B model" (Donaldson and Rosenbaum, 1969) vection term. 1 with an "invariant for the con- Hilst et al. then applied the model to the reaction of o and NO emanating from four cross-wind freeway line sources and 3 solved the resulting 12 coupled differential equations numerically. Due to their mathematical complexity and the lack of experi- mental support, none of the closure schemes mentioned above is particularly useful in practical applications to chemical engineering and atmospheric pollution. Closure schemes more suited to these needs have been developed by Toor, Brodkey, and their students. Toor (1962) noticed that by subtracting (1 .6) from (1 .5) the equation governing C = A - Bis identical to th e transport equation of an inert species: ( 1. 20) provided that VA= v = V. Since for very rapid reactions A and B do 8 not coexist, C is equivalent to A when C ~ 0 and it is equivalent to B when C < 0. The probability density functions for A and Bare thus the same as for C in the proper range, and 00 <A> = = J cp(C) dC 0 ( 1. 21) -15Following the work of Hawthorne et al. (1949), Toor assumed p(C) to be Gaussian, and an explicit result was obtained in terms of stoichiometric feed ratio S and the mean and variance of C. For stoichiometric feed (s=l) the result is <A(x,t)> = ( l . 22) at any point x in the fluid . In a one-dimensional tubular reactor, a more useful expression is _ <A(t)> _ <B(t)> <B(O)> J t) = <AToJ> - F ( l . 23) The quantity I 5 can be measured by an identical but nonreacting system . Toor (1969) further showed that for the case of initially unmixed reactants Is is equivalent to the intensity of segregation: <AI BI ( t) > I I <AII BI1 ( 0) > ( l . 24) where A1, BI are concentrations that species A and B would have if they were inert . Toor s theory has since received extensive experimental 1 support (Vassilatos and Toor, 1965; Keeler et al. 1965; Mao, 1969 ; Torrest and Ranz, 1970; Miyairi et al. 1971; Miyawaki et al. 1974) . Although his theory was based on the assumption of VA= v8 , a generalization to nonequal diffusivities has been reported recently (O'Brien, 1971 ). We shall discuss Toor s method in more detail in Appendix A. 1 -16For reaction of intermediate speed, Toor s model does not 1 apply because reactants do coexist . To handle these cases Toor (1969) proposed for reactions of any speed with stoichiometric feed ( s=l) that (1.25) This postulate was based on the observation that (1 . 25) is valid in the limit of both very slow and very rapid stoichiometric reactions . (The latter follows from (l.23-24) and the fact that <A B > -+ - <A> .) 1 1 Yieh (1970) provided additional theoretical support for (l .25) by assuming statistical independence of the product and reactant concentrations (although a negative correlation would be expected). Mao and Toor (1971) and Ajmera et al. (1975) tested (1.25) in both liquid and gas phase reactions in an idealized one-dimensional reactor * and found it to be reasonably accurate. Based on these findings Toor and Singh (1973) later proposed scale-up criteria for tubular reactors. In addi- tion to this model, several others have suggested schemes for describing reactions of intermediate speeds, (e.g., Kattan and Adler, 1967; Mao and Toor, 1970; Rao and Dunn, 1970; Evangelista et al. 1969; Kattan and Adler , 1972). Some of their results are discussed in Section II . 3. *More facts about this reactor and the experiments and models associated with it are given in Section II.3. and Appendix A. -17Lagrangian Approaches Lagrangian treatments of turbulence phenomena are relatively few in number because the particle statistics required by this approach are generally extremely difficult to determine . In a series of papers, Chung (1969, 1970, 1972) used probability density functions (pdf) of a model velocity field which obeys a Langevin equation to study reacting turbulent shear flow. The resulting Fokker-Planck equations were solved by the method of moments, assuming Chung s formulation has been questioned by 1 Maxwellian distribution. Bush and Fendell (1974) and 0 Brien (1974). 1 Different Lagrangian formu - lations were proposed by Hill (1970) and Riley (1973) for one-point concentration moments, but solutions were not given. Recently, Lamb (1972b) derived an expression which relates reactive decay of a species to the probability density function that describes the paths of inert fluid particles. Lamb s approach will be discussed further in Chapter IV. 1 Brief Summary Despite the number of previous studies, there does not exist a generally valid means of predicting the rates of chemical reactions in turbulent fluids. Toor s (1962) theory for very rapid reactions is 1 perhaps the most useful, but it has not been tested in systems other than chemical reactors. Schemes that appear to be more general, such as the IPCA s, are very difficult to apply due to their mathematical 1 complexity. -1 8- 1. 3 scope of this Thesis As noted earlier, there are two approaches available for turbulent phenomena. The Eulerian approach starts from t he familiar species continuity equation, but it is inevitably frustrated by the closure problem. Although a number of approximate methods can be in- voked for closing the hierarchy of equations, an exact solution cannot be obtained. However, an important issue is the development of simple, yet accurate, closure approximations for terms of the form <A'B' > which arise when chemical reactions take place. The development of such ap- proximations would enable the use of the Eulerian continuity equations in calculating turbulent, reacting flows. More importantly, perhaps, such a development would enable one to assess when predictions based on mean-value kinetics alone will not suffice. Chapter II is devoted to the development of a new closure approximation for turbulent, reacting flows. The closure model is then tested against laboratory data. The application of this model to atmospheric chemistry is the subject of Chapter III. A new formulation will be presented to predict the reac- tive decay of the species emitted from smokestacks. To assess the effect of turbulence on atmospher ic reaction rates, we employ the numerical planetary boundary layer flowfield of Deardorff (l974ab) as a source of the necessary turbulence data. The resulting prediction is then com- pared with some atmospheric observations. The Lagrangian approach to fluid mechanical problems describes Phenomena with respect to specific fluid particles. In the case of the -19turbulent diffusion of inert species, the Lagrangian formulation circumvents the closure problem which arises in the Eulerian description (Lamb, 1972a). However, the key va ria ble in the Lagrangian formulation, the transition probability density function (pd f ), is very difficult to predict theoretically or measure experimentally . When species may react chemically, the Lagrangian formulation requires that individual particle trajectories be known in order to account fo r species. 11 collisions of reactive 11 A numerical technique of measuring the transition pdf is demon- strated in Chapter III, while Chapter IV is devoted to a Lagrangian analysis of turbulent chemical reactions . Although the analysis is quite complicated, results of interest are obtained in some special cases . Conclusions and recommendations for future work are given in Chapter V. -20- II. A DIFFUSION ZONE MODEL FOR TURBULENT CHEMICAL REACTION The object of this chapter is to develop a closure scheme fo r second-order turbulent reaction systems . A so-ca l led di ffusion -zone model will be presented based largely on physical insight rather t han on mathematical assumptions. Applications of this scheme to si mu l ate reactor data are carried out, and agreement with data is shown to be ve ry good. A comparison with other related work is also given at the end of the chapter. II.l Development of the Diffusion Zone Model There are three phenomena influencing the concentration sta - tistics of reactive species in a turbulent fluid: turbulent mixing, molecular diffusion, and chemical kinetics . The relative roles of these three phenomena are illustrated in Figure 2.1, in which we depict the turbulent mixing of two fluid elements containing reactive species A and B. Turbulent mixing serves to stretch the patches of A and B into thin sheets and filaments as shown in Figure 2.1 (a) . However, the turbulence i s not capable of mixing the patches at the very smallest scales . The region of small-scale inhomogeneity can be estimated to have a length of h the order of (vV2 / E) 4 , where vis the kinematic viscosity of the fluid , V the molecular diffusivity of the species (for the sake of qualitative argument, assume VA= v8 = V), and Ethe rate of energy dissipation of -21- la (a) p (b) Figure 2.1 Schematics of turbulent mixing (a) without molecular diffusion and (b) with molecular diffusion. -22the turbulence (Tennekes and Lumley, 1972, p. 240). This length scale can be referred to as the contconinant microscale , first introduced by Batchelor. For chemical reaction to occur, it is necessary that the patches of A and B mix. Such mixing , on scales smaller than the con- taminant microscale, can only be achieved by molecular diffusion. Let us call the region in which this diffusion takes place the diffusion zone (Figure 2.1 (b)). If A and Bare inert, the diffusion zone is simply the region vAB in which A and B coexist. When A and 8 react with each other, the region in which they coexist is often called the reaction zone. It is normally smaller than vA 8 due to consumption of the species by reaction. For simplicity we develop the model considering only a single second-order chemical reaction of the form A+ B ➔ products and only a single spatial dimension. Subsequently, the model will be generalized to three space dimensions (Section II.4) Consider a volume V of a turbulent fluid into which species A and Bare released. At any given point in the fluid, the rates of change of the concentrations of A and B due to chemical reaction are given by dA _dB_ ctf - dt - ~ kAB. (2.1) We now define the volume averaged quantities A=} J Adv, V (2.2a) -23- = V l J Bdv, B (2.2b) V AB = ~ J ABdv (2 . 2c) V where A and B denote the concentrations of species A and B, respectively . Taking the volume average of (2.l) we have dA _dB_ (2.3) cit - dt - - kAB . If the volume Vis sufficiently large and if the joint mean values of A and Bare homogeneous (independent of location) and meet certain other weak requirements: then according to ergodic theory the volume averages defined by (2.2) are equivalent to ensemble averages and hence the former are invariant from one realization of the turbulent-chemical process to another . In the analyses that follow, we shall assume that all averages over the volume V possess this invariance property. Suppose that in the places of the reactive species A and B inert materials A1 and B1 are released instead. In this case we define vAB to be the total volume of fluid within Vin which A1 and B1 are found together. In practice, one might consider the edge of a particu- lar portion of vAB to be the point where a probe indicates that the concentration of A1 or B1(or both) is vanishingly small. In any event, VAB is defined only for the sake of convenience and it does not enter explicitly into the final mathematical form of the closure model equations. In a similar manner we define vA and vB to be the volumes within V --- * Such as the spatial correlation <A(x+6x)B(x)> goes to zero as 6x goes to infinity. See Monin and Yaglom (1971) for more details. -24in which the species AI is found alone and BI is found alone, respectively. (See Figure 2.lb.) It is important to keep in mind that vA, v , and vAB are defined in terms of the inert rather than the reactive 8 scalar species. Since AI(B ) exists only in vAB and vA (v 8 ), the volume aver1 ages analogous to (2.2) can be expressed in the following forms: - l - l A A A A (2.4a) AI= V (vAAI + vABal) BI= V (vBBI + vABbl) l A (2.4b) " AIBI = V µIalbivAB (2.4c) where l BI - -VB f B1dv , (2.5) = _l_ VAB f B1dv (2.6) A a1 = - l vAB A f Adv A bl I VAB VB VAB and (2. 7) Under the assumptions invoked earlier that Vis large and that the statistics of the scalar fields satisfy the conditions of the ergodic s theorem, the volume averages A1, 1 , and A1B1 are invariant from one realization of the turbulence to another. Hence, according to -25Eqs. (2.4), in general the quantities defined in (2.5)-(2.7) also possess this invariance 'property. Returning to the reactive species, we postulate that the mean values given by Equation (2 . 2) can be expressed in the following forms: - l " - l " " A= V (vAAI + vABa) l (2 .8a) ' (2.8b) A B = V (vBBI + vABb) - ' AA (2.8c) AB = V µabvAB where f Adv 1 a =- - VAB (2.9) (2.10) " " and where vA' vB' vAB' AI' and BI have the same meanings as before. The integrals in Eqs. (2.9) and (2 .10) denote integrations over the region vAB where, all things being constant, A and B would be found together if they were inert. Due to the depletion of both species A and B by the reaction, the total volume of fluid (within V) in which the reactive species are found together is always less than vAB" Substituting (2.8a) and (2.8b) into (2.8c) and combining .the result with (2.3), we obtain the rate equation (2.11) -26- on making use of Eqs. (2.4a,b) we can reduce (2.11) to the final form (2 .12) where (2 . 13) (2.14) (2.15) The parameters sA and sB' which we call the mixing parameters~ are the ratios of the quantities of A1 and BI in the domain vAB to the total quantities of these species present in V. Thus, sA and sB are strictly determined by the properties of the turbulence and molecular diffusion because they are defined in terms of the inert materials AI and B1. The parameter A, which we call the r eaction par ameter , is the ratio of the correlation coefficient of the concentrations of the reactive species in the domain vAB to that of the inert materials AI and B1 in the same region. In contrast to sA and sB' the chemistry plays a role in determining the functional form of A. Although rAB' AI, and BI are quantities that are readily measurable, the parameters A, sA' and sB are not and thus Eq . (2.12) is unsolvable until appropriate expressions for these parameters have been found. It is apparent, therefore, that we have transformed the closure -27problem associated with Eq. (2.3) into one of finding expressions for the mixing and reaction parameters sA, s8 , and \ . The latter task would seem to be more tractable because as it is formulated, the effects of the chemistry and those of the turbulence are to a large extent segregated. This permits a much easier use of physical insight in de- veloping expressions for the unknown parameters and it facilitates the use of the resulting rate equations in applied problems. In the next section we discuss in detail the determination of sA, s8 , and\. First, we note some special cases. premixed, sA = s8 = 1. If the reactants are fully Therefore (2.12) simplifies to (2.16) Eq. (2.16) is also the asymptotic form of (2.12) in the case t ➔ 00 • In practice, premixed reactants occur primarily in two cases. Firstly, when one species interacts with itself, e.g. A+ A ➔ products, (2.16) becomes dA _ 2 (2.17) dt - - kH A A where rA = Ai/~. Secondly, when A and Bare mixed under a condition un- favorable to reaction such as lack of radiation energy (in photochemical reactions), low temperature in an endothermal reaction, etc. the governing equation is again of the form of (2. 16). scopic mixing is considered complete. In both cases, the macro- If the microscopic mixing is also complete, then the reaction is much slower than the mixing process, and -28r and ss approach unity. In addition, we shall see subsequently r AB' 'A' that A also approaches unity. Thus, in this case dA = - kAB dt ( 2. 18) which is the familiar expression generally used for turbulent secondorder reaction. II.2 Estimation of Parameters Since both the mixing and reaction parameters are defined in terms of characteristics of the diffusion zone vAB' an understanding of the processes that control the generation and fate of this region is essential to the formulation of meaningful expressions for sA' s 8 , and A. Of course, the precise details of the evolution of the diffu- sion zone are not known, but sufficient knowledge of a heuristic nature is available (due mainly to the work of Batchelor (1959) and others) that a model of vAB adequate for the purposes of estimating sA' s8 , and A can be developed. The heuristic knowledge of which we speak is best described using an example. Suppose that marked fluid of molecular diffusivity Vis injected into a turbulent fluid of kinematic viscosity v in which turbulent energy is being dissipated at a rate E • Velocity fluctuations of L spatial scales that are large compared to (v 3/ E)~ serve to distort the initial cloud of marked fluid by pulling and squeezing it into convoluted sheets and ribbons of ever decreasing thickness. At scales small -29compared to (v 3/E)¼, the velocity field appears as a pure straining motion and consequently, sheets and ribbons of marked fluid produced by larger scale velocity perturbations are squeezed still further . If the molecular diffusivity of the marked fluid were zero, the stretching and squeezing of the sheets of marked fluid would continue indefinitely and mixing would occur only in the limit of infinitely large times. However, if vis finite, then according to Batchelor (1959), once the sheet widths 1 approach a thickness comparable to (vV 2 /E)~ the concentration gradients normal to the sheet are sufficiently large that molecular fluxes of marked fluid just balance the compression of the sheet by the straining motion . At points in the fluid where this occurs, thorough mixing of the marked and convecting fluids results, and hence, by our definition, portions of the diffusion zone vAB are created. At scales larger compared to 1 (vV 2 /E)~ concentration gradients are too weak to cause significant mixing and thus the generation of vAB is confined primarily to those portions 1 of the marked fluid where the sheet thickness is comparable to (vV 2 /E)~ . Based on these observations, we conclude that the evolution of vAB can be modeled in one dimension using concentration pulses of widths and Y. separations of the order of (vV 2/E) 4. We shall subsequently derive expressions for sA, s8 , and\ based on such a model. -30- 2 1 The Mixing Parameters II.. By definition sA is the ratio of the quantity of inert A brought into vAB by turbulent mixing and molecular diffusion to the total amount of A present~ Hence, it is a measure of the mixing efficiency. In principle, sA and sB can be determined by measurements of the mixi ng of inert species A1 and B1 in the turbulent fluid. When this type of measurement is not available, an approximation may be developed . If rAB is known, rA and r 8 in general are also known (see Appendix A). Hence , rA and r may be used as additional available information. 8 Let us define - l A A2 ) -2 rA = Ay/Ay = V (µAvAAY + µavABaI /AI (2.19) where µ -- -l A A2 I l µa - a2 I l VA f A2I dv l VAB f (2.20) VA a2I dv (2.21) VAB with similar definitions for r 8 , µB' and µB. Substituting (2 .4ab) and (2.14) into (2.19), we obtain A r = SA r _V_ [ µAAI ( l A vAB A a 1 - SA ) + µ~SA ] (2.22) -31(2.23) with a similar expression for rA 8;r 8 . Although (2.23) gives sA in terms of the measurable quantity rA 8/rA' the expression involves µ1 , µA' and µa' which are unknown. (The closure problem has, of course, not been resolved; it is implicit in the definitions of the second-order correlations µ1, µa' etc.) It is at this point that an approximation must be made to eliminate the closure problem. This is achieved here by considering some special cases. Let us say that in the first case each fluid element has a flat concentration distribution. Imagine that we measure the concentra- tions of A and B along the plane PP' in Figure 2.l(b) and that we observe a distribution as shown in case l of Table 2.1 (p.36) In reality, molecular diffusion will obliterate sharp discontinuities in concentrations. However, because sA and s8 represent the ratios of the amounts of inert A and Bin vAB to those in the entire system, the details of the concentration profiles within vAB are not of prime importance. We therefore assert that this approximation of concentration profiles is good, particularly for systems where v >> V. In this case, µa= µA= µb = µB = µI = l, and a1 = A1 and b1 = B1. Upon using these approximations in (2.23) we obtain -32- (2.24) When the reactants are non-premixed, sA and s8 start with a value of zero. rAB in this case is also zero. right prediction. In the limit oft ➔ 00 , reaction occurs in a confined environment. also go Hence (2.24) gives the sA and sB go to unity if the rAB' rA and r 8 in this limit to one, hence (2.24) also gives the correct prediction. How- ever, we expect that (2.24) slightly underestimates the true s 's. This is because when all the reactants are brought into the diffusion zones (tA = ss = 1), there still can be concentration fluctuations. Therefore, rAB approaches unity slower than s does, and rAB;r 8 (or rAB/rA) may be smaller than the true sA (or s 8 ). When the reaction is sufficiently fast, this estimation may predict a reaction rate slower than the true value. This will be seen later . Another important case one often encounters is that in which one of the reactants is emitted into an environment containing a uniform concentration of the other. Let us say that concentration B1 has a uni- form distribution as depicted in case 2 of Table 2.1. We have, in this case, A1B1 = B1A1. Consequently, rB = rAB = l, and µ8 = µb =µI= l. Using these relationships in (2.23), we obtain sA = l The correlation coefficient µa remains to be determined. (2.25) -33- We note that µa is a function of concentration statistics within the diffusion zone vAB" Since molecular diffusion is the most effective means of transport in the diffusion zone, µa is largely determined by molecular diffusion. We therefore propose that µa can be estimated to good approximation using the soluti?n of the molecular diffusion equation a2 a1 aa 1 _ (2.26) at - VA~ with initial condition (2. 27) otherwise. The width o is chosen to be the basic length scale of diffusion zones; i.e., (vV 2/£)¼. (This point is discussed further in the next subsection where we use essentially the same technique to estimate the reaction parameter A.) The solution of (2.26) and (2.27) is (see for example, Cars law and Jaeger, 1959). x+ - o x- - o 2l ) ] - -2l- ) - erf ( ( 2/VA t . 2IVAt a = - 0 2 [ erf Figure 2.2 shows a plot ofµ a (2.28) as a function of nondimensionalized time t VA/6 2 derived from the definition (2.21) of µa and (2.28). According to these results, for times larger than about 10 4 (6 2 /VA) sA and s8 are given by (2.29) --34- 10-5 10-4 10-3 10-2 10- 1 10° 10 1 10 2 10 3 10 4 10 5 2 T /C 2i /0) Figure 2.2 Temporal behavior of auto-correlation coefficient µa in a diffusion zone. -35- It turns out that these expressions are applicable to a wide range of problems of applied interest. The approximate formulas for sA and s8 developed in this section and their possible applications are summari zed in Table 2.1. As discussed before, these approximations may result in slight underestimates of reaction rates. It is hoped that a more universal formula for the mixing parameters will be developed in the future. For practical interest, however, these approximations appear to be useful and reasonably accurate (cf. Sections II . 3 and III . 4). (lb) -~-(2b) --B--- (2a) ,------, A r,B - rA Ila with µa given in Fig. 2.2 or taken as unity. r,A = l, B rAB rAB r,A = -r- ' ts = -r- Approximation Formulas Possible Applications One reactant is introduced into a uniform mixture of another reactant. In atmosphere, 11a may be taken as l. Liquid-phase reactors with non-premixed feeds. Better when v/V» l. Summary of Approximation Formulas for the Mixing Parameters --,-:--,1 ITT .---I (la) Di iD ,--, --------bI Legend: aI, t > 0 t = 0 Model Concentration Profile Table 2.1 I w I 0) -37- II.2.2 The Reaction Parameter The reaction parameter A is defined as the ratio of the concentration correlation coefficient in the reactive case to that in the inert case and depends on both molecular di ff us ·ion and chemical reaction . Thus, A cannot be determined solely from inert mix i ng experiments and must be predicted. We note that A is a function only of concentration statistics within the diffusion zone vAB' Since turbulence can only k* bring A and B within a distance of the order of (vV 2 /E) 4 of each other, only molecular diffusion is capable of producing the intimate mixing required for reaction. Let us first examine the limiting behavior of A. When the rate of reaction is much slower than that of molecular diffusion, the reactant concentrations do not differ appreciably from AI and BI. this case, A is approximately unity. In On the other hand, if the reaction rate is much faster than the mixing rate, reaction always tends to exhaust the reactants as soon as they are brought into the diffusion zone. Then, the concentrations of either A or B or both will be zero in vAB' In this limit A ➔ 0. Regardless of how fast the reaction is, there always will exist a period of time , 0 .::_ t .::_ TR' where TR is an initial reaction time scale, during which the concentrations of the reactants will not have k *This length scale is (vV 2 /E)~ if v>>D and (V 3/E) 4 if v<<D. (Batchelor, 1959) l -38- changed appreciably. In this period A; AI and B; BI and\= l. As the reaction and the diffusion proceed to completion, concentration gradients diminish and bothµ and should again approach unity. µI approach unity . Thus, as t ➔ In the range 'R < t < 00 00 , \ we expect\< l. For the purpose of estimating\ we propose to describe the concentration distributions in the diffusion zone by the following idealized representation (see Figure 2.3). Thin sheets of fluid containing A and Bare initially separated by a distance of one sheet width. The sheet width o is taken as the closest distance turbulence is capable of bringing k two fluid elements together, i.e. (vV 2 /s) 4 (assuming V = l/2(VA + VB) if the diffusivities are not equal). Let a, b, a 1 and b1 denote the instan- taneous concentrations of A and Bin the reactive and inert cases, respectively. The equations governing these concentrations are aa - V a2a kb atAw- a (2.25) (2.26) aa 1 a2aI at = VA ~ ( 2. 27) abI _ a2bI at- VB~ (2.28) with initial and boundary conditions 3/2 8 2 x 2 - l /2 o otherwise (2.29) ; --- ... CJ. ,.. ., ,, b t- S-t ~ --t- b 1 I I I I I I I I I I I I a.o r-----, bo ( b) r- 6-, a.o ~---- bo f*t Premixed Reactants Ol, b ~%1 t>o t-o Figure 2. 3 Schematic concentration profiles in estimation of\ (~ ) Non- -Premixed Reactants .,.. ., ,,. ,. ,. ... "t > 0 t: 0 I w I.D I -401/2 o < X < 3/2 o (2.30) otherwise a(x,t) = b(x,t) = 0 ± oo X+ (2.31) Defining the following dimensionless parameters, a* = ~ ao b* = Q_ bo a = :I 1 o b* = ~ I bo x* = -X0 ko2 a a 0 = ----=-- VA 8 (2.25)-(2.31) become aa*· = ~ a2 a* - aa a*b*• -at* ax*. :. µ (2.32) ab* a2 b* at*= Y ax* z - aa*b* (2.33) (2.34) (2.35) 3 011 a*(x*,0) = a1(x*,0) -~ j_ b*(x*,O) = bJ(x*,O) { a*(x*,t*) = b*(x*,t*) = 0 l - 2 2- x* 2- - 2 otherwise l2 -< x* -< l2 (2.36) ( 2. 37) otherwise x* + ± oo (2.38) -41 The solutions for al and bi are a* (x*,t*) = l 11erf (x*+ 3/ 2 ) - erf ( x*+l/ 2 )] 2 1 2/F" 2~ (2.39) b* (x*, t*) = l2 [erf ( x*- l 12 ) - erf ( x*- 3/ 2 ) ] 1 2/yt* 2/yf* (2.40) L Determination of a* and b* requires numerical solution of (2.32) and (2.33), the details of which are given in Appendix B. Once a*(x*,t*) and b*(x*,t*) have been determined, A is computed as follows. a*b*dx* aidx* bidx* LAB LAB LAB f f f (2.41) A -- I a*b*dx* L AB I I f a*dx* LAB f b*dx* LAB The boundaries of the one-dimensional diffusion zone LAB are chosen as the points where al or bj has a value of 0.001. (See Appendix B for details.) The behavior of A as a function oft* is shown in Figures 2.42.6 for the combinations of the parameters a, s, and y summarized in Table 2.2. Table 2. 2 Combinations of Parameters for A vs. t* Figure No. a s y 2.4 0.1-10 4 1 l 2.5 10 3 0.1-10 l 2.6 10 3 l 0.1-10 10-II I 10-2 10-3 2 A 10- 10-1 10° I I "'" I Figure 2.4 10-1 I 111111 ........ 10 1 I 1111111 l)At/b I 2 I I 102 I 111111 I 103 I 1111111 at 6 =y = 1. Temporal behavior of A for various values of a 10° I 1111111 I 4 3 10 - - - - - - 10 - - - - - 100. - - - - - - - - - - - - ,o. 1.0 O. I ct=- I 11 I ~ I N F I 10-1 I I 111111 I ' 101 - /)At/ ?.> 2 100 I I I 11111 - I I 102 1111111 103 I I 1111111 1. 10. 0. I Figure 2. 5 Temporal behavior of A for various values of~ at a = 1000 and y = 1. 1o-" , 10-2 10-3 2 A 10- 10-1 100 ~; I I +:> w r r- ~ r 1o-" I 10-3 2 A 10- 10-1 100 I 10-1 I I I I 1111 I I 10 1 I I I I 111 I ~At/ ~Z 10° I I I I 1111 .._ I ! 102 I I I I I 11 I 103 I I I I 1111 0.1 1. 10. Figure 2.6 Temporal behavior of A for various values of y at a= 1000 and j = 1. 10-2 I ~ ~\ '(: I -$::> -$::> I -45The temporal behavior of A is consistent with expectations, i.e. at t* = 0, A= l; fort*> 0, A decreases; and as t* becomes large, A again approaches 1. This limiting behavior is not shown in Figures 2.4 and 2.6 because it happened much later in the case B = l. In addition, the variation of A vs. t* with the parameters a, B, and y can be discussed on a physical basis. For a fl, A= 1 for all t*. When a~ l, reaction is slow com- pared to molecular diffusion, the concentration profiles differ little from the inert profiles, and A is always close to unity. As a increases beyond order unity, the reaction speed increases and A decreases from its initial value of 1. A is found to reach a minimum value over some period of time, the minimum value decreasing as a increases. Physically, as a increases, the reaction rate increases relative to the molecular diffusion rate, reactant concentrations within the diffusion zone decrease and become appreciably smaller than inert concentrations, and A decreases. (Figure 2. 4) When B = 1 (stoichiometric feed),A remains at its minimum value for a period of time of the order of 10 4 (8 2 /VA) (not shown in the figures). Forst l, the period of time over which A is at its minimum value is shorter than that for B = 1. This behavior occurs because in the case of non-stoichiometric feed, the reaction goes to completion faster than when B = l. The Amin (the minimum value of A) remains essen- tially the same in the case of B > 1 and the case of B = 1. For the -46- case _when s < 1, Amin differs considerably from the value at s = 1. closer examination reveals that Amin in the case of s < 1 has the same ko 2 b value as the Amin for the case where a= v~ , s ,:"._ 1. This is because whens< l, reactant Bis present in a lesser amount and the reaction is limited by the time scale (kb 0 )-l rather than (ka 0 )-l. (Compare the Amin . in the curve of a= 1000, s = 0.1 in Figure 2.5 with that in the curve of a= 100, s = 1 in Figure 2.4.) The effect of y is shown in Figure 2.6. When y < 1 (VB< VA) the behavior of A does not deviate appreciably from that for the case when y = l. But when y > l (V 8 > VA)' A reaches a minimum value which ko 2 a0 would be reached for the case where a= v:-(Compare the curve for B a= 1000, y = 10 in Figure 2.6 with the curve for a= 100,y = 1 in Figure 2.4.) It therefore becomes obvious that the controlling time scale for molecular diffusion is o2 /VB if VB> VA and is o2 /VA if VA> VB. (Note that A decays to its minimum in about 0.1 (o 2 /VA) = o2 /VB in the case of y = 10, rather than o2 /VA as in the case of y = l .) Of most interest in the temporal behavior of A is that region of time when A is at its minimum value. The early decay of A from 1 to Amin occurs on a time scale of the order of o2 /VA (if VA> VB). For a fast reaction (diffusion-limited) the reaction does not substantially commence until a time the order of 6 2 /VA, since this is the time required to bring the reactants into contact. For a slow reaction, the time scale o2 /VA required for A to decay to Amin is short compared to the reaction time scale. In addition, in the basic equation (2.13), A and rAB occurs -47- as a product and since rAB starts at a small value, the prediction should be insensitive to A at small times. increasing from Amin toward l as t* ➔ On the other hand, when A starts 00 , the reaction has essentially gone to completion, and this region is of little concern. of A of most importance is the minimum value Amin' Thus the value Based on the foregoing • arguments, we will take A to be a constant, Amin' for use in predicting the rate of decay of the reactive species. Based on the values of Amin in Figure 2.4 (s = y = l), it is found that Amm . can be described by a simple expression as a function of a: A =A . mm = .--:l=--=--=l+O. l6a (2.41) As we have already seen, the effect of sandy on Amm . is minor for s -> l and y .::_ l. Therefore we shall define a as the ratio of the shortest diffusion time scale to the longest reaction time scale. That is ko 2min[a 0 , b0 J a = max[VA,VB] (2.42) where a0 , b0 are feed concentrations of the reactive species A and B, respectively. Physically such a definition is justified because the diffusion process is dominated by the species with the larger diffusivity, while a second-order reaction is limited by the species with the smaller concentration. If reaction is much faster than diffusion, i.e. a ➔ a\= l/0.16. Thus, in this case (2.13) becomes 00 , -48- (2.43) which, as might be expected, is independent of the rate constant k. Physically, in this case, the reaction rate is solely controlled by the rate of turbulent mixing and molecular diffusion . Only inert mixing data is required in (2.43). (see Appendix A). A similar result was obtained by Toor (1962) In that model it is predicted that the rate of conver- sion of a very rapid reaction in a tubular reactor is dependent on inert concentration statistics only. In the next section we shall compare the predictions of Toor's and other models with the present one. It should be pointed out that in turbulence, quantities such as o are often known only in terms of their orders of magnitude. Thus one might wonder how sensitive the prediction of the model is to the value of a. It will be shown later (Figures 2.13 to 2.18) that an estimate of a accurate to within one order of magnitude is sufficient to predict laboratory data within experimental error. For the case of premixed reactants, A is calculated by solving (2.25-28) with the set of initial conditions depicted in Figure 2.3b, p.39. The temporal behavior of A in this case shows two distinct patterns. For the case of s = l, A stays atunityfor all times; but for cases wheres fl, A decays to the same minimum value as its non-premixed counterpart. Physically these behaviors can be understood with the aid of Figure 2.3(b).In the case of s = l, which is equivalent to the case -49- of A+ A ➔ products, bothµ and µI are autocorrelations of Gaussian type distributions. Consequently,µ= µI and A = l . When the feed ratio is not stoichiometric, the reactant in excess concentration depletes the other reactant and creates a valley in its concentrat i on profile (F i gure 2. 3b) • The concentration correl ation t he refore becomes si mil ar to the non-premixed case, and A shares the same minimum value . JI.3 Application of the Diffusion Zone Model to Laborato ry Data An extensive experi men t al program has been conducted over a long period in the laboratories of Toor (Vassilatos and Toor, 1965 ; Mao, 1969) and Brodkey (McKelvey , 1968; Zakanycz, 1971) on chemical reactions occurring in turbulent mixing. Both groups have employed a multi- jet reactor originally described by Vassilatos and Toor (1965) . The reactor and the nature of the experiments are described in Appendix A. In essence, in the reactor alternate jets of acid and base are ejected from a plate containing numerous small holes into a tube, creating efficient turbulent mixing . Conditions are such that the mean concentrations of A and B vary only with distance from the entrance . AI and B1 are constant throughout the reactor.) (In the inert case The reactor affords the opportunity to measure the rate of reaction between two species in the presence of rather well-characterized turbulent mixing. In this section we apply the diffusion zone model to several experiments reported on the two multijet reactors in order to calculate the degree of conversion with d.1st ance down the reactors . -50- Since Ar and Br are constant in the reactor, it is convenient to let FA= A/Ar and F8 = B/B 1, in which case (2.12) becomes (2.44) Since the reactions studied in the multijet reacto r were of the general form A+ nB ➔ products, we can write (2 A5) or (2.46) Substituting (2.46) into (2.44 ) gives (2.47) wheres= Br/nA 1 and the relationship dx/dt = U(x) has been used to convert time to distance x down the reactor. U(x) is the mean axial velocity at position x and is discussed in Appendix A. Values of rAB' rA, and r 8 required to evaluate ~A and ~B by (2.24) were obtained from conversion data for very rapid reactions in the case of s = 1. (see Appendix A.) late A. Equation (2.41) was used to calcu- Table 2.3 lists the experiments that we will consider. cally two types of reactions were studied. Basi- The first are those of .03331 .03405 .03387 .03357 .03363 .03335 .02910 .02830 . 02720 8320 8320 l .4xl0 7 NaOH NaOH OH- CO 2 CO 2 Nitril otriacetic Acid 3.0 1. 5 l.O l.O l.5 3.0 2 l l.O l. 5 3.0 l. 26 l.39 l. 71 2.53 3.88 s 2 2 n from Yieh (1970) * Reynolds number is based on nozzle measurements. .243 .243 .241 .241 .246 .245 .075 .067 .070 .071 .078 / 177. 1 172.3 165. 5 r a= kna 0 o2 /V, where v = VA= 113 = 2.3xlo- 5 cm/sec, o = 10- • cm .00699 .00625 .00655 .00658 . 00724 12400 NaOH CO 2 ao g-mole/£ £.,/g-mole-sec B k Mao (1969) Re= 2430 Mao (1969) Re= 2430 Mao (1969) Re= 1460 Vassilatos and Toor (1965) Re= 3700* Reference Experimental Conditions for Data Used to Evaluate Diffusion Zone Model A Species Table 2.3 I U7 I I-' -52so-ca , led intermediate rates, having a value of a of 0.01-1. type are verY fast reactions (a>> l) . The second The laboratory data reported of this type ar e independent of the rate constant k (diffusion-limited) , and some of them have been used to compute rAB etc. Slow reactions (a<<<l) are not of interest to us because they can be represented by the homogeneous law ( eq. ( l. 16)). In the calculation (2.46) was solved numerically to predict FA as a function of distance from the reactor inlet . The results of the calculations for the experiments listed in Table 2.3 are given in Figures 2.7-2.20. In each case, information on inert concentration correla- tions was obtained using Toor's model, as described in Appendix A. Hence, any errors inherent in rAB' rA, and r 8 will be carried into our calculation. In general, the predictions agree more closely with the data of Vassilatos and Toor (1965), Figures 2.7-2 .11, than with those of Mao (1969), Figures 2.12-2.20. In most cases it was found that the predicted concentration of A was somewhat higher than the observed for the data of Mao (1969). Since negative values of F(x) were reported in several sets of Mao's data, it is possible that his measurements were in error (Figures 2.14, 17, 19, 20). In addition, as we have indicated, the approxi- mation (2.24) will generally give a slower decay of the species and may have been the source of error. One of the most important features of our diffusion zone (DZ) model revealed in these calculations is the effect of inhomogeneous mixing on the rate of reaction. Comparison with the predictions of the -53- 0 • -----,-----.---r-----r-----.---.----------.---, DRTR CJF VRSSILRTCJS -- D . Z . MCJDEL ro 0 ® • (D ,..... 0 • (!) X "-,/ . LL ::t' 0 (\J 0 D 0 (!) • (!) . , _ _ _ - - - - - - - ' - - - - . . . J . . - - - L - - - _ _ . _ __ _---L.-_ _....___ 0.0 1.0 2.0 __,__ _.....J 3.0 X,CM K = 0. 12L!E 05 LIT /MCJL-SEC RO= 0 . 00699 MCJL/LIT BETR = 1.26 Figure 2.7 Comparison of DZ model prediction with reactor data. 4.0 -54- D · - - - - . - - - - . - -.--, r---r--~--r---.---, DRTR □ F VRSSILRTGS - 0. Z. MGDEL ~ 00 I D (.J) I D ......,, l.J_ :::1' I (!) D (!) N (!) I D (!) 0 I o..___ __.__ __.___ _..___ __.___ ___.___ _..____ __,__ ___, 0.0 1.0 2.0 3.0 X,CM 0. 12L!E 05 LIT /MGL-SEC RO::: 0.00625 MGL/LIT BETA ::: 1 .39 K ::: Fi~ure 2.3 Comparison of DZ model prediction with reactor data. 4.0 -55- 0 I ..----------,---,-----~----,----~---.---~-- DRTR □ F VRSSILRT □ S - D. Z. M □ DEL 00 © I 0 (D I 0 X ..._,, LJ..... ::f' I 0 C\I I 0 0 oL_ _L_ 0.0 _j__ 1 .o _L_ _J~_::::=t=:=:~===::=±:==-.J 2.0 3.0 X,CM K = 0.12L!E 05 LIT/M □ L-SEC RO= 0.00655 BETR = 1 . 71 M □ L/LIT Figure 2.9 Comparison of DZ model prediction with reactor data. 4.0 -56- D ......· ~ - - - - r - - . - - - - - - , - - - - , - - - , - - - - - - , - - - - - - r - - - ORTA l)f VRSSILRTl)S - 0. Z. Ml)OEL 00 © I D --- KRTTRN $ ADLER( 1967) (D I - ·- Hl)Ml)G. Sl)LN . ,..... D X '-../ LL ::r I D I I ' I \ ~ (!) \ I I ' ('\J I D I I I \ \ I \ \ 0.0 \: 1 .o 2.0 3.0 X,CM K = 0.12LlE 05 RO= 0.00658 BETR = 2.52 LIT/Ml)L-SEC Ml)L/LlT Figure 2.10 Comparison of DZ model prediction with reactor data and other models. 4.0 -57- D • -------r----,-------r---r---,------,-------r-------, ORTA ~F VRSSILRT~S - D. Z. M~DEL 00 D © • . (D D X '--/ lL::::f' D • . 0J D D . oL.-_ 0.0 _L_ _.J__---2:.J.---====:::±::==~====1:==--..J....-__J 1 .o 2.0 3.0 X,CM K == 0.121.!E 05 LIT/M~L-SEC RO== 0.00721.! M~L/LIT BETA= 3.88 Figure 2.11 Comparison of DZ model prediction with reactor data. 4.0 -58- homo gen eous 1aw, i . e. dA U( X) dx kAB are shown in several cases. (2.48) (Figures 2.10, 15-16). It is seen that the inhomogeneity in mixing has inhibited the reaction considerably. furthermore, the DZ-model predictions are very sensitive to the mixing condition in the reactor. This is seen particularly in Figures 2.12-14 where a break in calculated FA occurs around x = 0.2 cm. Comparing these results with the data used for rAB' rA' etc (Figure A.3), we conclude that this break is solely a reflection of the mixing condition. Physi- cally, this breaking point is where the jets from reactor head meet. Considerable backflow is induced here and the mixing is temporarily slowed. This effect is not shown in Vassilatos' data because data were not taken close enough to the head region. With these data it was necessary to extrapolate the data back into the region of the reactor head (see Fig. A.2). The data gathered by Mao (1969) from the laminar jets (Re= 1460) (Fig. A.3) reveal a much slower rate of mixing than the data from the turbulent jets (Re= 2430) (Fig. A.4). All these characteristics in the mixing process are reflected in the DZ model Predictions. The sensitivity of the model prediction to estimates of parameters such as a is shown in Figures 2.13 and 2.18. In the case of inter- mediate rate reactions, Figure 2.13 indicates that variations of an order of magnitude in the value of a result in negligible changes in the -59- 0 · - - - - r - - , -- - - - , - - - - r - - - - - - ; - - - - , - - - - , - -- - - , ...... . co 0 DRTR CJF MRCJ . RE=lij60 - 0. Z. MCJDEL . \ 0 \ \ \ \ \ \ 0 MRCJ $ TOCJRC1970) MR □ $ TCJCJRC1971 \ . (D \ \ X 0 . (\J 0 0 . o.__ __,__ _ 0.0 __J___ _ _ , __ _ _ J _ _ _ , __ _ _ . _ _ _J _ __ 100 200 3.0 XPCM K = 0.832E OLJ RO= 0.03331 BETR = 1.00 Figure 2. 12 LIT/MOL-SEC MOL/LIT Compar ison of DZ model prediction with reactor data and other models. _, 4.0 -60 - 0 · - - - - - , - - - - - - , - -- , - - - - ,,-----,-- - - , - ---,-------, ORTR OF MRO. RE=l~GO - 0. Z. MLlDEL ~ ro. 0 i 2 3 . (D .,... 0 X d.= 2A6 d- = 0. 214- 6 (Predi.ded) d = 0,02'+6 '-../ l.J.... ::t' . 0 . (\J 0 0 a L_l___ _j__ 0.0 1.0 _l_~t::::::::::::f=::======4=......JL...-d 2.0 3.0 4 .0 X9CM K = 0. 832E OW LIT/MOL-SEC RO= 0 . 03W05 MOL/LIT BETR = 1 .50 Figure 2. 13 Comparison of DZ model prediction with re ac t or data . The model is not very sensitive to estimates of its parameter ~- -61- 0 • - - - - . - - , - - - - -r--....... . DRTR OF MRO r- RE=lWGO - 0. Z . MODEL ~ OJ 0 - , - -----,-- -,---- - - - r - - - , . (D 0 X 0 . . 0J C) 0 . a t___ 0.0 _L_ _..1__-----,,i.,_---=:::::::t:=====6=~.L..,_----'----1 1 .0 3 .0 K = 0.832E Qij LIT/MOL-SEC RO= 0.03387 MOL/LIT BETR = 3 .00 Figure 2. 14 Comparison of DZ model prediction with reactor data. L!.O -62. ted mean concentration. The prediction of A in the case of very pre d1c fast reactions is more sensitive to the value of a. However, in this . ·t the reaction is diffusion-controlled and the influence of the 11m1 reaction parameter A is small . We therefore see that again variations of one order of magnitude in o: change the convers"ion rate by only a small factor (Fig. 2.18). Some previous work is discussed here. Kattan and Adler (1967) developed a so-called stochastic mixing model in which fluid elements containing reactive species are assumed to interact through a random process of coagulation and redispersion. The rate of this coagulation- redispersion process is adjusted to match the inert mixing data. Chemi- cal reaction is incorporated into the model by changing the concentration of the species in each fluid element in accordance with the governing kinetics. Kattan and Adler (1967) compared the prediction of their model with Vassilatos' data and obtained results such as those plotted in Figure 2.10. It is seen that their model cannot describe the entire range of the chemical conversion process . The model of Mao and Toor (1970) (see also Mao, 1969) divides the fluid into slabs of thickness o which are alternately arranged with reactants A and Bat constant initial concentrations. Since only mo- lecular diffusion is allowed between the slabs, the set of model equations is very similar to (2.25-28). The thickness o is first adjusted so that the solution of the system fits the data of a very rapid reactio • the same turbulence . • n in The resulting value of o is then used to -63- predict conversions of other reactions in the same turbulence. Rao and Edwards (1971) have shown that the model of Mao and Toor (1970) is in essence the same as the model of Kattan and Adler (1967). One of the calculations of Mao and Toor (1970) is given in Figure 2.12 in comparison with our DZ-model. The deviation of their model from the observed data is quite obvious. Mao and Toor (1971) later proposed a second model in which the reactor is divided into two regions. The first region extends from the head of the reactor to the point where the jets meet. Because of the considerable backflow in this region, they treat this domain as a homogeneous mixer, (2.48). The second region covers the re- mainder of the reactor and is modeled using the hypothesis of Toor (1969) A'B 1 = A B1. Under this hypothesis, 1 udA (2.49) dx where A1B1 is derived from very rapid reaction data as in Appendix A, and U is the constant bulk average velocity of the flow. One of the cal- culations of Mao and Toor (1971) is given in Figure 2.12. The prediction is not in good agreement with experimental data. Moreover, since the division between the two regions of the reactor cannot be determined a Priori, this model is awkward for many purposes. Yieh (1970) discovered that the model of Mao and Toor (1971), eq. (2.49),is improved if U is measured as a function of x rather th an as a constant bulk velocity. Yieh assumed that the product concen- tration in a multijet reactor is uniformly distributed across the -64reactor. That is, if C is the instantaneous concentration of the product, C' = C - C = 0 (2.50) Since mass conservation requires that A+ B + C = 0, we must have A' + B' = - C' = 0. This result causes the equation governing A'S' (eq. (l.15)) to reduce to a form that is independent of the reaction kinetics. Yieh (1970) thus concluded that A'B' = A'B' I I and arrived at (2.49). Yieh's predictions based on the same velocity profile U(x) as we used give good agreement with observations. ures 2.15-16). (Fig- However, assumption (2.50) is not generally valid and therefore Yieh's model is greatly limited in its scope of applications. Applications of the diffusion zone model to very rapid reactions (a>>l) are shown in Figures 2.18-20. Relatively good agreement between the predicted and observed conversion rates results in all three cases. Toor (1962) derived a model for very rapid reactions in - 65 - 0 ......· - - - - . . . . . - - - - . - - - - - - - r - - - - r - - - - - . - - - ~ - - - - - . - - - - . (X) I (!) \ D \ (!) "' . (D D DRTR CJF MR □ . RE=2ij30 - D. Z . MCJDEL ~ \ \ (!) YI EHC 1970) \ \ HCJMCJG . SCJLN. \ \ \ X \ \ '-.../ \ . " LL :::j< D \ \ " '-;: ~ .., . \ (\J D --·-·- · - D Q . . . _ _ _ __ · -. -----L-_ ____.___ _. . _ __ _ . . . . j ._ _ 0 .0 1.0 - ·- - · - - - _ . . . . j ._ __ _ _ . __ _~ - - - - 3.0 K = 0.832E Qij LIT/MCJL-SEC RO= 0 .03357 MCJL/LIT BETR = 1 .00 Figure 2.15 Comparison of DZ model prediction with reactor data and other models. 4 .0 --66- 0 • ~ - - - - . - - - - . -- DRfR OF MRO. RE=2ij30 - 0 . Z. MODEL ~ (XJ . 0 Q) \ \ w. -,-----r-------r------r---~----, \ Q) ' '\ \ Y IEHC 1970) \ HCJMOG \ 0 \ b SCJLN. \ \ \ \ \ \ \ 0 \ \ \ \ (\J . 0 LO K = 0.832E OLJ RO= 0.03363 3.0 LIT/MOL-SEC MOL/LIT BETR = 1.50 Figure 2.16 Comparison of DZ model prediction with reactor data and other models. 4.0 -67- 0 . ~---.--~--,-----,---,---,-------.----, DRTR OF MR □. RE=2ij30 - 0. Z,, M[jOEL ~ co. 0 . (1J 0 r'- x ...__,, lJ_ :::r . 0 . (\J 0 0 Q) . Q L __ _j__ ___L__ _L___:=::=:=±==='=--L----1.._ 0.0 3.0 K = 0.832E 0~ LIT/MOL-SEC RO= 0.03335 MOL/LIT BETR = 3.00 Figure 2.17 Comparison of DZ model prediction with reactor data. __J 4.0 -68which the conversion depends only on inert concentration statistics (Appendix A). Since we have used the experimental data in Figure 2.18 and Toor's model to estimate the values of rAB etc. required in the DZ model, a comparison with Toor's model for this case (B = l) is not possible. However, the comparison of predictions of both models in the case of non-stoichiometric feed are shown in Figures 2.19-20. culations of Toor ' s model were done by Mao (1969). The cal- The difference be- tween the two models seems to be minor. We have therefore demonstrated that the diffusion-zone model provides a realistic description of the rate of decay of chemical species reacting in a given turbulence . This model applies to systems of any feed concentrations, any reaction speed and any turbulence, provided that the inert concentration statistics in that turbulence are given . The DZ model does not require adjustment of any of its parameters and estimate of its parameters can be obtained from turbulence data such as rate of energy dissipation E, etc, that are usually available. We also note that the solutions of the model equations are simple to obtain numerically, in contrast to some currently available models. Having demonstrated the applicability of the DZ model to reactor data, we consider in the next chapter applications to nonlinear chemical processes in the atmospheric boundary layer. -69- 0 · - - - - r - - - . - - - . - - - . - - - - - - - . -- DATR OF MAG. RE=2430 - 0. Z. MODEL . ~ Cl) 0 i . CD ,..... 0 2 3 X '-/ --;---.---.-,, oZ::::1771 0\ = 177.1 ( Predicted) CA = t 7. 71 . IJ_ :::f' 0 . ('\J C) . C) 0.0 1.0 3.0 0. 1LJOE 08 LIT /MOL-SEC RO= 0.02910 MOL/LIT BETA= 1.00 K = Figure 2.18 Comparison of DZ model prediction with reactor data in the case of very rapid reactions. The model is not very sensitive to estimates of its parameter a. -70- 0 . ...... ' . ~ co 0 I . ___:_ TCJCJRC 1962) (D 0 DRTR 0F MR0. RE=2ij30 D. Z. M0DEL X 0 . (\J 0 0 . QL..._ 0.0 ___.L_ ____L__,.,.,____::::i;:;;;::::::=:i:::====:1.--..J..--....J_-_J 1 .o 2.0 3.0 LLO X,CM K = 0. lijOE 08 LIT /M0L-SEC RO= 0.02830 M0L/LIT BETR = 1 .50 Fi gure 2.19 Comparison of DZ model prediction with reactor data of a very rapid reaction. Toor's model was calculated by Mao(l969). -71- 0 ·--~----r-----r---..------.----r-----,..------- . DRTR OF MR □, RE=2~30 - 0. Z. MODEL co 0 ~ ,,,-... 0 0 \ - - - TOOR( 1962) (!) \ \ \ X \(!) '-" \ . \ LL. :::r 0 \ \ \ \ N 0 ~ . \ \ \ \ \ 0 0 . \ (!) '- '- L__ 0.0 _j__ _--1..:----::..::,___....C:====:b::==ad.-----L--...L..-__J 1 .0 (!) 2.0 3.0 Ll.O X,CM K = O.l~OE 08 LIT/MOL-SEC RO= 0.02720 MOL/LIT BETR = 3.00 Figure 2.20 Comparison of DZ model prediction with reactor data of a very rapid reaction. Toor's model was calculated by Mao(1969). -72- Generalization of the Diffusion Zone Model to Three Spatial Jl,4 - Dimensions For conceptual simplicity, the diffusion zone model was developed for the case of a spatially homogeneous sy stem. ,t he model is generalized to these dimensions. In this sect"ion The basic equation govern- ing the ensemble mean concentration of a species A in a turbulent fluid is (see Chapter I) ( l. 9) Assume the process is stationary so that the ensemble average can be· replaced by a time average, A. If K-theory is used to represent the turbulent mass flux and molecular diffusion flux is neglected, we obtain ~ +- at LI; ax . ( K., ~ ax. ) - k AB aA = a ax. 1 , 1 (2.51) where at a given location x, the time average is defined as t 0 +T A(x) = lim } T-+<x> f t A(x,t) dt ( 2. 52) 0 Now let us define TAB' similar to vAB' as the sum of periods of time within the sample period T when A and Bare found together at x. Similarly, TA and TB are the sums of the time periods in T when A and B, respectively, exist alone at x. inert measurements. All these quantities are defined by Thus, we can decompose (2 . 52) into -73- where AI= -l TA A I A dt 1 TAB f A dt a =- - A A with similar definitions for B, B, and b. Also, we let l AB= rµab TAB AA where Likewise we can define the corresponding inert quantities: where µI _rLTAB1 I TAB Substituting the foregoing expressions into (2.51), we see that the local expression for AB shows the same form as its spatially homogeneous counterpart, -74- ( 2. 53) where ~ I'; A - :: bITAB ss:: - - (2.54) BIT (2.55) Combining (2.51) and (2.53) yields the three-dimensional form of the diffusion zone model, ~ + u. aA at 1 ax.1 == _a_ ax.1 (2.56) Similar expressions can easily be obtained for B. The parameters r,A' r, 8 and A defined in (2.54-55) bear the similar significances as their one-dimensional counterparts (2.14-15 ). Intuitively they can be estimated by the same method, even the same formulas, given in Section II.2. We shall not pursue this matter further because the full equation (2 . 56) is rather difficult to apply. Applica- , tion of (2.56) to predict concentration decay of reactive species in a ~th ree-dimensional turbulent motion can be achieved if the inert concen, tration statistics are given in a three-dimensional form . Since inert :concentration statistics as such are not always available and since Kth • eory is not generally valid, we consider in the next chapter an alterna. the DZ model to three-dimensional problems. f tive meth Od of applying -75- III. APPLICATION OF THE DIFFUSION ZONE MODEL TO CHEMICALLY REACTING PLUMES The principal aim of this thesis is to develop means of accounting for the effect of turbulence on the rates of chemical reactions. As we noted in Chapter I, when the characteristic time for chemical reaction is comparable to or shorter than the characteristic time for turbu lent mixing, the overall rate of reaction can be expected to be influenced by the turbulent mixing process . In atmospheric chemical reactions, such a situation is most likely to occur immediately after gases have been released into the atmosphere. The usual case is that the gases, once emitted, react rapidly with other gases present in the ambient air . It is possible that the overall rate of conversion in such a case is strongly dependent on the rate of mixing of the effluent and the ambient air. An important example of this phenomenon is the reaction of nitric oxide and ozone which occurs when NO is emitted into ozone-containing air , The rate constant for this reaction at 25°C is 29 ppm- 1min-l (Johnston and At an o concentration of 0.1 ppm, the time constant for 3 NO decay is 5 seconds. However, observations such as those of Davis et Crosby, 1954). al. (1974) indicate that the time scale for 50% conversion of NO to N0 2 has a value of 12 to 60 minutes in a power plant plume: Thus , it is suspected that in this instance the rate of turbulent mixing may limit speed and extent of conversion • of No to No . the 2 * The------:---be con~ersion of NO to NO in a power plant plume may not , of _course, attributable solely to lhe reaction of NO and o3. All possible reactions must be considered . ----- -76Two major examples of pollutant sources in the vicinity of which chemical reactions may be important are highways and power plants. Because the characterization of turbulence on a highway containing vehicles is extremely difficult, we shall confine our attention to chemical reactions occurring in a plume from an elevated point source. We consider in particular a continuous, elevated point source of A in an ambient atmosphere containing B . The objective of the study is to develop means of predicting the conversion of A as a function of downwind distance and turbulence characteristics. In principle, the steady-state ensemble mean concentrations of two species A and Bin a turbulent flow, when A and B react according to A+ B+ P can be determined by solution of the continuity equations. a<A> -= k axk a <UkA I> - k<AB> axk (3. la) a axk a <UkB'> - k<AB> axk (3.lb) <Uk> where --- A' = A <A> B' = B , u'k = Uk - <Uk> The derivation of these equations was given in Chapter I. If appropri- ate expressions are available for the turbulent flux terms <ukA'> etc. -77and the kinetic term <AB>, (3.la) and (3.lb) can be solved to predict <A> and as functions of location in the plume. As discussed in detail in Chapters I and II, the product moment <AB> cannot always be approximated by <A>, and the diffusion-zone model was developed to represent <AB>. Since the diffusion-zone model requires knowledge of the concentration statistics of inert species, it is necessary to obtain these statistics for the given problem. Consequently, a considerable portion of this chapter is devoted to the generation of inert concentration statistics in atmospheric turbulence. If the necessary inert concentration statistics are known, <AB> can be evaluated by the diffusion-zone model, and (3.la) and (3.lb) can be solved. However, we generally do not require all the information pro- vided by <A> and as functions of x1,x 2 , and x3 . We are often interested only in the overall level of conversion (integrated over the plume) as a function of distance from the point of release. Furthermore, in this study our concern is with concentration correlation terms like <AB> and not with turbulent flux terms. Thus, we seek a form of (3.la,b) that is free of the flux terms <ukA'> and <ukB'>. It might be noted that such a formulation for computing overall chemical reaction rates in plumes has not been available heretofore. -78- JII.l Diffusion Zone Model for a Chemically Reacting Plume Let us consider a continuous, elevated point source in the atmospheric boundary layer with the x-axis aligned in the direction of the mean wind. If the slender plume approximation can be made, i.e. that turbulent diffusion in the x-direction can be neglected relative to advection by the mean flow, (3.la) becomes a<A> a { , A. } a { •A. } k AB ax - = -ay - <V > + -az - <W > - < > (3.2) Since the equation for is similar, we shall consider this equation only. We define an arbitrary area Ac in the y-z plane. If zi is the vertical extent of the atmospheric boundary layer, and Ly is a distance in they-direction from the plume centerline large enough so that the spread of the plume in they-direction does not exceed Ly for all values of x of interest, Ac can be defined as the area bounded by z = 0 to zi and y = - Ly toy= Ly· We now define the spatial averages of <A> and <A 2> over AC as A(x) - l - Ac Ly Jz.1 <A(x,y,z)> dzdy I (3.3) JLYJzi <A2(x,y,z)> dzdy (3.4) -L y o _.2(x) -- Ale A -Ly o - 79- The main quantity of interest to us is A(x)Ac which is the mean value of the total amount of A at location x. As we shall see later, this quantity does not depend on the choice of Ac, neither does its rate of change due to chemical reaction. If there is little wind shear and no surface flux of A or B, (3.2) can be integrated with respect toy and z over Ac to yield ~ dA dx = - ~ k AB (3 .5) L fz.1 y <AB> dzdy (3 . 6) where I-L y o Equation (3.5), together with its counterpart, ~ dB dx k AB (3. 7) describes the overall rate of conversion as a function of downwind distance in the plume. ,-,I computing AB. To solve these equations we must develop a means of The diffusion-zone model developed in Chapter II can be ,v modified to represent AB, and we now carry out the modification . In a stationary process the ensemble average of a variable can be obtained in principle by averaging the instantaneous values of the variable over a sufficiently long period of time. ~ of time be T, then we can express A as Let this period -80- 1 A(x) -- Ac provided that T + dtdzdy (3.8) Thus, A(x) can be envisioned as a space-time aver- 00 • age carried out at any position x. The time averaging is necessary to produce the ensemble average <A(x,y,z)>, and then the space averaging - is required to yield A(x) from <A>. At any x we can envision an infinite ensemble of plume cross sections as depicted in Figure 3.1 . Let S be the hypothetical volume of space defined by - Ly -< y -< Ly , 0 -< z -< z.l and t 0 -< t -< t 0 +T. Then (3.8) can be rewritten as -, I A(x) - S A(S,x) dS (3.9) where dS = dydzdt and S = AcT. At any downwind location x, this hypo- s thetical volume of space Sis defined and used for averaging the concentrations. Within this volume, the diffusion zone can be defined as the space where A and B coexist as inerts. Let sAB be the total volume of the diffusion zones, and let sA and s8 be the total volume of space occupied by A alone and B alone, respectively. Then sAB' sA and 5s, varying with x, are simply the counterparts to vAB' vA and v8 of Chapter II. (See Figure 2.1 .) Once we have defined the diffusion zone, the rest of the formulation will be essentially the same as that in Chapter II. Let the average concentration be decomposed into two parts: - l ~ ~ A = S ( sAAI + sABa ) ( 3. 1Oa) Ac -Plume Cross-section zi~------- •z of the Time History Figure 3.1 SA 0 Sa f::7 ~ SAB Schematic of diffusion zones in the hypothetical space S. • 0 ----------~~ 0-------------~t - L~ o Lj t ti T z l nstantaneous -Plume. Cross-section I (X) I f-J -82~ l A A B =- (S B . S B I +sABb) (3.1Gb) where l A= -I SA A I l a -SAB Ads SA J Ads SAG b = -slAB A J B ds r-.J The concentration product AB is rewritten as (3.11) With (3.9) and (3.11), (3.3) becomes (3.12) Let A1 and B1 be the concentrations of A and B if they were inert, ~ l A t sAB a I ) A AI = S ( SA A ~ l A ( 3. l 3a) A BI = S ( 5s B I+ sAB bI ) (3.13b) and let (3.14) A where a 1, etc are defined similarly to (2.14)-(2.16). -83- Invoking the hypothesis that the average concentration outside the diffusion zone are the same for both inert and reactive cases, we obtain from (3.12) ~ ~ dA dx ~ ::: kHAB ~ ~ ~ (A - (1-r;A)A 1)(B r;Ar;B (1-r;B)BI) (3.15a) (1 -r;B)BI) (3.15b) ~ dB ::: kHAB dx r;Ar;B (A ~ ~ (1-r;A)A 1)(B and and We note that (3.15) is similar in form to (2.13). (2.13) are defined in different vector Although (3.15) and spaces, the physical interpre- tations of the two equations are basically the same. Equations (3.15a) and (3.15b) constitute the basic equations ~ ~ for A and Bas a function of x in a plume in which A and Bare undergoing a second-order reaction. The equations apply to the situation in which the reaction occurs among emitted species in the plume as well as between a species emitted and one in the ambient atmosphere. Subse- quently we apply (3.15) to the computation of the conversion of NO in a plume through reaction with o3 in the ambient air. To employ (3.15) we must obtain inert concentration statistics. The next section is therefore devoted to the problem of simulating the dispersion of a plume of inert material in the atmospheric boundary layer. -84- III.2 Simulation of a Plume in the Atmospheric Boundary Layer The diffusion zone model requires that one have inert concentration statistics for the turbulent field in question. In principle, such data can be obtained experimentally in the field or in the laboratory, or through numerical simulation. Because of the difficulties in obtain- ing field data on atmospheric diffusion, there is growing interest in direct numerical simulation of turbulent flows (Deardorff, l970a,Orszag and Patterson, 1972). These simulations provide a turbulent field 11 within which diffusion experiments may be carried out. 11 11 11 Given the tur- bulent velocity field,diffusion studies can be performed in either a Lagrangian or Eulerian fashion. In the Lagrangian experiments, particles are released from a designated source, and their movements are traced. This technique has been used to study Lagrangian velocity statistics (Deardorff and Peskin, 1970; Riley and Patterson, 1974) as well as the mean concentration from a point source (Lamb et al., 1975). No studies currently exist in which the Eulerian continuity equation has been solved with the fluid velocities derived from a numerical turbulence simulation. In the present work, the Lagrangian particle approach was used in conjunction with the numerical turbulence simulation of Deardorff (1974ab) to generate the statistics required for use of the diffusion-zone model. III.2.1 Description of the Velocity Field The velocity field chosen for this study is that generated by Deardorff (1974ab) simulating the wind field on a particular day of the Wangara study conducted in Southeastern Australia. The simulated -85- turbulence statistics and the temperature profiles have been shown to be in good agreement with the observations of Clarke et al (1971). Briefly, Deardorff's simulation involves the numerical solution of the three-dimensional equations of motion, mass continuity, heat and moisture balances. The calculations use a grid size that falls within the inertial subrange of the turbulence. The grid-volume-averaged equations are solved using a second-order closure scheme (Deardorff, 1973). Initial conditions for horizontally averaged values of velocities temperature and moisture were taken from DAY 33 of Clarke et al (1971) starting at 9 A.M. This particular day of the Wangara experimental data was chosen because of clear skies, very little horizontal advection of heat or moisture, and lack of any frontal activity within 1000 km. To start up the motions on a resolvable scale, a temperature fluctuation was introduced initially at levels below z = 200 m. The simulation covers an area of 5x5 km 2 with a vertical depth of 2 km. This volume is divided into a 40x40x40 grid, each cell of which measures ~x=~y=125 m and ~z=50 m. Du~ .to the vast amount of data generated, only data for the period from 12:29 to 13:21 P.M. local time were stored and available to us. We found this set of wind data favorable in study- ing our integrated plume model because there is little wind shear in most of the mixed layer and because the turbulence is in a quasi-steady state. The major difficulty in using Deardorff's velocity field to compute second order concentration statistics is the need to parameterize -86- the effects of subgrid-scale (SGS) velocity fluctuations. All velocity fluctuations of a spatial scale smaller than the grid size, namely 125xl25x50 meter are not resolvable in Deardorff's simulation and this 1 ~, leads to problems in estimating the mean square concentrations in plumes released from a chimney stack. To be more specific, the mean square concentration of an inert material in a turbulent field is expressed theoretically as (Lamb , 1972a) p (XI XII t IX I t I XII t II) • ~ ' ~ ' _o • o' _o' o <c 2 (x,t)> 1 11 11 11 11 11 S(x' _o' t o )S(x_o' t o ) dx'dt' _o o dx_o dt o dx'dx (3.16) where pis the probability density function that two particles released, one at:~ at time t~, the other at x~ at time t~, are found at x' and x 11 , respectively,at time t; S(: 0 ,t 0 ) is the source distribution and 6V is the sample volume centered at x. Thus, second-order concentration sta- tistics are determined by the statistics of the displacements of particle pairs. If two particles are released several meters apart, only SGS turbulence can effectively separate them and contribute to two-particle statistics . Therefore, to obtain the mean square concentration of an inert species for the present study of chemical reactions in a plume, it is necessary that a simulation of SGS velocity fluctuations be developed for use with Deardorff's field. -87- 111.2.2 Sub-Grid Scale Velocity Simulation The most direct solution to the problem of subgrid scale velocity simulation is a reduction in the grid size in the direct calculation. However, the enormous amount of computing time required pro- hibits this avenue . One must resort to an alternative simulation which produces particle statistics that are as realistic as possible. The two most important statistical quantities in determining p are the one-particle mean square displacement and the two-particle mean square separation. The former is mainly controlled by energy-containing eddies and hence is determined by the supergrid -scale motions. The l atter , however, is controlled by turbulent eddies with sizes comparable to the distance between the two particles. If the particles are initially close to each other, their separation rate is controlled by SGS motions . As the particles are separated farther and farther, large scale motions become more and more important in influencing their rate of separation . Our criterion for an acceptable SGS velocity simulation is that it gives a mean square separation that agrees with available theory. In the inertial subrange, the rate of separation of two particles can be predicted by similarity theory (Batchelor, 1950) ~ (3 . 17) where <£ 2 (t) > is the mean square separation of two particles at time t , or using the notation of (3 . 16); £0 is the initial -88- separation Ix~ - x~I; £ is the rate of energy dissipation; and c1 , c2 are constants of the order of unity. It has been confirmed by observa- tion in the atmosphere (Gifford, l957ab) that (3.17) applies as long as the root-mean-square separation is much smaller than the characteristic length scale of the energy-containing eddies. In Deardorff's calcula- tion, the grid size was chosen so that it falls within the inertial subrange (Deardorff, 1973). Thus, an acceptable SGS velocity parameteri- zation is one which provides a separation rate consistent with (3.17) kc up to a point where <£2 > 2 is comparable to the grid size. There are two approaches to generate two-particle statistics l) Lagrangian - one assigns random velocities to fluid particles and tracks an ensemble of particles 2) Eulerian - one creates an Eulerian velocity field within which particles are imbedded. Let us consider each approach. Consider first the Lagrangian approach . Let up (t) and up 2(t) 1 be one-dimensional velocities of particles 1 and 2. Then the mean square separation at any time can be written as <£ 2 (t)> = <(x = <( = £2 t I 1 2 [up 1 (t') - up 2(t')]dt' + x10 - x20 ) > O t o 2 - x2 ) > t + ff < upl (t' )u pl (t") > dt'dt" t t OO t t + <up2 (t' )u p2 (t")> dt'dt" - 2 <upl (t' )u p2 (t")> dt O O oo (3.18) ff Jf -89- Obviously, the simulated mean square separation depends entirely on the Lagrangian velocity correlation functions chosen. Unfortunately, except for a few simple cases such as statistically independent particle velocities, a systematic approach to obtain the mean square separation from known correlation functions is lacking. Furthermore, there is no general numerical algorithm from which one can generate Lagrangian random velocities with desired correlations. We therefore turn to the Eulerian approach. Eulerian methods have been used by a number of previous investigators. Patterson and Corrsin (1966) constructed a binary field with zero autocorrelation to study one-dimensional dispersion. Kraich- nan (1970) dealt with two- and three-dimensional homogeneous isotropic velocity field using a random wave simulation (see below). Riley and Corrsin (1971) used a similar method but included a weak mean shear in an otherwise homogeneous flow. Thompson (1971) simulated turbulent dif- fusion by constructing an Eulerian field with given correlation. Lind- berg and Thompson (1974) introduced a method of generating a random field with correct Reynolds' stresses. The random wave simulation used by Kraichnan (1970) and Riley and Corrsin (1971) is considered here. This method does not simulate the wave number interactions that constitute a fundamental feature of turbulent flows but it does possess the means by which the energy spectrum can be arbitrarily chosen to give desired turbulent statistics. Since in this study we are not as concerned with the exact features of -90the turbulence as we are with the Lagrangian statistics, and since the SGS spectrum is well-defined, (see below) the random wave method is chosen to simulate the SGS velocity field. Let u(x) be the three-dimensional SGS Eul erian velocity. Following Kraichnan (1970), we have N u(x) = I (3 . 19) n=l where ~n = (Knl, Kn2 ' Kn 3 ), and the random amplitudes A and Bare calculated by (3.20a) ~( .,. X ~n _ ~n) _ = :n (3 . 20b) The vectors s , ~ are independent random variables. ~n ~n sistent with the cont i nuity equation, (Kraichnan , 1970). Eq . (3.20) is con- The choices of K ~n , s_n and ~ _n are discussed below. In Deardorff s simulation, the energy spectrum has been trun1 cated at the wave number corresponding approximately to the grid size. Since this wave number was chosen to be in the inertial subrange, the SGS spectrum follows the well-known 35 power law, namely ( 3. 21 ) -91- where a is a universal constant. the value 1 .81) (In Deardorff's calculations a has The SGS energy is therefore given by f 00 = aE 2/3 K-5/3 dK (3.22) 1T I',. where t,. = (t,.xt,.yt,. z) 113 . The SGS energy ESGS is one of the variables included in Deardorff 1 s data set. With this information and a relation- ship between i:: and ESGS given by Deardorff and Drake (1975), the amplitudes of the number vectors K_n may be chosen. Thus, let the elements of s_n and ~ _n be generated from a Gaussian independent random number generator with zero mean and standard deviation an (to be determined). The components Kn1• of the vector K_n are picked randomly from a spherical surface of radius Kn= IK ,_,n I, with the restriction that IKn,.J -> 2'.... "" This gives (3 . 23) From (3.19 -20), the first component of u(x) is "1(x) = n~l {(sn2Kne - 'n3Knz)cos(~n·~) + (sn2Kn3 - ' n3Kn2)sin( ~n·~)} Since sn,., ~ n1. and Kn1• are statistically independent, each with zero mean, we find that -92- Similar expressions can be found for <u2> and <u~>. By the definition of ESGS and by (3 . 23), we see that our simula t ion i mpl i es N I (3 . 24) n=l In order that (3 . 24) be equivalent to (3 . 22) , we must have N N I n=l I i =1 Kn +~ 2 Kn I CX E: 2/3 K -5/3 dK (3.25) 1 K - .;,,/\ K nc n Integrating (3.25), we have (3 . 26) Using relationship (3.26) , random vectors ~n' ~n' ?n can be chosen and the SGS velocity (3.19) can be calculated . In principle , the number of terms Nin (3 . 19) should be as large as possible . However , in practice si nce we are interested only in particle pairs with initial separations £0 of the order of several meters , a cut -off wave number KN~ ; is sufficient . Considering the 0 computing time and storage requirements of this SGS simulation, we have chosen N to be 15 . -93A trial calculation was carried out using an energy dissipation rate of 0. 00195 m2 /sec 3 . Pairs of particles initially located 1 meter apart were advanced by the velocity field (3.19) through the simple integration: x(t + 6t) = x(t) + 6t u(t). The accuracy of the trajectories produced by this expression is contingent upon 6tu being small compared to the grid size (because the latter is a measure of the smallest spatial variations resolvable in Deardorff 1 s velocity field). To meet this requirement, we chose 6t to be l second . The particles are followed until t sl/3 ~ i 213 - 600 when becomes larger 0 than ~2 • Figure 3.2 shows that the simulated mean square separation is in good agreement with Batchelor 1 s theory. The one-dimensional SGS spectrum (x-component) given by (Hinze, 1959) 0 and the spectrum of the super-grid scale velocity field computed by Deardorff are shown in Figure 3.3 . It can be seen that the SGS simula- tion complements the energy spectrum of Deardorff s velocity field. 1 therefore conclude that (3.19) is an adequate description of the SGS velocity field . We -94105 = - - - - - r - - - - - ~ - - - - - - - r - - - - - - , 103 10 1 + MCJDEL - BATCHELCJR(1950) Cl= 0.70 C2= 0 .18 10- 1 ..._____._._..,&.J..l......._____,....................... 10- 1 J..UJ_ __,___._.-L.L.UJ.L-..___.-J-J...LU.U 102 103 Figure 3.2 Comparison of two-particle mean square separation in a simulated velocity field with Batchelor's prediction. -95- 102 SUPER-6-RlD SUB -GR1D 10 1 "\ \ \ -5/3 \;~· \ 10- 1 \ \ \ \ \ 10 -3 I ..._._...................____.___.__.L-L-L..LJ..J..J.--"--L-L.L..l-LJI.J..L__J.._J__.L.LJ_J...1UJ' 1 IL__...L_.L1..L.11 10° Figure 3.3 One-dimensional energy spectra for supergrid and sub-grid scale velocity fields. The grid scale is denoted by 6. -96- III.2.3 Diffusion Experiment in Simulated Atmospheric Turbulence Having parameterized the SGS velocity field, we can now con- struct an ensemble of particle trajectories, from which Lagrangian statistics can be computed, by releasing particles into the flow field and tracing their motions. A particle located at x has the instantaneous /V velocity, where ~o(:) is the supergrid scale velocity calculated at~ through linear interpolation of the velocities at neighboring grid-point values and uSG is calculated from (3.19). time of the simulated Wangara data. Experiments start at 12 :31 local Sixteen hundred pairs of particles, each initially separated by a distance of 4 meters in they-direction, were released from the center of each grid square on the horizontal plane at z = 125 m. The particles were tracked until 12:55 resulting in a total simulation time of 1432 seconds. In Deardorff's original calculation, the governing equations were solved using a time step of 4 seconds, but velocities were recorded only at alternate time steps . step 6t After some tests, it was found that a time of 8 seconds was sufficiently small, compared to the time scale of variation in the super-grid scale velocity field, that interpolation of the velocities between the recorded time steps was unnecessary. As mentioned earlier, the SGS field requires a much smaller time step . As a compromise to the computing speed we used 68t for the first 1000 seconds and 6 for the remainder of the simulation . The trajectory is then f computed as follows: -97t+llt X ~P (t+llt) = ~p(t) + ~D(~p(t))llt + f t The trajectories of the 3200 particles were stored on tape for use in subsequent calculations of the concentration statistics. Some typical particle trajectories are shown in Figure 3.4. Several checks were performed to test the qualitative features of the particle displacements. First, the mean rise and the mean vertical spread of the particles were calculated for comparfson with other data. (Figures 3.5-6) The non-dimensional quantities are defined as fo 11 ows: (Deardorff, 1972 ) t* = tw*/zi <Z 0 z >* = <Z >/Z. p p 1 * where zi is the mixed layer depth, zp is the z-component of particle displacement, and<> denotes average over the 1600 pairs. The convec- tive velocity scale w* is defined by w? S z.] 1 1/3 (3.27) where g is the gravitational acceleration, em is the mean potential temperature in the mixed layer, and w?s is the kinematic heat flux close to the ground. During the simulation period (local time 12:31 to 12:55). We find by interpolating the published values of zi and w* (Deardorff, 1974b) -98z. = 0.039t + 1134.4 (3.28a) 1 (3.28b) where tis in seconds, zi in meters and w* is in units of m/sec. Deardorff and Willis (1974) observed both numerically and experimentally that <z >* and 0* vary with t* in a manner resembling that p z of a damped oscillator. The oscillation is attributed to the dominance of large vertical velocity eddies whose depths are comparable to that of the mixed layer. The numerical result of Deardorff and Willis (1974) is plotted with our data. (Figure 3.4-5) The agreement of the two appears to be good . The difference in the magnitude of overshooting is probably due both to differences in the particle release heights and differences in the aspect ratios of the two experiments. (Aspect ratio is defined the ratio of the horizontal extent of the simulated or measured fluid to the vertical extent. See Willis and Deardorff, 1975.) As a second check of our particle trajectory calculations, we computed the single particle mean square displacement at and the two particle mean square separation <l 2 >-lt . (Figure 3.7) Two necessary conditions are (1) 2ot ~ <1 2 > - l~ , and (2)<l 2 >- 1o2 should be prop- ortional to t 3 when <l 2 >½reaches a value comparable to the grid size~It is shown in Figure 3.7 that both conditions are satisfied. -99- . >- ~ 1. / . o. o. 1111111l1111.LL f. 2. 3. 4. s. X, KM Figure 3.4 Typical particle trajectories as observed in the numerical diffusion experiment. Simulation period: 1000 sec. Arrow indicates starting point. Particles leaving the field are inserted back at the opposite boundary. - 100- . (D 0 . ::r 0 * I\ a.. N . OEARD~RFF $ WILLISC1974) - THIS STUDY V C\I 0 0 + . o--------'----~~---_..______,__________, 0. 1. 2. 3. 4. 5. Figure 3.5 Comparison of simulated mean rise of particles with that of Deardorff and Willis(1974). -101- . (D a + . + +++++++ ++ +++ + ----+++ ++++ ++++++++++++++++++ ::::r a *bN N a . DERRD~RFF $ W1LLISC1974) - THIS STUDY + a.,_._ _ _____._ _ _ ___.__ _ _ __.___ _ ___.__ _ _ 0. 1. 2. 3. lL ~ 5. Figure 3.6 Comparison of mean vertical spread of particles with that of Deardorff and \·J illis(l974) . -102- 10 6 r - - - - - - - - r - - - - - . . . . - - - - - - - , , - - - - N ~ 10"' b.2 ... ('Jo C( 0 10 3 NO QI( I A ~ V 10 2 10 1 10° .____ 10° -L-,L..L.£..JU.U..._._.,......L...L~U-_._....L..JLI....1...I..UL-.J.........J:.-J...JL...LM.&.I 10 1 10 2 10 3 1011 T, SEC Figure 3.7 One- particle mean square disp l acement and two-partic l e mean square separa tion as observed in the rn1111<'ric,1 l diffusion <'XJJ('Y' illl<'nl. -103Tests have also been performed by setting the SGS particle velocity to zero and observing the changes that result in the mean square separation. It was found that super-grid scale motions begin to domi- nate at about 300 seconds after release and contribute more than 98% of the separation fort> 700 seconds. This result and the calcula- tions shown in Figure 3.6 indicate that the t 3 - region is created partly by the SGS simulation and partly by the supergrid scale velocities. Some Eulerian and Lagrangian super-grid scale velocity correlation coefficients were also computed and are shown in Figures 3.8-11 for future reference. As expected the Lagrangian autocorrelation co- efficient is in general larger than the Eulerian one-point correlation coefficient. Having established confidence in the accuracy of the simulated particle trajectories, we turn next to the estimation of concentration statistics in a point source plume. III.3 Concentration Statistics in a Simulated Plume In applying working equations (3.l5a,b) of the DZ model, we need information on the inert concentration statistics. ~ ~ ~ cally, we need to know A1, s 1 , rAB' sA and ss• More specifi- Throughout the rest of this chapter, we shall concentrate on the problem of a plume containing species A emitted into a uniform environment of B. In this case, we ~ ~ shall see in Section III.4 that rAB = sA = l and ss = 1/rA, where ,_ A2 I -104- 0 • < u.;c o, o> u;., x, o> > < u.'. 2 ( 0 O ) > I { ui = ti •II IAZ, : 1J' I U3::: W' a: C) • OL...-_ ___.:::...::::::::::i:,:_~::::...,c---~....L---~~::S~-~-=:::t:::...,.,=,,,_...-_J o.o -.:.=~ Figure 3.8 Note: 2.5 Spatial velocity correlations for Deardorff's flow field at 12:31 local time (t=O in the simulation.) In Figures 3.8-11, (u' ,v' ,w') and (up,vP,,~p) are Eulerian and Lagrang ian super-grid scale velocities respectively. -105- 0 ......• .( u' <o. o> t.<' co, t n co 0 • R~Ct) = <t.<.'1co,o))~<U'2<0,-t)f.l 'RL (t) = >< Up ct>> , en Zco E) ....... a • , ,, t l/1 ,% l- a: _j W::::t' a: . a:O ID u N 0 0 • • OL------J.---+--....._------'------..1.----__J o.o 0.2 0.6 0.8 1.0 SEC Figure 3.9 Eulerian and Lagrangian velocity correlations of Deardorff's flow field as observed in the diffusion experiment: u'-component. -106- 0 ......• (X) < v'o,o, v'co,t) > 0 <11 •2 lo, o>>r/:z.. <1)'•l co. t n•f1 • ( 'tip' CO > Vp' ( t) '> er, Zco ~ ..,_.. 0 . 1 12 < tfp :i(ol) (V,,'\{J'}h 1- (I -' w:::, a: . a::o 8 u N 0 ·• O L - - - - - - - ' - - - - - - ~ - - - - - - ' - ~ - _ . . . : : : : : . . . . _ - - - - 1_ _ _ ___J o.o 0.2 0.ij 0.6 T, 10 Figure 3. 10 3 0.8 1.0 SEC Eulerian and Lagrangian velocity correlations of Deardorff s flow field as ob served in the diffusion experiment: v -component. 1 1 -107- a • co • a (t..!'/ ( 'l) WplO)) (f) 1/ :(. IA ~ N"r 101) Zco E:) • .__a 1'2 ,,,_ .(W'r ct>> 1(I _j w ::!' a: . a:o E) u N a 0 • . o,____~;:-::::,---~--~--:::::,,--=----::;;,--..----~r--=--=::'.----' o.o -~ Fi gure 3.11 .o Eulerian and Lagrangian velocity correlations of Deardorff s flow field as observed in the diffusion experiment: ~ - com~onent. 1 1 -108- Therefore the minimum information required to apply the DZ model will ~ be A1 and Af· By definition, ~ (x) = -1 A I Ac JL y fz.1 ( 3. 2 ga) -Ly o A2I(x) = _Ale JLYJzi <Ay(x,y ,z)> dzdy -Ly o (3.29b) ~ When x-dispersion is ignored and Lyis sufficiently large, A1 is inde- --- is not. pendent of x but Ay To evaluate Af(x) one must have the distribution of <Ay(x,y,z)>. The general expression for <Af> is given by (3 . 16) in terms of a twoparticle probability density function (pdf). If the sampling volume 6V is sufficiently small * that spatial variations of the pdf are negligible within 6V, (3.16) reduces to the form 000000 (3.30) *The sample volume 6v can not be arbitrarily small or sample fluctuations become seriously large (Lamb, 1972a.) It will be shown later that 6V must be much smaller than <£ 2>3/2. -109- where ~o - = (x o ,yo ,z 0 ), dx~o = dx ody odz o . Assuming the source is a circular area perpendicular to the x-axis, with area as and uniform strength QA,we then have <Ay(~,t)> = Q2A JtJtJ J p(x,x,tjy t ;y' ,t') dy dy'dt dt' ~~ ~o o ~o o where io = (y o ,z o ) and dy~o = dy 0 dz 0 . ~o ~o o o (3 . 31) In stationary turbulence, pis a function of the time differences t - t 0 and t - t~ only, i . e . (3 . 32) where 1 = t - t o and 1 1 = t - t'o· If the plume is slender, x-dispersion can be neglected and (3.32) can be approximated by (3.33) Substituting (3.33) into (3.31), changing variables, and integrating, we obtain <Af(:) > = o: f f P(i,~,tlio•i~) diodi~ (3 . 34) as as where t = x/. Substituting (3.34) into (3.29b\ we then find that for a slender plume emitted from an area source with uniform strength into a stationary turbulence f J P(t,io'i~) d~~dio as as (3 . 35) -110- where P(t ,~0 ,~~) = f p(~ ,~, tl~ 0 ,~~) (3 . 36) dy Ac Judging from the mean square separat ion (3 .1 7), we note that th e dependence of Pony and y0' vanishes as t c113 ; i 213 becomes large . There~o ~ o fore, in the limit of large t , we can write (3 . 35) in the form q2a 2 A S 2 A1(t) ~ zA P(t, ~, E: ') C ~ (3 . 37) where ~ and ~• are any two points in the source area a . When tis small , ~ s ~ in particular t c 113 i~ 213 .2. 5, Pis not independent of ~o and~~ and the general expression (3.35) must be evaluated numerically . This poses quite a computational burden, however, because the evaluation of P for a single set of initial release points (y 0 ,y~) requires a large amount of computer time. Therefore, we seek a way of reducing (3.35) to a form like (3.37). This can be achieved if we can find the set of vectors (y. . . o ,y') for which ,.,,o f f P(t ,~o'~~) d~~d~o (3.38) as as because with this expression we can write (3 . 35) in the desired form (3 . 39) Let us attempt to estimate (~ 0 ,~~) • -lll- III.3.1 Simplification of Equation (3 .35). For an arbitrary function ~(x,x 0 ,x~), we note that co VV where t co = x - x~, and g(t ,x;) is the joint pdf that any two points chosen at random in the space V will be separated by t point will be located at x~. and that one Applying this relationship to the RHS of (3.38) and letting t 0 =yo- yl , we have ( 3. 40) co co To proceed we need a description of P but unfortunately little is known about the form of the two-particle density p(y,y 1 ,t jy ,y~) of which P ~ ~ is a function (see (3.36)). _o ~ It is known only that in homogeneous sta- tionary turbulence the marginal density p(~•!ly 0 ) is Gaussian (Batchelor, 1949). Therefore, following the study of Lamb (1972a~,we shall assume that p(y,y 1 ,tjy ,y~) is jointly normal. (It turns out that the estimate - ~ ~o ~ of (~ 0 ,i~) is only a weak function of the form assumed for p. ) We then have for the case of isotropic turbulence ( 3. 41) -112- p = <(y-<y>)(y•-<y•>)>/a2 = <(z-<z>)(z 1 -<z 1 >)>/a 2 It is easy to show that p = l - <t2>-t2 0 4cr 2 (3.42) where t2 = (y _ y•)2+ (z _ z')2 o O O O 0 With this assumption, P(t) is obtained by setting y' = y and integrating y over Ac. By virtue of the definition of Ac, integration over this domain is equivalent to an integration over an infinite region. l (2/ncr) 2 (l-p) Thus we have (<y>-<y' >)2+ (<z>- <z' >) 2 ] 4cr 2 (l- p) (3.43) In the diffusion experiment, there is no mean motion in the y-direction. Hence ' <y> -- y0 and <y' > = y0' . However, the mean verticlc position of each particle does change with time (see Figure 3.5). This mean motion is a weak function of height of release; but since the area of the source is small, the difference in mean rise due to different release heights of particles is negligible. z0 - z~, and (3.43) becomes We thus have <z> - <2 1 > = -113- (3 . 44) where !l = !l ~o o I . Substituting (3 .44) into (3 .40) we get P(t' ;:;,., 0 ) = J 1 ~ 4cr 2Jl~(l- p) l.J g( fl o ) dil o 4,ro 2 (l -p ) exp [ (3 .45) 00 where 00 00 The remaining task is to evaluate the right hand side of (3 .45) and to find an initial separation I 0 such that (3.45) holds. The function g( il ) can be obtained numerically by selecting 0 points at random from a circle of radius rand counting the frequency of occurrence of each value of £o . The function rg( il 0 /r) has a universal form for any value of rand is shown in Figure 3. 12. Note that it has a maximum at £0 /r = 0.66. The parameters cr and p required in (3.45) are taken from the diffusion experiment data. Fitting a polynomial to each of the curves in Figure 3. 7, and applying (3 .42), we obtain 02 =} cr f = 0. 646t 2 + 5.513t (3.46) (3.47) -114- 10 1 7: = 50 sec --·--·--1/) en (lJ ~ 0 ..... 10-1 T/2 150 sec --•-. -·---·--- C QI 250 sec .....E -·--· -·-· p .. -- -·--. ------ 10-2 0.66 rg (Qo/1) Legend: as Pct, Qo) ~ r= 5.56 m 0 0 .Qo(t) Joi r = 0.66 Figure 3.12 Calculated distribution functions g(l0 /r) and P(t,l 0 ). The location of I'0 of the mean value of Pis strongly associated with 9max· -115These relationships and the function g allow us to determine I0 as a function of travel time for various values of r. Figure 3.12shows a family of curves of aS P(t,£) for the case of r = 5.56 m. The locus 0 of aSP(t,i) is indicated by the circled curve in Figure 3.12. It is 0 seen that£ 0 moves from 0.59r at t = 50 sec. to a limiting value 0.66r as t gets larger. Clearly the location of i 0 is strongly associated with the behavior of the weighting function g{ i 0 ). has a maximum at£ 0 = 0.66r.) (Recall that g(t) 0 These results suggest that except for small values oft, the functinn P can be approximated from an ensemble of trajectories of particle pairs of initial separation of i 0 = 0.66r where r is the radius of the source area. - Once P(t,I) is obtained, 0 Af can be computed from (3.39) with P(t,~ 0 ,~~) replaced by P(t,i 0 ) . - In summary, we propose the following approximation for Ay(t): I Q2a2 A2(t) = As p2{~,t) dy I 2A A C (3.48) C where p2(~,t) = p(~,~,tlt 0 ) is the pdf that two particles released a distance i 0 apart will be found together in the sample volume centered at y at time t. The task now is to measure p2({,t). III.3.2 Sampling Procedure and Statistical Error of the ,,_ Estimate of Af· ~ As noted earlier, the sample volume ~v used to estimate Af must be sufficiently small that spatial variations of p2 within ~v are negligible. To estimate the size of ~v, we deduce from the joint-normal -116k form (3.41) that the length-scale of variation of p2 is cr(l-p) 2 • By (3.42), we have Since p varies from Oto l, this length-scale is approximately 11y = 0. l £ y( t) , ( 3. 49) 11z=O.l £z(t) where £y (t), t k 2 2 z(t) are they- and z-components of <£ > at time t. (The adequacy of this sample volume size was confirmed through a series of test runs.) Using this sample volume, we computed the mean square concentration in the following manner; First, the coordinate system of the recorded particle trajectories was rotated to align the x-axis with the mean wind direction. The mean wind speed is 2.74 m/sec at an angle of 156.88° with respect to the east-west axis. Next, the coordinates of each pair of particles were translated to a new system (n,~): -117- Each pair of particles was then considered to constitute one realization of the ensemble of particle pair trajectories. After the above opera- tions were performed, the sample volume ~v was determined from (3.49) and the averaging area Ac, defined by (3.3), was divided into cells of size ~V. Within this framework p2 may be approximated by (3.50) where NR is the number of realizations and l t,,VZ ' 0 if both particles of the i-th realization are located in t,,1/ centered at y at time t; otherwise. On substituting (3.50) into (3.39) we obtain (3.51) ,- That (3.51) provides an unbiased estimate of Aj(t) can be seen by taking the ensemble average of this expression and using the definition p2 = <ci>. We find (3.52) -118,...., This result insures us that we can make the error in the estimate E{Ay} ,-., of Ay arbitrarily small by making NR sufficiently large. However, for all practical values of NR, there is an error associated with E{Af}. Let the mean square error in (3.51) be defined by ,..,..,, X = <(E{Af} Applying (3.53) and (3.44), we obtain 1 w ~R l R i=l ~R ~2 1 l <c.(y)c.(y )> dydy 1 - Aj 1 j=l ~ J ~ (3.53) In our experiment, each particle pair released from every grid is considered to be one realization, but these realizations are not statistically independent because of the proximity of the particle pairs to one another at their time of release. Suppose that for each i, ci is perfectly correlated with mother realizations but is statistically independent of all the rest, i.e. <c. (y)c. (y ) > 1 l ~ J ~ if the i-th and the j-th r ealizations are correlated; (3.54) otherwise. Substituting (3.54) into (3.53), we find that (3.55) -119- Recalling the definition of ci(y), this function can have non-zero value at one location, namely if y' = y ; if y' 1 y. Applying this relationship to (3. 55) we have 2 2 2 _ m+l QAaS > dy_ - A2I 2] -N [ ( 2A) tiV J <c?(y) 1 _ R C Ac X - The mean fractional error£ is then tiV where = y I <cf (r) >dt Ac ( 3. 57) To estimate£, we must first approximate y. To this end we assume that in each realization the particles are equally likely to be found anywhere within a box of size ntiV. Two extreme cases are considered: 1) If the two particles are totally correlated, they are always observed together and consequently -120- <c. (y )> = l ~ l -n { { (k,)2 otherwise, 0 ' l ) nl ( w <c?(y)> = l ~ if yE ntN; 4 , if y -E ntN, otherwise 0 ' It follows that in this case y = l and£= 0. 2) If positions of the two particles are statistically independent, then 2 l <c.(y)> = l ~ { n2"N if Y-E ntN, 0 otherwise l <c?(y) > = l ~ { (1 ) ( l) ~ ~v 4 if y -E ntN otherwise and it follows that y = n and We note that the error is larger in the statistically independent case. This is because when the particles are correlated, they are more likely to be found together and fewer realizations are needed to detect all possible positions of the particle pair. For the same number of -121realizations, the error should be smaller in the correlated case. In reality, particle pairs are correlated initially and then become independent gradually. Hence the true error is expected to lie somewhere between the two error bounds given above. Let us assume that (3.58) Based on the spatial velocity correlations in Figures 3.8 - 11 , we estimate that initially m has a value of 10 and decays to zero in about 500 seconds. Using NR = 1600, ntN = 4oyoy and taking p from (3.42), Eq. {3 . 58) yields the following values of E : t, sec lO 100 200 400 800 1000 1400 % 1.1 2.3 2.6 3.2 3.4 5.0 6.4 € ' 0 These results indicate that the statistical error is negligible during the period of simulation. III.3 . 3 Numerical Results Let us assume that on N of the NR realizations the function ci is non-zero somewhere in Ac. (Recall that by definition ci can at most have non-zero value at one location.) ,.,J _ Af (t) - Eq. (3.51) then reduces to Q2Aas2 rN N l 2 AN tN 2 C R ( 3. 59) ~ Once N has been determined from the numerical experiments, Ay(t) can be calculated from (3.59) and the result combined with -122- (3. 60) to obtain the desired quantity (3,61) The computed values of rA normalized by its initial value, i.e. is plotted in Figure 3. 13 (Note that rAo = Ac/as, which follows from (3.60) and the fact that at t = O,Af = Qjas/ 2 Ac . ) The figure shows that the measured quantity does approach the correct analytical limits ~ ~ ~ ~ ~ both at t = O (rA/rAo = 1) and at t ➔ (rA ➔ 1). As expected, rA/rAo 00 decays continually with travel time because the plume is gradually mixing with its environment. ~ ~ We note that rA/rAo has decayed to one- tenth of its initial value after about five minutes. Deardorff (1972) found that mean properties of the unstable planetary boundary layer can be made universal if all lengths are scaled by zi' the height of the inversion below which convection is confined, and times are scaled by zi/w*, where w* is the convective velocity scale defined in ( 3.27). The travel time t* nondimensionalized in this manner is given in Fig~.13. It is found through a least squares fit that -123- NUMERICRL DRTR - CURVE FITTED T~ DRTR + ,..__, r;, r / AO 10-2 10-4 .___ _ __.___ _ __.__ _ ___.__ _ ___.'--_ _____, 0.0 0.5 1 .0 1 .5 2.0 2.5 0 5 10 15 20 25 T, MINUTES Figure 3.13 Second- order concentration statistics of an inert species in a simulated plume. Release height: 125 m. Source area: circular, radius 5.56 m. Inversion height Zi:1100 m. Convective velocity scale: 1.96 m/sec. Stability: -Z/L~911, where L is the Monin-Obukhov Length. -124- ( 3. 63) where a1 =-1.736, a 2 = - 0.111, a3 = 0.359, a4 = - 0.095. In the next section we utilize these results in the diffusion zone model to simulate a reactive plume in the planetary boundary layer. III.4 Study of the N0-0 3 Reaction in Simulated Atmospheric Turbulence As we noted previously, the two major air pollutant sources for which turbulent chemical reactions may be important are highways and power plants. A pollutant emitted in significant magnitude from each source is NO. Since NO, once emitted, then participates in a com- plex series of reactions leading to photochemical smog, understanding the fate of NO from these two types of sources is important in understanding photochemical smog chemistry. Hilst et al. (1973) studied a case of NO emitted from a highway source. We concentrate here on NO emissions from an elevated point source in a simulated unstable planetary boundary layer containing a uniform concentration of ozone. The specifications of the problem that we shall treat were formulated during an early phase of this project when observational data were not available to us. Just recently, however, we discovered measurements by Davis et al. (1974) of NO, o3 , and so 2 concentrations in plumes produced by the power plant in Morgantown, Maryland. It turns out that the conditions under which these data were measured are sufficiently similar to those of the problem that we simulated that these data provide -125a partial check of the validity of the diffusion zone model and our implementation of it. After describing the numeri cal simulation, we pre- sent a comparison of the predicted concentrations with these measured by Davis, et al. It has been observed, e.g. by Davis et al. (1974), that in the initial position of power plant plumes, the reactions that dominate ozone concentrations are ( 3. 64) (3.65) The value of k1 at 25°c is 29 ppm-l min-l (Johnston and Crosby,1954). A typical value for k2 is 0.3 min-l, depending on the intensity of sunlight . Let A= [NO], B =[0 3], and C = [N0 2], the local instantaneous rate of reaction is In this case, the diffusion-zone model for the extent of reaction in the plume takes the form ~ ~ dA _dB_ dt - dt - - ~ ~ - ~ ~ [A-(1-sA)AI][B-(l-ss)BI] +k2C (3.66) -127- weakly sensitive to the value ofµ a . We therefore assumeµ a = 1 as an approximation. Figure 2.2 indicates that this assumption is more criti- cal at small than at large t. Hence, initially the predicted reaction rate may be underestimated slightly. The value of the reaction parameter A was determined from the formula developed in Chapter II (Equation 2.41 ). According to (2.42), proper value for the parameter a that enters in (2.41) is based on the minimum of the initial concentrations and the maximum of the diffusivities. Since NO is always released in a higher concentration than the background o3 concentration, we have 1 1+OJ 6a A=-,,--~ (3.69) The background concentration of o3 assumed in this simulation (and typical of the observations of Davis et al.) is 80:!:_3 ppb. 0.003 and A; 1. ~ This value gives a = Reactions with a values of the order of 0.01 to 1 are considered to be of intermediate speed and their rates are influenced by the rate of mixing. Applying (3.67) and the condition of A= 1 to (3.66), we find -dA = - dt ~ k r~ A B-(1 - -)B 1 ~J + k C 1 A~ I 2~ r (3.70) A Since (3.64-65) are the only reactions among NO, N0 2 and o3 , stoichiometry requires that -128- substituting (3 . 71) into (3.70), normali zing by the constant A1, we have where FA= A/A 1. Applying the relationship (3.60) and the fact that rAo = Ac/as, (3 . 72) becomes _ k ( ~\ l u ~ A ) FA [FA _ l + Bo ( Q r~ Ao ) r Ao AA l (3 . 73) Equation (3 . 73) describes the decay of total amount of NO emitted due to chemical reaction with o3 . The rate of decay is a function of the kinetic rate constant, background ozone level and turbulent mixing. Con- sistent with the observations of Davis et al . (1974) , this model pre dicts a higher conversion with a higher background ozone level (larger B0 ) . In order to eventually compare the results of this simulation with the plume data of Davis et al . , we need to know the emission flux QA of NO in the Morgantown plume. Unfortunately, the value of QA is not available because the reported data were normalized by the local -129NOx concentrations * . Therefore some estimate has to be made. In Appendix C we estimate the NO-emission of the Morgantown power plant to be between 6500 and 11000 lb/hr at the stacks. Because of the tre- mendous amount of thermal energy exhausted by the power plant, the plume near the stacks behaves like a thermal jet rather than a passive plume carried by_a turbulent fluid. Thus, only the portion of the plume dominated by ambient turbulence can be compared with our simulation. For this reason we treat the data as representing a plume emanating from a fictitious source perpendicular to the plume axis and located at the point downwind where the plume first becomes horizontal. (see Figure 3.14) Suppose that at this point, the plume radius is four times the radius of the stack. Assuming the twin stacks have radii of 5 m each, the average emission flux QA at the fictitious source is estimated to be 2.5-4.4 lb/hr.~. Using a mean wind speed of 24 km/hr as observed by Davis et al., we find the fictitious source concentration QA/ to be 39-67 ppm (assuming that each cubic meter of air is 0.089 lb-mole). Estimation of the location of this fictitious source so defined requires more information such as plume temperature etc. than we have. However, Davis and Klauber (1975) reported that at 0.8 km downwind, NO concentration was~ 2 ppm. This suggests that the fictitious source must be rather near to the stack, and so we shall *One of the authors indicated to us that the normalization NOx was not recorded. Figure 3.14 ~ Source Fictious X NO 03 Schematic of a power plant plume. The fictitious source is located at the point downwind where the plume first becomes horizontal. z t I 0 w I-' -131- assume that no appreciable chemical conversion has taken place before NO reaches the fictitious source location . Admi t tedly , these estimates are rather crude but without them we have no basis for judging the accuracy of our simulation. Using these estimates and (3.63) , we solved (3.73) for FA with the initial condition FA(O) = l . The solution FA is plotted in Figure 3. 15 as a function of in (3.28). t* = tw*/zi, where w* and zi are given The data of Davis et al. (1974) are shown on the same graph for comparison . To connect these data with the t* coordinate, we first transformed the distance-scale to real time-scale by using a uniform wind speed of 24 km/hr . Then we normalized t using z.l = 1000 m, supplied by Davis (private communication), and w* = 1.8 m/sec est imated from the observed cloud cover, wind speed , and approximate surface roughness of the Morgantown area. The error bar on each of the data points indicates the range of observed values . Due to the uncertainty in the estimate of NO-emission rate, we computed FA for both the maximum and the minimum of the estimated QA/ values . Although there are too few data points to make a definite judgement, the computed and observed conversion rates are in good agree ment. The fact that the observed conversion is greater at small travel times than that predicted by the model suggests that the buoyancy generated turbulence in the initial portion of the plume is more effective in mixing the NO and ambient ozone than we had assumed . Recall that since buoyancy effects are oot included in our simulation, we treated -132- D . ..--, r---ac:==:::::--- co • D . (D D :::f' D • DRVISETRLC197L!) 0 . Z. MClDEL c:i . ('\J D 2 0 QA/ = 67 PPM QA/ = 39 PPM . o.____ _ __,____ _ _- - - - ' - - - - - - L - - - - l - - - - - - - ' 0 .0 0.5 1.0 1.5 2.0 2.5 15 20 25 T* 0 5 10 T, MINUTES B0 = 0 .08 PPM = 2.71..! M/SEC Fi qure 3. 15 Reactive decay of the total mass of NO emitted fro m a power pl ant stack into an ozone laden atmosphere. So lid curves are the DZ model prediction, wh il e circl es are obs ervations. -133the observational data as though the observed plume were neutrally buoyant and issuing from a fictitious source some distance from the stack. Part of the underestimate of the computed conversion may also be due to our assumptions that µa in (3.68) has unit value. There may also be some effect of the differences in the release heights of the simulated and observed plumes . An important aspect of the computed and observed data that are in close agreement is the magnitude of the time scale of conversion of NO. As is evident in Figure 3.1~ in both the simulation and the observations 50% of the NO released has been converted to the N0 2 during about the first 17 minutes of travel. Davis et al. (1974) re- ported that the conversion times they observed ranged from about 12 to 60 minutes depending on the meteorology and background ozone concentration. In view of the rather large computational effort that is necessary 'to derive the parameters (such as rA) required to implement the diffusion zone model, it is appropriate to inquire whether a plume simulation of comparable accuracy could be achieved using a simple, conventional approach. To investigate this matter, we constructed a Gaussian type model of the NO plume treated above and excercised it under the same conditions as those implicit in the data shown in Figure 3.15. The model equation is (cf.(3.2)) -134- (3.74) where Ky and K2 are turbulent eddy diffusivities . For si mplicity we assume Ky= K2 =Kand is a function of time on l y. If it is also as- sumed that the plume is axisymmetric , the n (3. 74) becomes (3 . 75) where t = x/ and r = ✓y 2 + z 2 . i~e assume that K(t) is related t o the mean square displacement of a single particle by K(t) = ~f (3 . 76) and use 0 2 (t) as given in (3.46). These assumptions are not critical, because since the conversion FA is defined by rX)2nr<A(r,t)>dr 0 (3 . 77 ) the effect of cross -sectional diffusion is reduced . Test runs indi cate a wea k dependence of FA on the eddy diffusivity K(t). The nonlinear kinetic term in (3.75) can be decomposed into two ~arts : <AB>= <A> + <A 1 B1 > (3 . 78 ) -135The fluctuation term <A'B> represents the degree of segregation of the reactants, and for fast reactions it limits the rate of reaction . It was for the purpose of describing concentration fluctuation effects that the diffusion-zone model was developed . In keepi ng with the pres - ent practice of neglecting fluctuation terms, we shall omit <A'B> from the Gaussian plume model equation . With this assumption and the st oich i- ometry condition + <C>= B0 , we obtain a<A> _ a _ K(t) ( a <A> + l a<A> ) k A k (B ) a""t - a-r- ar 2 r ar - 1< > + 2 o- <8> • 2 (3.79) The initial conditions are if O < r < r - = ci, 0 <A(r,0)> (3 .80a) \ otherwise <B(r,O)> = B 0 for all r (3.80b) and the boundary conditions are a = 0 ar ;J<A >- -ar <A> = 0, where r 0 = Bo at r =O (3.81a) at r ➔ "° is the source radius (=5 . 56 m). (3 .81b) We choose a value of QA/= 47 ppm and solve (3.79) with (3.80-81) numerically. is then used to compute FA by (3.77). The solution <A(r,t)> The results are plotted along with the DZ model prediction in Figure 3.16. -136- 0 .--: r , - - - 1 I I I . co 0 I I I I I I I I . (D 0 I I I I I I . \ :::r 0 \ ORVIS ET RLC1974) - D. Z. M0DEL -- GAUSSIAN M0DEL ~ \ \ \ \ \ ('\J 0 • --------D . o...____ _ _....____ _ ___.__ _ _~ - - - - L - - - _ _ J 0.0 o.s 1.0 1.5 2.0 2.5 15 20 25 T* 0 5 10 T' MINUTES QA/ = 47 PPM Bo = 0.08 PPM = 2.74 M/SEC Fi gu re 3. 16 Comparison of predictions of the DZ model and a Gaussian model. The differe nc e in the predic~ tion s is attr ibuted to the ne1 lect of <A' l3 ' > in the Gaussian 111odel. -137It is immediately clear that the Gaussian model (3.79) is unacceptable: It predicts a 50%-conversion time of NO of three minutes and a much more rapid initial rate of reaction than is observed. If the value of QA/ is increased, the prediction of FA comes into closer agreement with the observations . However, the NO-emission at the Morgan- town plant would have to be 250,000 lb/hr for the predictions of the Gaussian model (3.79) to agree with the observations. Since the depend- ence of FA on K is weak, the difference between the prediction of the DZ model and that of the Gaussian model can only be attributed to the neglect of the fluctuation term <A1 B1 >. When NO is first emitted, the o3 in the plume itself is quickly depleted by the reaction, and the excess NO must wait for turbulent mixing to supply fresh o3 from outside the plume for the reaction to proceed. In this initial period, reactants are basically segregated and <A1 8 1 > has a negative value comparable in magnitude to that of the product <A>. Thus , neglecting <A1 8 1 > causes an overestimate of the reaction rate as is evident in the rapid decay of FA in the Gaussian model. By contrast, the DZ model, taking into consideration the effect of <A1 B1 >, predicts a minimal decay of FA during the initial period. III.5 Discussion In this chapter we have developed a special version of the DZ model for simulating a point source plume in the planetary boundary layer. The parameters required to implement the model were derived -138- from data generated by Deardorff's numerical turbulence model. Because the turbulent structure of the unstable boundary layer possesses certain universal characteristics (Deardorff, 1970b), the parameter values we derived may be used in an extremely wide range of simulation studies provided that variables are expressed in terms of the length and time scales zi and zi/w*, respectively, introduced earlier. The large amount of computing time required to derive these parameters is therefore justified. Within the limitations of the available data, we have shown how the DZ model provides a realistic description of chemical reaction in the planetary boundary layer. Unlike many of the· currently available c1osure schemes, the DZ model does not introduce additional differential equations. Also, the difficulty in determining eddy diffusivities was circumvented in our boundary layer plume simulation by computing the total material in cross-sections of the plume as a function of downwind distance. This approach may prove to be of considerable value in air pollution studies. At the same time it renders the computation on applying the DZ model to be extremely simple. The effect of inhomogeneous mixing on the rates of atmospheric chemical reactions can be studied using the DZ model. We compared the prediction of the DZ model with that of a Gaussian model in which the concentration fluctuations were neglected. The comparison indicates that the rate of the N0-0 3 reaction is greatly reduced by inhomogeneities in the concentration fields produced by atmospheric turbulence. For this reason the common practice in air pollution modeling of neglecting -139- concentration fluctuation effects can produce large errors when applied to reactions of intermediate and very rapid rates. We have therefore shown that the concept of diffusion zones is extremely versatile and useful in studying chemical reactions in turbulent flows . Further development of the DZ model should focus on extensions of the model to treat complex reaction systems and applications to sources of various types . There is also the need to develop expressions for the mixing parameters sA' s8 that are applicable to situations other than those we have considered here. -140- IV. LAGRANGIAN MODELS OF TURBULENT CHEMICAL REACTIONS As noted in Chapter I, there are two general approaches to describe the phenomenon of turbulent chemical reacti on : an Eulerian approach, in which the problem is formulated in a coordinate system fixed in space; and a Lagrangian approach in which the behavior of the individual fluid particles is treated. In the previous chapters we proposed a hypothesis to handle the closure problem that one encounters in the Eulerian approach. In this chapter, we view the phenomenon of turbulent chemical reactions from a different perspective by analyzing it with Lagrangian* methods. Since the Lagrangian description addresses the motions of the individual fluid elements that constitute the fluid flow, it is perhaps more natural physically than the Eulerian description. However, a factor that hinders the wide use of Lagrangian descriptions of turbulence problems is that in order to measure Lagrangian properties, one must somehow track individual fluid particles. This is often difficult to achieve experimentally; and consequently, while the mathematics of Lagrangian methods are inherently simpler than those of Eulerian closure methods, the required multi-particle statistical properties needed for Lagrangian analysis are virtually impossible to determine experimentally . Due to these difficulties very few studies have been based on this approach . (See review in Chapter I . ) *Fora detailed discussion on the Lagrangian description, see Monin and Yaglom (1971), Chap. 5. -141Recently Lamb (1972a) has succeeded in deriving expressions for the concentration statistics of inert foreign particles in a turbu1ent fluid from a Lagrangian property that we sha 11 refer to as the occupancy probability density function (opdf). (See Eq. (3.16).) The opdf contains information on the trajectories of the particle motion and in some cases it can be measured experimentally. For example , a demonstration of the measurement of the two-particle opdf was given in Chapter III . However, in cases where a number of different species of particles interact chemically, it is generally infeasible to determi ne such an opdf for each type of particle. Therefore Lamb (1972b) further derived an approximate relationship that enables one to determine the opdf of a particular type of species which is undergoing a second order reaction with itself from the opdf for the inert particles. Lamb's approach provides in effect both a circumvention of the Eulerian closure problem and a different insight into the nature of nonlinear, turbulent reactions. In particular, this approach can provide concentration statistics of a turbulent plume near its source where many Eulerian closure approximations break down . It is the purpose of this chapter to 1) demonstrate this technique, 2) show that it can be extended to multicomponent chemistry, and 3) derive several simplified models from the general description. These extensions of the approach are necessary in order to facilitate its application to air pollution chemistry. -142IV.l A Lagrangian Description of a Second-Order Chemical Reaction in a Turbulent Fluid The problem under consideration in this section is the follow- ing. Given the initial distribution of, say, N inert particles and the opdf of the trajectories of these particles, what would be the marginal densities of the opdf if the particles were chemically reactive? Let the N particles be initially located (t = 0) at where m y., m z.m) r.m = ( x., 1 1 1 1 is the position vector of the i-th particle at time t = m~t where ~t is some small time interval to be defined subsequently. Let the trajectory vector that the i-th particle has followed in m steps in a discrete time-continuous space random process be m m-1 , ... r.,r. 2 1) R.m = ( r.,r. 1 1 1 1 1 The opdf that gives the probability density of the trajectories of all N particles can then be defined as (4.1) The probability density of the trajectories of all N particles, p(R~, R;, ... , R~!r 1,r 2, ... , rN), is defined such that p(R~,R;, ... , R~I -143- dR~ is the probability that the j_th particle followed a trajectory such that after the first time step the particle occupied the differential volume of space rt to rt+ drt, after the second time step occupied the differen tial volume of space, rt to rt+ drf, etc., form steps . Thus, dR~ denotes dr~ dr~-l . .. drt. The opdf defined in (4.1) is strictly a Lagrangian property of the turbulence because it is defined in terms of the positions of inert particles. It has been normalized such that (4.2) co co where and m m m dr.m = dx.dy.dz . l l l l m is a differential volume element centered at ri. If the opdf in (4 . 1) is known, the probability that particles follow a particular trajectory can be generated from it simply by integrating over the given trajectory. For example, the probability that the i-th particle will follow a paths in m time steps and the other particles will follow arbitrary paths is I Is 00 (4.3) -144where an overbar denotes the complement, i.e. , m m m Rf!ll = Rlm'R2' ••• ' Ri - 1 'Ri +1 ' ••• ' RmN (4 .4) m m m m r~l = r m 1 ,r 2 , ... , ri-l'ri+ 1 , ••• , rN (4 .5 ) It is commonly assumed that chemical substances can react only when they are sufficiently close to each other . Therefore, we ma intain that in a turbulent fluid, particles have a chance to react only when their trajectories pass within some specified distance of each other at any instant . In particular, we propose the following Lagrangian model of chemical reactions. A+ B + C. Consider the bi-molecular chemical reaction Suppose that each particle has a uniform chemical composi- tion, that is, it can be either of species A or B or the product C. We envision that associated with each particle there is a reaction volume defined as where <v' 2 > kc2 is the r.m . s. amplitude of the particle velocity fluctua - tions, and a is a characteristic cross - sectional area of reaction for the particle. kc (If the turbulence is stationary and homogeneous, <v' 2 > 2 is a function of neither time nor position . ) Within this framework we postulate that if the trajectories of two particles are such that their reaction volumes intersect , then reaction occurs at the point of intersection with probability one. Therefore, during any interval ~t, the -145- probability that the i-th particle (say of species A) does not react is the probability that its reaction volume does not intersect that of any particle of B during that interval. Since we are approximating a con- tinuous process by a discrete one, the smallness of tt is necessary to insure that the probability of reaction is not overestimated. For this reason we stipulate that tt must obey the condition, or (4.6) Here <Cmax> can be interpreted as the average volume occupied by each particle in the region of maximum concentration, and TR is a reaction time scale. After two particles collide, it is possible that they form either a single new particle or two new particles. We assume that in the first case , the product species can be represented as two individual uniform particles moving together. In this way the total number of particles under consideration is conserved. In a situation where NA particles of species A and (N-NA) particles of species Bare initially distributed at r 1,... , rN and A rNA+l , ... , rN' respectively, the probability of finding that the i-th particle has identity A and is in the differential volume dr~ centered at r~l at time t = mtt can be written as -146( 4. 7) According to the reaction model outlined above, (4.7) is equivalent to the product of the probabilities of two independent events, namely that: 1. The i-th particle was originally an A. 2. The i-th particle never collided with an 8-particle, i . e. it managed to stay outside the reaction volumes of all 8-particles on its route torim. Let the probability of the first event be pi. Clearly P; = Cl (4.8) 0 Let bt be the region of cartesian space (x,y,z) outside the reaction volumes oV of all B-particles during the £-th time step, that is J (· ) ctr = J (·) dr - ~ J (·) dr bt oo J (4. 9) oV. J where j sums over all existing 8-particles in the £-th step. Then if we restrict attention to passive reactions, i.e. those that do not perturb the fluid motion, then the probability of the second event above is easily established from (4.3) to be -147- ( m m- l , ... , r li , mRi l ) { J fbl fb2 . . . fbm-1 p r i , r i oo , , , , . dr.m- 1 . . . dr.'dRm-l)d . r.m (4 .1 0) where pis the opdf of a collection of inert particles introduced in (4.1). The condition (Jr 1 , ... , rN) is not explicitly stated in (4. 10) but it is understood that all densities under consideration are conditioned on the initial distribution. Let Bt be such that J p dRt = I I ... I p dr.t , , , dr~dr '. (4.11) then from (4.7), (4 .8) and (4.10) , (4.11) we have Hence , , , , , m m-1 ,R. m-1 ) dRm.-1 dRm.-1 p(r.,R. (4.12) This expression gives us the marginal density of finding the i-th par- , ticle at r~ with identity A at the end of them-th step . In order to make (4.12) useful, we have to establish a rela- , tionship between p(r~,A) and the density of previous steps. tions (4.11) and (4.9) we have By defini- -148= p. 1 II II ( 00 8m-2 00 m m-1 ,r.m-1 ,R.m-2 ,R. m-2 )dr.m-1 dr. m-1 dR.m- 2dR m p r.,r. .-2 11 1 1 1 1 1 1 1 bm- 1 m m-1 , r.m-1, R.m-2 , R.m- 2)drm-1 m-2 p( r.,r. . dr m . - 1dR.m-2 dR. 11 1 1 1 1 1 1 1 -1 ,R.m- 2 , R.m- 2) dr m -1dR.m- 2dR.m- 2 p( r.m , r.m-1 ,r.m . -1 dr.m 11 1 1 1 1 1 1 1 (4. 13) where j sums over all existing 8- particles in the (m-1)-th time step. Further reduction of (4.13) requires an expression for the probability that the j - th particle be of species Bat the (m-1)-th step. For this to occur we must have: l. The j-th particle was initially a B. (This event occurs with probability (l- pj).) 2. The j - th particle has a trajectory in Am- 2 that lies out- side the reaction volumes of all A- particles in the time domain O < t < (m- 2) 6t. Therefore we see that the second term in (4.13) is essentially the joint probability density that the i-th particle was an A until the (m-l)th step , that it lies within the reaction volume of the j-th particl e at (m- l) 6t , and that i t is at r~ at time t = m6t . Concurrently the j -t h part i cl e had to be a 8-particle initially , it must not have suffered any collision with an A particle during its travels prio r to time t = (m-2)6t , and it must lie at rj-l at time t = (m-l) 6t. Since any of the particles -149- with subscripts NA+l , .. . , N can be a B-particle , this joint probabil ity density can be defined as (4 . 14) where r m-1 .J. = complement of ( r.m-1 , r.m-1) and Rm-2 . . = complement of 1 J l ( R.m-2 , R.m- 2) . l J lJ Substituting (4.14 ) into ( 4.13 ) , we have p(rm . ,r.m-1 ,r.m- 1 ,m r ..-1 ,R.m- 2 ,R.m-2 ,Rm-2) .. drm.-ldrm.- ldrm_-. ldRm.-2dRm.-2dRm.-.2 l l J lJ l J lJ l J lJ l J lJ (4 . 15) In essence all we have done is decompose (4.12) into two parts : one associated with the event that the i-th particle suffered no collis i on with Bin the first (m-2) steps (but its path in the (m-1)-th steps is arbitrary) and the second associated with the event that the i-th particle suffered no collision with Bin the first (m-2) steps but defin i tely suffered a collision with Bin the (m-1) -th step. -150- In the reaction A+ B + C, particles of species B behave i n every respect like particles of A. Therefore the marginal density of finding the i-th particle at rr and of species B i s read i ly seen to be p( r~ , B) = ( l-p.) l I IAm-2 I I oo .- ldRm.-2dRm.- 2 p( r m . , r.m- 1 , m r . - 1 , R.m-2 ,Rm. - 2) drm.- ldrm 11 pj(l-pi) II 00 I l J lJ l l l l l l l J III 8m- 2 Am-2 m m-1 ,r.m-lm-1 p ( r.,r. ,r . . ,Rm-2 . ,R.m-2 ,Rm-2) .. l l oo lJ 00 oVj drm.-ldrm.-ldrm.-_ldRm. -2dRm.-2dRm.-_2 l J lJ l J lJ (4 . 16) By the hypothesis of uniform composition of particles, the i-th particle found at r~ can be either of species A, B or C (provided that A+ B + C is the only reaction possible). Therefore it .is clear that the probability den sity that the i-th particle is of species Candis located at r~ at time t = m6t is l Eqs . (4 . 15- 17) are our working equations . It should be kept in mind that they are conditioned on the initial distribution of particles (r 1, ... , r N) . Once the opdf is given , the concentration statistics of reactive materials undergoing a second-order reaction in turbulence can be evaluated. -151The formulation given above can easily be extended to the cas e where particles are not released at the same moment . Th i s is done by redefining the inert opdf as m-1 rN m- 1 ' ••• ' (4. 18) where Nm i s the total number of particles present i n them- th time step . Obviously, among the Nm particles present, Nm- Nm-l are released in th i s current time step . The corresponding formulas based on this opdf are essentially the same as (4 . 15- 17) . Only the formula for p(r~,A) is presented below : for i = 1 , ... , Nm- l for i = Nm- l +1 , . . . , Nm (4 . 19) where pis defined in (4 .18). Because the opdf's defined in (4.1) and (4.18) are virtually impossible to determine, their forms must be approximated . This is done in the following section . Several Lagrangian models of turbulent chemical reactions are derived thereby. -152IV.2 Lagrangian Models of Turbulent Chemical Reactions Before we begin to derive the models, it is appropriate to introduce the following formulas from Lamb (1972a) which relates the opdf s to the concentration statistics. 1 Let c(r,t) denote the instan- taneous concentration of some species measured at location r (r is a position vector) at time t, and let angle brackets denote ensemble average, then <c(r,t)> (4.20a) <c 2 (r,t)> = lim (4.20b) tN+O where ~Vis the sampling volume centered at r. Since the sampling volume is small, we can assume that p(r~) and p(rr,rj) do not vary much in ~V. This leads to (4 . 21a) <c 2 (r,t)> = N N I I (4 . 21b) i =l j=l Hi If the total number of particles is large, the summation over all the particles can be replaced by an integral and the formulas then read: -15300 <C(r,t)> = It 0 J p(r,tjp,T)S(p,T) dp dT (4.22a) -oo -Itf f I 00 <c 2 (r,t)> 0 (4.22b) -oo where p(r,tjp,T) is the pdf that an inert particle released at location p and time Twill be found at location rand time t. duced earlier in (3.16) with similar notation. Eq. (4.22b) was intro- These formulas will be used frequently in the subsequent analyses. Two-Step Model In what follows we shall restrict outselves to particles of reactive materials immediately after they are released. Setting m = 2 in (4.15), we have r.+½oV N - I p.(l-p.) J JJ p(rf,A) = p.p(r?) l J l l • l J= j;t2 00 p(r?,r!,r~!r~,r~) dddr~ l r--½oV J lJlJ lJ (4.23) where the dependence of p on ri for i f i, j has been removed by integration. Unlike the general formulas (4.15), the information required 2 • equa t.10n, v,· z p( ri,ri,rj 1 1 Iri,rj , f or th 1s O O ) ,· s more rea d,· 1y acces,•bl e and hence the concentration statistics of reactive materials can be obtained. Summing the whole equation of (4.23) over for i = l ,2, ... N, and applying (4.21), we obtain -154- N N I I i =l j=l (4 . 24) i;j where 2 l ,r.11 r.,r. o o) p -- p( r.,r. 1 1 J 1 J r=r?1 t = 2£1t and AI is the concentration that A would have if it were inert. Since (4.24) is valid only for O < t < 2£1t and the magnitude of £1t is restricted by condition (4.6), namely £1t << (k<cmax>- 1 ), applications of this model to fast reactions are severely limited . However, in these cases (4 . 24) might prove to be of great value in resolving the initial phases of reaction that many closure schemes are unable to describe . For this reason we expect that the two-step model is useful only in conjunction with closure schemes that normally do not apply near the source . Homogeneous Model The second model derived from the Lagrangian description is the homogeneous model . This model considers the case of a dilute mixture of substances A and B reacting to form C(A + B + 2C) in a well stirred, continuous reactor . It is assumed that the cell is continually stirred so that i f any foreign substance is introduced into the cell -155at any time, it becomes thoroughly mixed in the reactor within a time Tm. Suppose that the time scale TR in (4.6) of the reaction is suffi- ciently large that an interval 6t exists such that (4 . 25) This value of 6t enables us to use formula (4. 15-17) because it satisfies (4.6). Since 6t is much larger than the mixing time Tm, it is expected that the position of any particle at time n6t is statistically independent of its position at any other time m6t, m 'f n, and that the positions of the particles at any given time are statistically independent of one another. These two results yield the following probability relationships: = p( rim) p( rim- l ) p ( rim-2) (4.26) and ( 4. 27) respectively. Furthermore, when the above conditions are satisfied it is logical that the marginal densities p(r~) of each particle will be equal and will be constant in space and time. Thus, we assume -156- p(r~) = q = { l r.m V ' l inside V, i = 1,2, ... , N (4. 28 ) m=l,2, ... , r~ outside V, 0, where V is the volume of the reactor . Hence 9., ctr.l 00 We now turn to the evaluation of (4.15) . By applying (4.26-28) it is easy to show that (4. 15) reduces to ( m- l ,A;rjm-1 ,B) oV. pr; J (4 .29) where ( m- 1 m- 1 ) .. ) ( m-1 ,r.m-2 ,R.m-2 ,Rm-2 pr; ,A;rj,B = I I I pr. l J J lJ oo Am- 2Bm-2 dr~- 2ctRr~-ZdR~~ 2 l J lJ (4.30) f which in this case is simply J J Am-2 8m-2 00 Summing (4.29) over i for i = 1 to NA, and applying (4.21) we have <A(r,t+~t) > = <A(r,t)> - oV<A(r,t)B(r,t) > ( 4. 31) -157where t = (m-l)6t and the formula (4.21b) has been generalized to include the case of two species. k Since oV = a<v 12 > 26t by definition, we can rearrange (4.31) to give (4.32) Taking the limit of (4.32) as 6t approaches zero, we have a<A(r,t)> = - k<A(r,t)B(r,t)> at (4.33) L where k = a<v' 2 >~ One interesting point is observed here is that in the case of turbulent mixing of inert substances, the concentrations will be uncorrelated if the particle motions are statistically independent. This can easily be seen by applying the condition (4.26-27) to (4.20-4.22). One finds <c 2 (r,t)> = <c(r,t)> 2 Yet in the cases where the materials are reactive, we see from (4 .33) that statistical independence of particle motions (a condition of the turbulence) does not necessarily cause uncorrelated concentrations . The concentration correlation here is caused by the chemical reaction itself. When two species react to form another species at certain point, the concentrations of these two reactants instantly become smaller than -158concentrations at the vicinity of this point. tuations or correlations are created. Hence concentration fluc- Since reactions are occurring at the molecular level, only molecular diffusion can destroy concentration correlations of this kind. Let T0 be the time scale associated with molecular diffusion. (Note that T0 could be greater than Tm, the turbulent mixing t ime scale in (4.25) . ) If TR>>T0 , molecular diffusion destroys the concentration correlation much faster than the reaction is able to create it, then we have a<A (r,t) > at = a at -k <A(r,t) ><B(r,t)> (4.34) Eq. (4.34) is the familiar expression often used in chemical kinetics in describing the rate of a second-order reaction A+ B ➔ C in a well-stirred tank reactor. Many applications of this formula can be found (e .g. Vas silatos and Toor, 1965). For reactions with higher speed, this formula often fails to apply_ (See Glassman and Eberstein,1963.) Markov Process In deriving the homogeneous model, the particle positions at any time t were assumed to be statistically independent of their positions at all previous times. In real turbulence, however, particles "remember their past histories. To account for this characteristic 11 of turbulent motion, we may assume that particles remember their -159- trajectories during the period 6t immediately prior to any given instant but not anything earlier. tion. This is commonly known as the Markov assump- If Tis the Lagrangian time scale of the turbulence, the Markov assumption implies 6t > T (4.35) With this assumption, we have m m-1 rri-=-T ,R.m-2 ,R.m-2 ) -_ pr. ( ml r.m-1 ,r. m-T) pr. ( m-1 ,r ITI-1 p(r.,r . ,r. . ,R.m-2 ,R.m-2) 11 1 1 1 11 1 1 1 1 1 (4.36) Since particles do not interfere with one another's motion, 1 ,r.m- 1 )=pr.r ( mI mp(rm . Ir.m.1). 1 1 1 1 1 ( 4. 37) Applying (4.36) and (4.37) to (4.15), integrating, and using definitions (4.12) and (4.30), we obtain = . )p(r.m-1 ,A) dr.m-1 f p(r.m rm-1 1 1 1 1 00 N 1 )p(r.m- 1 ,A;r.m- 1 ,B)dr.m- 1dr.m- 1 - ·=N I +l J f p(r.m1 Jr.m1 1 J 1 J J A oo (4.38) o . J We note that the transition probability density p(r~lr~-l) is independent of chemical reaction . If the turbulence is homogeneous and station- ary, we have m m-1) = ( mIr.m-1 ) = p ( r.-r. pr. l l 1 1 -160- where 6xk(k = l ,2,3) are components of the vector 6r = r~-r~-l. l l The function p( 6x1 , 6x , 6x3 ) is determined by the state of the turbulence. 2 Substituting ( 4.39) into ( 4.38), summing the equation over i = l to NA and applying (4.21), we have (prov ided that the reaction volume 6V is sufficiently small) <A(r,t+6t) > = f p<A(r-6r , t) > d(r-tr) 00 - av f p<A(r- 6r,t)B(r-6r,t)> d(r- 6r) 00 where t = m6t and r = r~. Eq. (4 .40) can be expanded into a Taylor series about the point (r,t): <A(r,t)> + 6t a<A(r,t) > + 0(6t) 2 = at f p{ <A(r,t) > 00 f --o V p{ <A(r,t)B(r,t) > - 6x k a<A(r,t)B(r,t) > + 0( 6x ) 2 }d( 6r) (4.41) 3Xk 00 where d( t r) = d( t x1 )d( 6x2 )d( 6x3 ). The moments of the probability den- sity p = p( tx 1 ,6x 2 ,6x 3 ) are defined by l = f p d( 6r) (4.42a) 00 (4.42b) -161(4.42c ) 00 Substituting (4 . 42) into (4 . 41), we have a2 <A(r,t) > axkax i (4. 43 ) The moment <6x k> is the k-th component of the mean displacement of a particle in a time 6t. If the turbulence is homogeneous, only a mean fluid motion can produce a non-zero <txk >. Let the mean fluid velocity be a constant and have components <Uk>. Then The second moments <tx kt xt> are influenced by the turbulent velocity fluctuation and can be written in the form (4.45) where Rk ,., 0 = lim Jt{ < v 1 (t)v 1 (t-t 1 ) + v'(t-t 1 )v 1 (t)>} t -+<xi k t k t dt 1 (4 . 46) 0 and vk is the Lagrangian velocity of a particle measured with respect to the coordinate system movin g with velocity (<u1 >,< u2>,< u3>). -162Substituting (4.44-45 ) into (4.43) and taking the limit of 6t ➔ O, we get l 3<A(r,t)> 32 <A(r,t)> - k <A (r,t)B(r,t)> 3<A(r,t)> + <Uk> = 2 RkQ, 3t 3Xk3XQ, 3Xk (4 .47) h:2 where k = ov/6t = a<v' 2> as in (4 .33). In actuality, the limit 6t ➔ 0 cannot be achieved because of the restriction (4 . 35). To approach this limit, we should first normalize the coordinates tot*= t/T and xk = xk/Lk' where T and Lk are characteristic time and length scales, respectively, of the concentration <A>. We then let 6t* ➔ 0 in the limit of large T and Lk. That is to say, (4.47) is valid when T » 6t (4. 48a) (4 .4gb ) Combining (4.35) with (4.48), and assuming the Lagrangian length scale L of the turbulence to be <v' 2 >' 2T, we have T » T (4 . 49a) (4 .49b) Also, conditions (4.35) and (4.6) can be combined to give T k<c max >« l Equation (4.47) is similar to the Eulerian diffusion equation of reactive species which we introduced frequently in the previous -163chapters . The LHS of ( 4.47) is the convective t r ansport of spec i es A, while the effects of turbulent and molecu l ar di ffu sion are both i ncluded in the first term of the RHS . the reaction decay. The second te rm i n t he RHS i s Because this term repres ents an add i tion al un- known quantity , Eq. ( 4.47) is not closed . Through an ext ention of Prandtl 's mixing-length hypothesis , Lamb (1973) showed that ( 4 . 51) <A(r , t)B(r,t) > = <A(r , t) >< B(r , t) > provided that conditions ( 4.49-5~ are satisfied and that the r eaction time scale TR is much larger than the diffusion time scale TD . A special case of (4 . 47) with assumption (4.51) was solved in Chapter I II . In general, the Markov model can only be used for slow reactions at distance far from the source. IV.3 Discussion We have demonstrated the Lagrangian techniques of describing turbulent chemical reaction. Three simplified models were presented but the applications of each appear to be limited. Although the La - grangian description gives a better picture of the physical phenomena than t he Eul er ian methods , its use is greatly hindered by the l ack of information of the opdf. Only in a limited number of cases the opdf be measured or hypothesized (e . g. , Chapter III . ) can Efforts were made by the author to numerically measure the complete opdf in a one-dimensional turbulence simulated by the Burgers' equation (Bur gers, 1948) . Because Burger's model pertains to compressible flow, -164- it was found that particles tend to cluster together rather than separate. Consequently the result of such studies are of little direct value to problems associated with incompressible fluids. It therefore appears that the Lagrangian approach to the phenomena of turbulent chemical reaction has rather limited application at present. Experimental techniques to measure Lagrangian properties of turbulence will probably remain very difficult. However, with the present speed of growth in computational technology, it may soon become possible to obtain a complete opdf through numerical simulations. Further development of the Lagrangian method therefore will rely on the advancement of numerical simulation of turbulent flows. -165- V. CONCLUSION AND RECOMMENDATIONS FOR FURTHER WORK In the thesis we have demonstrated the importance of the effect of inhomogeneous mixing on the rates of chemical reactions in turbulent flows and have developed what we call a diffusion zone (DZ) model for describing these effects. Turbulent mixin g (macromixing), molecular diffusion (micrornixing), and chemical kinetics have been incorporated into the diffusion zone model. It was shown that the model gives a realistic description of the phenomenon of turbulent chemical reactions in chemical reactors. The model was also applied to the study of the N0-0 3 reaction in a powerplant plume in the planetary boundary layer. The prediction of the model showed qualitative agreement with observations in a case study. We also developed a Lagrangian model of turbulent chemical reactions but its applicability is much more restricted than that of the diffusion zone model. Following is a list of conclusions drawn on the work of this thesis and recommendations for future work. 1. The rates of decay of the mean concentrations of species involved in second-order isothermal reactions in a turbulent fluid can be predicted using the so-called diffusion zone model developed here. The model requires as inputs the mean and mean square concentration of inert species released from the same sources in the same turbulent fluid. -1662. The diffusion zone model can be applied to turbulent chemical reactions of any speed . However, in the limits of very slow and very rapid reaction, some simplication may be made. In the case of slow reactions, the mean concentrations alone are adequate to determine the reaction rates. In the case of very rapid reaction, Toor's model may be used for reactor design purposes . The Damkohler number a, defined in (2.42), provides a measure of the reaction speed. The ranges of Damkohler number to which each of the models just mentioned is applicable are summarized in Table 5.1. 3. The inert concentration statistics required by the DZ model can be obtained in one of the three ways: direct measurement, measurement of the conversion rate of a very rapid reactions (Toor, 1962) or by numerical simulation. A numerical technique for measuring the mean-square concentration of an inert species was introduced in this work. 4. The effect of inhomogeneous mixing on the rate of the NO-O 3 reaction in the atmosphere appears to be large. It was shown that the segregation of NO in a point source plume limited considerably the rate of reaction of this species with ambient ozone. 5. The 50% conversion time-scale of NO emitted from a stationary point source into an ozone laden atomosphere was pre- Intermediate or Rapid Very Rapid (Diffusion-Controlled) 0.01 - 50. 50. < * See (2.42) for details. Very Sl O\'J (Kinetics-Controlled) <0. 01 Definition of the Reaction Regime ( *) max{ OA,VG } o: = k /min{a 0 ,b 0 } Toor(l962) DZ model - k< A> Reactor Analysis DZ morlel DZ model - k<A> Atmo spheric Reactions Recommended Model for -k<AB> Table 5.1 Recommended Use of Models for Turbulent Chemical Reactions of the Type A+ B ➔ P I --...J CJ') I ,_. -168- dieted by the DZ model. The prediction is in good agreement with observations. 6. The Lagrangian approach is hindered by two main problems: the lack of information regarding the trajectories of fluid particles in turbulent fluid, and the closure problem that arises when materials that undergo nonlinear chemical reactions are treated. The first problem can be alleviated to a large degree by simulating particle trajectories in turbulence numerically. A technique for simulating the trajectories and joint displacement probability density function of two particles in a three-dimensional turbulence was developed in this work. The closure problem associ- ated with the Lagrangian method is subject to treatment by any one of the Eulerian schemes described in Chapter I. In summary, we have studied turbulent chemical reactions in general and have developed a useful tool, the diffusion zone model, for treating this phenomenon mathematically. The capability of the DZ model to predict decay of reactive species in turbulent flows has been tested using both laboratory and atmospheric observations. It is recommended that the possibility of applying the DZ model to more complex kinetics be investigated; with an expanded -169- capability of this type it would be possible, for instance, to examine the effect of turbulent mixing on the selectivity *or the product distribution in a reactor (Brodkey, 1974), or in the area of air pollution modeling, to study controversial issues such as ozone build-up in the far end of a power-plant olume (Davis et al , 1974). Some ideas of this extension are discussed in Appendix D. A few other studies may give improvements of the pre- dictions of the DZ model. Experimental measurements of the diffu- sion zones is one example. By choosing proper sample sizes , con- tinuous measurements of the instantaneous concentrations of two inert species mixing in a turbulent flow would provide a quantitative picture of the diffusion zones . The mixing parameters sA, sB could then be computed from observed data. Another problem of interest is the trentment of mixing fluid elements of various 1 ages 1 • This problem arises in studying atmospheric reactions among pollutant species emitted at different times. This problem is also briefly discussed in Appendix D. * The selectivity problem refers to the fact that in a complex reaction system, a species is often involved in more than one reaction, and that depending on the reaction condition, such as temperature or turbulent mixing, the species may be consumed by certain reactions but not others. -170- APPENDIX A. MULTIJET REACTOR DATA The purpose of this appendix is to give a brief description of the multijet reactor data, and how we used the data to calculate correlations like rAB in Chapter II. As noted in Chapter I, Toor (1962) developed a theory which relates the concentration statistics of two species undergoing a very rapid irreversible reaction in a trubulent mixer to that of an inert species in an identical mixer. To test the theory experimentally, Vessilatos and Toor (1965) designed an ideal one-dimensional turbular reactor the head of which is made of some 100 small nozzles. (Figure A.I) Reactants are fed through alternate jets to simulate a cross-sectionally uniform concentration profile. work (Mao, 1969) . The reactor was modified in a subsequent Both Vassilatos and Mao made extensive measurements in the reactor, covering a wide range of reaction speeds. However, their results only offer a partial test to Toor's theory because no inert mixing data were taken. This inert mixing data as well as axial velocity profile were later supplied by Brodkey's group (McKelvey, 1968; Zakanycz,1971). Toor's theory has since been shown to be good by these works and some others (Keeler et al, 1965; Torrest and Ranz, 1970). The reactor is then used to test other hypothesis (e.g., Mao and Toor, 1971). It appears that data relating to this reactor are the only set of data available which contains the minimal required information to test models for turbulent chemical reactions. Figure A.I •c 0 I I I I I I STREAM B I ◄ ~ ~ f."¼ -'.•· )II,, t = X/U ~{) r \ ! } ~ 0& J-~ · ·-·::· •::•.•..•. ·.•·•.•.· . ... ~~~ p~ /'- START REACTION Schematic diagram of the multijet reactor. The reactor head is composed of some 100 small tubes. Reactants are fed through alternate tubes to simulate a cross-sectionally uniform concentration profile . (Source: Hill, 1976) REACTOR HEAD A I ~ I --.J ...... 1 ...... -172- Ideally, we can directly apply the inert mixing data of McKelvey (1968) to our diffusion-zone model and test the model prediction against reaction data . Unfortunately , MeKelvey reported his mean square concentration fluctuation data in a normalized form without specifying the normalizing factor. An alternative is to use Toor s 1 theory to convert the data for very rapid stoichiometric reaction to the inert mixing data. This conversion is reliable because Toor 1 s theory has been well established in this reactor. However, it rules out the possibility of comparing the DZ model with Toor 1 s model in the case of very rapid stoichiometric reactions. As noted in (1.23), the theory of Toor (1962) says that for very rapid reactions with stoichiometric feed of two species F (x) 00 = <A(x)> k <AI BI ( x) > 2 - <A(xo) > -- •- - - ~ ! , : <A B1 (xo)> 2 1 (A . l) where x=<u(x)>t is the axial position along the reactor length and xo is some reference point where the reactor has reached the state of cross-sectional homogeneity . Toor and his coworkers have designed the multijet reactor in such a way that this state is approximately reached at the inlet of the reactor, and hence x =O . Toor (1969) further derived a relationship between the concentration fluctuations of two unpremixed inert species fed into this reactor: <A I BI ( X) > I I <A 1 8 1 ( 0) > I I (A.2) -173Since the concentration field is assumed cross-sectionally homogeneous, the initial concentration fluctuations can be obtained theoretically (Keeler et al, 1965) by averaging over the cross-section. If half of the nozzles feed A and half of them B, which is true in Mao's reactor and approximately true for Vassilatos' ,, it is easy to show that <A 2(o)> = <AI> 2 1 = <BI> 2 (A.3) Since the species are not premixed we have <A 1B1(o)> = -<AI><BI> (A.4) Combining (A.I) through (A.4), we obtain <AIBI(x)> rAB(x) - = l+ <AI><BI> <AI><BI> 2 <AI(x)> rA(x) - <Al> <A I1 BI1 (x)> = 1- Foo 2(x) (A.5) 2 = 1+ <AI (x)> <Ar> = 1+ Foo 2 (x) = r B (A.6) where <AI>' <BI> are independent of x, due to homogeneity. Relationships (A.5) and (A.6) are used to convert the concentration decay of very rapid stoichiometric reactions F (x) to the inert 00 concentration correlations rAB' rA and r 8 . Conversion data of reactions with rate constants greater than 10 7 t/g-mole-sec are not distingushable from one another, and are considered to be diffusion-limited (Mao, 1969). We shall use the mean concentration data of these reac- tions for 00 • Figures A.2-A.4 show the data of F for three different F flow conditions. 00 A cubic spline function has been used to interpolate data between measurements. in the model calculation. The interpolated curves of F are then used 00 -174- 0 ......." C B - - - - - - - , - - - - - - - , - - - - - - - . - - - - - - - . VERY RAPID REACTION STOICHIOMETRIC FEED CD 0 C9 (D 0 N 0 • VASS ILATOS CURVE FITTING . o.___ ____._ ____.__ __...._ __.___ ___.___ __,___ _.____ __, 0.0 1.0 2.0 3.0 4.0 X,CM Figure A.2 Extent of reaction of a very rapid stoichiometric reaction: Data of Vassilatos. -175- a •c-------.-----~-----~----- VERY RAPID REACTION STOICHIOMETRIC FEED (X) a • © (.D a • MAO , RE= 14:60 CURVE FITTING X a 0J a D • . o.____ ___...___ ____._ ___.._ _____.__ ___.___ ____.__ ___,__ __ 0.0 1.0 2.0 3.0 4.0 X,CM Figure A.3 Extent of reaction of a very rapid stoichiometric reaction: Data of Mao, lamiGar jets. -176- 0 - "ot---------.-------r-------r----------, VERY RAPID REACTION STOICHIOMETRIC FEED ~ 0 • (!) (0 0 • MAO, RE=2ll30 CURVE FITTING X C\I 0 0 • . o..___ ___..____ ___.._ ___.__ ___.___ _-L-..._ _ .L___ _L - -_ 0.0 1.0 2.0 3.0 ___. l!.0 X,CM Figure A.4 Extent of reaction of a very rapid stoichiometric reaction: Data of Mao, turbulent jets. -177In applying the data of F (x) to the Dz model calculation 00 where the independent variable is time, we have to know the axial velocity. Because of the jetting effects of the fluids coming out from the head, the velocity measured at the center of the nozzles is initially different from that measured between the nozzles (Yieh, 1970). Hence the use of a bulk-average velocity for <u(x)> will lead to considerable error. Yieh (1970) argued that in modeling Vassilatos 1 reactor, the velocity measured between the nozzles is more representative of the mean velocity because the area between the jets consitutes a major part of cross-sectional area of the reactor. In Mao's reactor, however, the reactor head is refined by using a larger number of nozzles, the velocity is hence more uniform. Following Yieh (1970), we shall use the velocity profile measured between the jets in our calculations. The axial velocity profiles for different flow conditions are shown in Figure A.5. -17 8- a a--------r---------.------~--------. N MEAN AX1AL VELCCITY -e- V/\SSILATOS a ...... CD -A-- MAO, Re= 14-00 -+- MP.O, Re= 2430 • "O xoo '-./ :::, a ____....____...____.____._____.___........._____,___ 3.0 L!.O 0.0 1 .0 2.0 X,CM Figure A.5 Mean axia l velocity prof iles of the multijet reactor in three different flow conditions. Data were measured between the jets. (Yieh, 1970) -179APPENDIX B. NUMERICAL CALCULATION OF REACTiON PARAMETER In this appendix the numerical method that was employed in calculating the reaction parameter \ will be discussed . We will also show that the minimum value of , is approximately independent of the choice of initial conditions for the model equations (2 . 25-2 . 31) . Three basic steps involved in calculating \ are: 1) Solve (2.32-33) numerically for a* and b* , 2) Determine the one-dimensional diffusion-zone LAB using the profiles of a 1 and b{ given in (2.39-40), 3) Calculate \ by (2 .41). Numerical solutions to a system like (2.32-33) have been studied previously by Brian et al (1961) and Mao and Toor (1970). Both of these calculations involved explicit schemes and a small time-step had to be used to insure numerical stability. In this work the quasilineari- zation sGheme of Lee (1968) is used since it allows a larger time-step This scheme uses O'Brien- to cut down the amount of computation time. Hyman-Kaplan implicit scheme for integration and Newton-Ralphson interation scheme for linearization and is briefly described below . Let anr,s, bnr,s denote the concentrations a*, b*, respectively , at t*=s~t, x*= r6x for then-th iteration. n+l ar,s+l n+l ar+l,s+l - a r,s - We represent (2.32-33) by n+l 2a n+l r,s+l + a r-1,s+l = n+l bn+l + bn+l br+l,s+l - 2 r,s+l r-1,s+l bn+l - b r,s+l r,s -------= n+l bn . - aB ar,s+l r,s+l y 2 6X for r= 1,2, .... , N-1 - a -180- where N is the number of x-increments and a r,s , br,s denote known concentrations of the previous time step . The x-increment increases with t such that the boundary conditions 0 ' for all n and s are always maintained. To proceed we first let a~,s+l = ar,s ' b~,s+l = br,s 2 2 and solve for ar,s+l' br,s+l using Thomas method for inverting a tridiagonal matrix . The iteration continues by using anr,s +l' bnr,s +l to n n I n+l n I I n+ 1 n calculate ar,s+l' br,s+l until ar,s+l - ar,s+l ~ E and br,s+l - br,s+ll ~ E, where E is any desired accuracy. The boundaries of the diffusion-zone LAB' in principle, in determined by where either aj or bj assumes a zero value, however, mathematical expressions (2.39-40) of aj and bj go to zero only in the limit of x*++ 00 • Hence in practice one has to determine approximate boundaries for LAB · Friedlander and Keller (1963), in treating a similar problem, assumed the boundaries as where Jdaj/dx*I~ 0.05. We found, however, to our satisfaction that LAB can be simply defined as where aI 1 =-- ~ C , a ao b* > c/s I ""' where ca 0 corresponds to some lower limit of a concentration measuring unit. Tests indicated that the minimum value of>-. varies very little for c < 0.01. A value of c=0.001 was used in our calculation. In solving (2.32-33), the initial conditions (2.36-37) have been used. However, it was discovered that the value of interest, namely >-.min' only depends weakly on the choice of initial conditions . In Figure -181B.1 A vs. t* with a=700, s= y= 1 is plotted for three different initial conditions. The solid curve is the result of y using (2.36-37) as initial condition. If initially the two pulses are next to each other, y decays earlier from unity as indicated by the dot-dashed curve. Lastly, the pulses are broken in half and alternately arranged. The dashed curve of A first decays to some pseudo Amin . which corresponds to the value Amin would assume if the two pairs of pulses were far apart, or only one of them is present. ( Since the pulse width in the last case is half of the original, the pseudominimum of A is the Amin for the case where a= 700/4. ) As the pulses gets wider, the interaction of neighboring pairs brings A to its true minimum. The important aspect of these calculations is that the initial condition influences the behavior of A fort ~o 2 /D, but only weakly affects the value of Amin" This supports our assumption in Chapter IIthat Amin is a function of a only. A . Figure B.1 10· 3 10 -If II I I I I 1111 b \\ 10-1, I I I I II I II t: I ( <a2/D) 10° I I I I 11111 J]_ 0.. I 10 1 I I I II Ill --~ I I I II 1111 Y=t r=1 cf...= 100 •• I I I 1111 II Temporal behavior of the reaction parameter\ for various hypothetical initial concentration distributions . The minimum value of\ is only weakly dependent of the choice of initial concentrations . 10-2 I I I I I II II lITli a. 10-3 \ "-.__ "-...... - ~'\. ... - - - - - - - - - "'- . . Do ill LEGEND 0.. ."'- "- -- --,~~-\-··-·· -~,, .... -~-- o.. • ·--·--....... 10-2 1 o- 1 10° I N co I ,_. -183APPENDIX C. ESTIMATION OF NO-EMISSION FOR A POWER PLANT This appendix describes how we estimated the NO-emission for Morgantown power plant . Efforts made to obtain direct measurements of NO-emission for the plant were in vain . We therefore use the following methods to estimate the NO-emission: 1. According to a Morgantown pl ant engineer, the plant consumes 225 tons/hr of coal and 40.5 tons/hr of fuel oil currently (October 1975). Using a power rating of 1000 Mw and typical values of heat of combustion, 13000 Btu/lb of coal and 18000 Btu/lb of fuel oil, we calculate the plant efficiency to be 23.4% . During the time of flight observations (Summer 1974), the plant burned a fuel mixture of 25% coal and 75% oil. required. Thus 109 tons/hr of coal and 327 tons/hr of oil were The plant engineer estimated a 1.9% sulfur content for both fuel oil and coal. If the burning was complete the estimated so 2 - emission is 16.5 tons/hr. D. D. Davis (private communication) reported that the concentration of NO was always about half of that of so 2 during the flights. 2. We therefore estimate that NO was emitted at 16500 lb/hr. R. C. Flagan (private communication) indicated that power plant NO-emission is about 0.7-1.2 lb (as N0 2 ) per Mega Btu (Thermal). Using the estimated efficiency of 23.4%, Morgantown plant inputs 4 Mega Btu/sec. The N~emission is hence 10000- 17000 lb tas N0 2) /hr or 6500-11000 lb (as NO)/hr. The figures, although still high, appear to be more reasonable.(For reference, California power plants output NOx ranging from 60 to 2000 lb/hr while the current goal is 15 lb/hr for 1976.) -184- APPENDIX D. CONSIDERATIONS OF MULTIPLE KINETICS In many reaction systems, a chemical species of interest is often involved in reactions with more than one species. The application of the diffusion zone model to such systems is discussed here. We assume that the chemical kinetics of the reaction system (often determined by well-mixed reactor data) are given. The problem is to apply the diffusion zone model to studies of the effects of inhomogeneous mixing on the rates of turbulent chemical reactions in complex systems. Consider a system of N chemical species: A1 , A2, combination, each species could be involved in a maximum of N bimolecular chemical reactions. That is, for the i-th species Ai the N reactions could be A; + A.J kij N 4 J-c 1 at At j = 1, ... ' N (D.1) ~;ti,j .. is the rate constant and a Q, for i = 1, ... ,N. Where in (D.1), klJ is the stoichiometric coefficient of the reaction. If Ai is inert to Aj k;j is zero, and if Ai is not a product species of the i-j reaction at assumes a value of zero. Let Ai also denote the concentration of species, the instantaneous local rate of change of Ai can be written as ~N .. A.A. + R. dt - - E klJ 1 J 1 j=l (D.2) The term R1 in (D.2) denotes all other kinetic terms influencing the decay of A1 , e.g., production of Ai by other species and decay of Ai due to -185- reactions of orders other than two. The average concentration of Ai 1s defined as -A. = -1 1 V f A.dv (D.3) 1 V The rate equation of Ai is then dA. dt N --1. = - L j=l k .. A.A. lJ 1 J (D.4) + R. 1 Assuming the diffusion zone vij is defined as the region of space where A1 and Aj coexist, (to be further discussed) we can write -A. =-(v.A. 1 /\ + v . ."a.) 1 V 1 1 lJ 1 f\ " where v is the region where A exists excluding v,..J I and A. and~1 1 . 1 . 1 are average concentrations in v. nad v .. respectively. Similarly, lJ 1 letting the subscript I denote inert concentrations, we have -A. 11 1 " . + v .. a. /\ ) =-=-(v.A 1 V 1 1 lJ 1 Following through similar arguments as in the case of a single reaction, we arrive at ~= - dt N L k .. >.. .. j=1 z; .. z; .. lJ lJ lJ Jl r lJ .. [A.-(1-z;; .. )A. J[A.-(1-z;; .. )A. J 1 lJ 1 1 J Jl J 1 + R. 1 (D.5) -186- where >.. . . = lJ I\ Vi/jI 1;; .. = - 1 J V Aj I Although (D.5) gives us a general relationship for the decay of concentration of a species involving in N second-order reactions and other reactions, the parameters contained in the equation are rather difficult to estimate. The difficulty arises mainly because the diffusion zone is not as well-defined as in the case of a single reaction. To elucidate: the problem, let us consider two special cases: Case 1: Case 2: Al+ A2 ~ A3 A3 + A4 ~ Products Al + A2 Al+ A3 --- Products ~ Products In the first case, the diffusion zone v34 is defined as the region of -187space where A3 and A4 would coexist if they were inert; but since A3 is a product species such a region of space is not defined in the inert mixing problem. In the second case, the application of the diffusion zone model is hindered by the lack of a generalized model of the reaction parameterA. The problem here is that the diffusion zones v12 and v13 share the volume of space v123 and consequently the parameter A12 and A13 are not necessarily describable in terms of the material distributions in v12 and v13 as assumed by the model of A that we develop~d ' in Section II.2.2. It is easy to image other situations where the diffusion zone is not well-defined. We hence conclude that the application of diffusion zone model to complex reaction systems depends on the generalization of the diffusion zone concept. For example, we may note that in the first case the product species A3 can exist only in v12 . Thus, as a first approximation, one might assume that where v124 is the region of space where A1 , A2 , A4 coexist if they were inert under the same mixing condition. In general, however, the problem remains to be a challenge to the future workers. A somewhat similar problem arises in treating reaction systems where species are released at different times. One example is the emission of fresh reactive species into an atmosphere laden with both product and reactant species. In cases where ambient fluid is void of the emitted -188- species, say A , then the diffusion zone model may be applied in the manner developed in this work . However, if the emitted species A is already present in the background fluid, then a problem arises in the description of the reaction parameter A. For example, the fluid into which the fresh species is emitted can possess diffusion zones in any of the following three stages: (a) Initial per i od, >..= 1 A _c-1 -_ct _B (b) Reaction period, A= A . mm (c) Completion per iod , A= 1 B -- -- - - - ·- - A Here A and B denote ambient species concentrations. As we showed in Ch. II, in cases where reactive species are released only at some initial instant t o , the parameter has a value Amin . during most of the period of reaction. It was based on this observation that an expression for A was developed. Needless to say, the model of A is not necessarily valid if reactants are injected into the system at some time after t 0 . In view of the difficulty of developing a A-model applicable to all possible conditions of reactant releases , it is suggested that the diffusion zone model be used in the cases where the ambient fluid is well-mixed at each time that fresh species are introduced. Inasmuch as the class of reactions that we have considered in this work are prevalent in both engineering and atmospheric pollution studies, the diffusion zone model developed here should find widespread uses. -189- REFERENCES Ajmera, P. V., Singh, M., Toor, H. L. 1975 . Chem. Eng . Commun. Submitted for publication . Batchelor, G. K. 1949. Aust. J . Sci . Res . 2:437-50 1950. Quart. J. Roy. Met. Soc . 76:133-46 1959. J. Fl. Mech . 5:113-33 Brian, P. L. T., Hurley, J. F. , Hasselti ne, E. H. 1961. Borghi, R. 1974. Adv. Geophys. 18B :349-65 Brodkey, R. S. 1974. ed. AIChE J . 7:226 Turbulent Mixing in Non-Reactive and Reactive Flows , S. N. B. Murthy, 225-26 . Project SQUID Workshop, Purdue Univ., West Lafayette, Indiana, May 1974 (To be published). 461 pp. Burgers, J.M . 1948 . Adv. Appl. Math. 1:171 - 99 Bush, W. B., Fendell, F. E. 1974. Carslaw, H. S. , Jaeger, J. C. 1959. Acta Astronaut . 1:645-66 Conduction of Heat in Solids. Oxford, 510 pp. Chung, P. M. 1969 . AIAA J. 7:1982-91 1970. Phys . Fluids 13 :1153-65 1971. AIAA J. 9:2282-84 1972 . Phys . Fluids 15:1753-46 Clarke, R. H. , Dyer, A. J., Brook, R. R., Reid, D. G., Troop, A. J. 1971 . Boundary Layer Data. CSIRO Div . of Meteorol. Phys., Tech. Paper No. 19, 340 pp. -190Corrsin, S. 1958. Phys. Fluids 1: 42-47 1961. J. Fluid Mech. 11:407-16 1962 . Proc. 1961 Symp. on Fluid Dyn. and Appl. Math . , ed. J. B. Diaz, S. I. Pai, 105-24 . 1964. New York : Gordon and Breach. Phys. Fluids 7: 1156-59, 8:1920 Davis, D. D., Klaaber, G. 1975. Proc. Symp. on Chem. Kinetics Data for the Lower and Upper Atm. Warrenton, Va., Sep. 1974 (in International J. Chem. Kinetics 1975) Davis, D. D., Smith, G., Klauber, G. 1974 . Science Nov. 22, 1974. 186: 733-36 Deardorff, J. W. 1970a. J. Fluid Mech. 41:453-80 1970b. J. Atm. Sci. 27:1211-13 1972. J. Atm. Sci. 29:91-115 1973. J. Fluid Eng. , Sep. 1973, 429-38 1974a. Boundary Layer Meteorol. 7: 81-106 1974b. Boundary Layer Meteorol. 7:199-226 Drake, M.A. 1975. Boundary Layer Data from a Numerical Integration in Three Dimensions . NCAR report. Peskin, R. L. 1970. Phys, Fluids 13:584-95 Willis, G. E. 1974. Adv. Geophys 18B:187-200 Donaldson, C. du P., Hilst, G. R. 1972a. ~-- 1972b. Environ. Sci. Technol. 6:812-16 Proc. 1972 Hear Transfer Fluid Mech. Inst., ed . R. B. Landis, G. J. Hordemann, 253-61 Donaldson, C. du P., Rosenbaum, H. 1969. Compresible Turb. Boundary Layers, 231-53. NASA SP-216. -191Dopazo, C., O'Brien, E. E. 1973. Phys, Fluids 16: 2075-81 Evangelista, J. J., Katz, S., Shinnar, R. 1969. AIChE J. 15: 843-53 Flagan, R. C., Appleton, J. P. 1974. Combust. Flame 23: 249-67 Friedlander, S. K., Keller, K. H. 1963. Gifford, F., Jr. 1957a. --- 1957b. Chem. Eng. Sci. 18:365-75 J. of Meteorol. 14: 410-14 J. of Metheorol. 14: 475-76 Glassman, I., Eberstein, I. J. 1963. AIAA J. 1: 1424-26 Hawthorne, W. R. , Weddell, D. S., Hottel, H. C. 1949. Flame, Explos, Phenom . , 266-88. Hill, J. C. 1970. 3rd Symp. Combust., Baltimore: Williams and Wilkins. Phys. Fluids 13: 1394-96 Hill, J. C. 1976. Homogeneous Turbulent Mixing with Chemical Reaction. Annual Rev. Fluid Mech. Volume 8 Hilst, G. R., Donaldson, C. du P., Teske, M., Contiliano, R., Freiberg, J. 1973. A Coupled Two-Dimensional Diffusion and Chemistry Model for Turbulent and Inhomog. Mixed Reaction Systems. U.S. Environmental Protection Agency EPA-R4-73-016c.(PB-234193/1Gl) Hinze, J. 0. 1959. Turbulence. McGraw-Hill, New York. 586 pp. Johnston, H.S., Crosby, J. H. 1954. Kattan, A., Adler, R. J. 1967. --- 1972. J. Chem. Phys. 22: 689-92 AIChE J. 13: 580-85 Chem. Eng. Sci. 27: 1013-28 Keeler, R. N. , Peterson, E. E., Prausnitz, J.M. 1965. Kraichnan, R. H. 1970. AIChE J. 11: 221-27 Phys. Fluids 13:22-31 Lamb, R. G. 1972a. Stat . of Marked Particle Cone. in Turb. Fl ., Part I: General Theory; 1972b. Unpublished Note s. Stat. of Marked Particle ConG. in Turb. Fl., Part II: The Case of Chemically Reactive Particles: Unpublished Notes. -192Lamb. R. G. 1973. Atmos. Environ. 7: 257-63 Chen, W. H., Seinfeld, J. H. 1975. Numerico-Empirical Analyses of Atmospheric Diffusion Theories. J. Atm. Sci. 32:1794-1807. Lee, J. 1966 Phys. Fluids 9: 363-72 --- 1973. Quart. Appl. Math. 31: 155-76 Lee, E. S. 1968. Quasilinearization and Invariant Imbedding. Academic Press, New York. Libby, P. A. 1972. Combust. Sci. Technol . 6: 23-28 Lin, C. -H. 1971. Stochastic Chemical Reaction with Molecular Diffusion. MS theses. SUNY, Stony Brook, New York. 73 pp. 0 1 Brien, E. E. 1972. Astronaut. Acta 17: 771-81 Lindberg, W. R., Thompson, R. 0. R. Y. 1974. Quart. J. Roy. Met. Soc. 100: 608-23 Mao, K. W. 1969. Chemical Reactions with Turbulent Mixing. PhD thesis. Carnegie-Mellon Univ. , Toor, H. L. 1970. AIChE J. 16: 49-52 ---, Toor, H. L. 1971. Ind. Eng. Chem. Fundam. 10: 192-97 McCarthy, H. E. 1970. Second-Order Chemical Reactions with Turbulent Mixing. PhD thesis. McKelvey, K. N. 1968. Univ. of Delaware, Newark, Delaware. Turbulent Mixing with Chemical Reaction. PhD thesis. Ohio State Univ., Columbus, Ohio. 123 pp. Miyairi, Y., Kamiwano, M., Yamamoto, K. 1971. Miyawaki, 0., Tsujikawa, H., Uraguchi, Y. 1974. 7: 52-56 168 pp. Int. Chem. Eng. 11: 344-52 J. Chem. Eng. Jap. -193Manin, A. S., Yaglom, A. M. 1971. Statistical Fluid Mechanics. MIT Press, Cambridge. 0 1 Brien, E. E. 1966. 769 pp. Phys. Fluids 9: 1561-65 1968. Phys. Fluids 11: 1883-88 1971. P~ys. Fluids 14: 1326-31 1974, Advan. Geophys. , Eng, R. M. 1970. Phys Fluids 13: 1393-94 Francis, G. C. 1962. Lin, C. -H. 1972. 18B: 341-48 J. Fluid Mech. 13: 369-82 Phys Fluids 15: 931-33 0rszag, S. A., Patterson, G. S., Jr., 1972. Phys. Rev. Lett. 28: 76-79 Pao, Y. H. 1964. AIAA J. 2: 1550-59 Patterson, G. S., Jr., Corrsin, S. 1966. Plasmas. S.I. Pai, ed. Rao, D. P., Dunn, I. J. 1970. Edwards, L. L. 1971. Riley, J. J. 1973. in Dynamics of Fluids and Academic Press, New York. p.275 Chem. Eng. Sci. 25: 1275-82 AIChE J. 17: 1264-65 Phys. Fluids 16: 1160-61 , Corrsin, S. 1971. Preprint, Conf. on Air Poll. Meteorology, Rayleigh, N. C., p.16. Patterson, G. S., Jr. 1974. Romberg, G. 1974. Phys Fluid 17: 292-97 Int. J. Heat Mass Transfer 17: 1163-79 Seinfeld, J. H. 1975. Air Pollution: Physical and Chemical Funda- mentals. 523 pp. McGraw-Hill, New York. -194Shea, J. R., III, 1976. PhD thesis. A Chemical Reaction in A Turbulent Jet. Calif. Inst. of Tech., Pasadena) Calif. (Submitted July 1, 1975) Tennekes, H., Lumley, J. L. 1972. A First Course in Turbulente. MIT Press, Cambridge . Thompson, R. 0. R. Y. 1971. Toor, H. L. 1962. 1969. Quart. J. Roy . -Met. 97: 93-98 AIChE J. 8: 70-78 Ind. Eng. Chem. Fundam . 8: 655-59 Singh , M. 1973. Ind. Eng . Chem . Fundam. 12: 448-51 forrest, R. S., Ranz, W. E. 1970. AIChE J. Uspensky, J. V. 19 3i 16: 930-42 Introduction to Mathematic Probability. McGraw-Hill, New York. p.265 Vassilatos, G., Toor, H. L. 1965. AIChE J. 11: 666-73 Willis, G. E. , Deardorff, J. W. 1975. A Laboratory Model of Diffusion into the Covective Planetary Boundary Layer. Submitted for Publication. Yieh, H. -N. 1970. Turbulent Mixing with Chemical Reaction. PhD thesis. Ohio State Univ., Columbus, Ohio . 135 pp. Zakanycz, S. 1971. Turbulence and the Mixing of Binary Gases. thesis. Ohio State Univ., Columbus, Ohio, 172 pp. PhD