Computational fluid-particle dynamics for the flame synthesis of alumina particles

Chemical Engineering Science 55 (2000) 177}191

Computational #uid-particle dynamics for the #ame synthesis of alumina particles Tue Johannessen , Sotiris E. Pratsinis
, Hans Livbjerg * Department of Chemical Engineering, Technical University of Denmark, DK-2800 Lyngby, Denmark
Institute of Process Engineering, Swiss Federal Institute of Technology, CH-8092 Zurich, Switzerland Received 10 November 1998; received in revised form 29 March 1999; accepted 30 March 1999

Abstract A mathematical model for the dynamics of particle growth during synthesis of ultra "ne particles in di!usion #ames is presented. The model includes the kinetics of particle coalescence and coagulation, and when combined with a calculation of the temperature, velocity and gas composition distribution in the #ame, the e%uent aerosol characteristics are calculated. The model is validated by comparison with an experimental study of the synthesis of alumina particles by combustion of Al-tri-sec-butoxide. Two parameters of the coalescence kinetics are estimated by regression of the model predictions to the measured speci"c surface area of the product particles. The estimated kinetics can be used to predict the surface area and shape of the particles for a wide range of synthesis conditions.  1999 Elsevier Science Ltd. All rights reserved.

1. Introduction Aerosol routes to materials synthesis are well established for manufacture of several commodities, and may have even greater potential for the future development of new materials with unique properties (Gurav, Kodas, Pluym & Xiong, 1993). Aerosol particles with high speci"c surface areas and primary particle size in the nanometre range can be generated by the combustion of volatile precursor vapours in #ames. Due to their size and shape, such particles have numerous applications in the manufacture of pigments, porous ceramic membranes and catalysts or catalyst carriers and other materials (Pratsinis, 1998). Their generation at high temperature may endow them with a much higher thermal stability than nano-particles made in aqueous solution by, e.g. the sol}gel techniques. The latter particles require additional post-processing for ceramic membranes or catalyst carriers in high temperature applications because their micro-structure is altered by recrystallization (Tranto, 1997). The mechanisms of particle formation in #ames comprise the combustion of the precursor, the nucleation of

* Corresponding author. Tel.: #0045-45-88-32-88; fax: #00454588-22-58.

particle clusters followed by the growth of particles by condensation and coagulation, while the coagulated particles sinter or coalesce towards compact shapes. It is often assumed that the precursor combustion is virtually instantaneous, and that the supersaturation of the condensing product molecules is so large that the gas to particle conversion is a purely collision limited coagulation of entities, which initially are single molecules but subsequently grow to multimolecular particles of continuously increasing size (Pratsinis, 1998; Wu, Windeler, Steiner, Bors & Friedlander, 1993; Lehtinen, 1997; Kruis, Kusters, Pratsinis & Scarlett, 1993). In this case the kinetics of coalescence and coagulation determine the ultimate shape and size of the particles in interaction with the complex #ame environment. If the rate of coalescence is slow, dendritic aggregate particles are formed, while a fast coalescence rate leads to compact, spheroidal particles. In this paper we present a mathematical model for the particle dynamics including the kinetics of coalescence and coagulation. The model is combined with a calculation of the detailed temperature, velocity and composition distribution in a di!usion #ame which allows a calculation of the e%uent particle characteristics to be made for di!erent vapour pressures and feed rates of the precursor and for di!erent #ame conditions. The model is compared with the results of an experimental study of

0009-2509/00/$ - see front matter  1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 0 9 - 2 5 0 9 ( 9 9 ) 0 0 1 8 3 - 9

178

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

Nomenclature a a Q c N D D D D GK d N E  f " g h jR G j G k k  k @ k A Kn m G m .+ N n N p q qR r R r A r N

particle aggregate surface area, m surface area of a completely coalesced (spherical) aggregate, m speci"c heat capacity, J/mol/K, or, J/kg/K particle di!usion coe$cient, m/s aggregate fractal dimension molecular di!usion coe$cient for species i, m/s primary particle diameter, m activation energy, J/mol dilution factor gravitational acceleration vector, m/s enthalpy, J/mol, or, J/kg turbulent mass #ux vector of species i, kg/m/s molecular di!usive mass #ux, kg/m/s turbulent kinetic energy, m/s parameter in expression for the characteristic coalescence time, m\/s boltzmanns constant, 1.308;10\ J/K thermal conductivity, J/m/s/K nudsen number mass fraction of species i mass fraction of the PM pseudo gas component particle number concentration per unit mass of gas, kg\ number of primary particles per aggregate pressure, Pa molecular di!usion #ux of heat, J/m/s turbulent di!usion #ux of heat, J/m/s radial coordinate, m gas constant particle/aggregate collision radius, m primary particle radius, m

the synthesis of alumina particles by combustion of Altri-sec-butoxide vapours.

2. Experimental 2.1. Apparatus The laboratory apparatus for the synthesis of Al O ,   shown on Fig. 1, consists of the burner and a saturator. The saturator is a 250 ml bubble #ask with the liquid precursor Al-tri-sec-butoxide (ATSB), through which nitrogen is bubbled via a porous disk distributor, which provides very small bubbles. They rise slowly through a liquid height of 10}12 cm, which ensures that the gas becomes saturated with ATSB. The ATSB is not as unstable and hazardous as other volatile aluminium compounds, e.g. AlCl and tri-methyl-aluminum, but its 

s S  S F S G t ¹ ¹  ¹ " u u V v v N w .+ z

Greek a b e j k k J o o N p s sR q D

position in (r,z) space speci"c surface area, m/g heat source term (chemical reactions, radiant heat absorption) net rate of production of species i per unit volume by chemical reactions time, s temperature parameter in expression for the characteristic coalescence time, 1400 K decomposition temperature for precursor, K velocity vector, m/s #uctuating component of the instant velocity in the x-direction, m/s particle/aggregate solid volume, m primary particle volume, m total particle mass #ow, kg/s axial coordinate, m the di!erential operator nabla Letters fraction of total particle mass #ow collision frequency function, m/s dissipation rate of turbulent kinetic energy, m/s mean free path of the gas, m gas viscosity, kg/m/s turbulent viscosity, kg/m/s gas density, kg/m particle density, kg/m surface tension, N/m stress tensor, kg/s/m turbulent momentum #ux tensor, kg/s/m characteristic coalescence time, s

vapour pressure is lower, and the saturation temperature must be 130}1703C to yield adequate production rates of alumina. The liquid temperature is controlled by a thermocouple inserted in a glass pocket extending into the centre of the bottle and connected to a temperature controller. To prevent recondensation of the precursor, the temperature of the burner and the tube line from saturator to the burner is held at a temperature well above the dew point. To avoid thermal decomposition the ATSB must be kept below 2003C and in an inert atmosphere (Haanappel, van Corbach, Fransen & Gellings, 1994). The ATSB decomposes in the presence of water vapour or oxygen, even at low temperature. The overall oxidation reaction is given by 2Al(C H CH(CH )O) #xO      PAl O #yCO #zH O.    

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

179

Fig. 2. Flame con"guration B and C.

Table 1 Operating conditions for the burner. Runs 1}3 di!er from 4}6 with respect to the N -#ow, while the methane #ow varies within each series.  The air #ow rate is 1.8 l/min, and the saturator temperature 1603C, corresponding to a mole-fraction of 0.0046 for ATSB in the carrier gas

Fig. 1. The experimental set-up. A gas ejector is used to continuously dilute and cool a small sample of product gas before measuring the size distribution of the diluted gas in the SMPS system.

The burner is similar to that of Vemury and Pratsinis (1995). It consists of three concentric quartz tubes. The innermost tube has an outer diameter of 4 mm and the spacing between the following tubes is 1 mm. All tubes have a wall thickness of 1 mm. Fig. 2 illustrates the two di!erent burner modes used. They are labelled B and C, following the notation of Pratsinis, Zhu and Vemury (1996). In mode B there is nitrogen and ATSB in the central tube, air in the second tube and methane in the third (double di!usion #ame) and in mode C methane in the second tube and air in the third (single di!usion #ame). To stabilize the #ame, the burner is surrounded by a shield open only at the top, and the whole system is placed in a fume hood. The #ow rates of the feed gases to the burner are measured by calibrated rotameters. Particles are collected on 47 mm glass "bre "lters (Gelman Scienti"c) placed in a "lter holder "xed 20 cm above the burner. The product gas is withdrawn through the "lter by a vacuum pump. The speci"c surface area, S , is  measured by multipoint nitrogen adsorption (Gemini 2360, Micromeritics) at 77 K using the BET equation. Assuming spherical particles, one can calculate a BET equivalent primary particle size, d , from #2 d "6/o S (1) #2 N  where o is the solid density. Particles for transmission N electron microscopy were collected above the #ame on TEM-grids by exposing cold, clean grids to the hot product gas at the burner outlet for "ve seconds, using thermophoresis as the dominant deposition mechanism.

Flame

B1

B2

B3

B4

B5

B6

N (l/min)  CH (l/min)  Flame

0.38 0.38 C1

0.38 0.27 C2

0.38 0.145 C3

0.26 0.38 C4

0.26 0.27 C5

0.26 0.145 C6

N (l/min)  CH (l/min) 

0.38 0.38

0.38 0.27

0.38 0.145

0.26 0.38

0.26 0.27

0.26 0.145

The e%uent particles consist of large aggregates, each containing numerous, smaller primary particles. The aggregate size distribution is measured with an SMPS (scanning mobility particle sizer) system consisting of a di!erential mobility analyser (DMA 3071A, TSI) and a condensational particle counter (CPC 3022, TSI). The gas ejector, shown on Fig. 1, uses a high #ow rate of particle free dilution gas to withdraw a small sample #ow of the product gas through a capillary tube by the ejector action of the dilution gas. This sampling technique has the advantage of a low residence time ((20 ms) of the hot and concentrated aerosol in the capillary tube, and the subsequent, immediate cooling and dilution of the sample stream e!ectively quenches further particle sintering and strongly reduces the coagulation rate of aggregates. Several distributions were measured for each con"guration. The experimental design of twelve measurements, shown in Table 1, comprises a constant air feed rate, three levels of methane feed rate, and two levels of nitrogen, pairwise identical for the two burner modes. The feed rates were chosen so that a stable #ame, based on visual observation, is obtained. For the two highest methane #ows, the #ame is slightly sub-stoichiometric with respect to oxygen from the air feed. Additional oxygen is provided from the surrounding still air. The #ames are all blue in colour, indicating total fuel combustion.

180

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

For #ame C, there is a single #ame front between methane and the supplied air. Flame B has two #ame fronts, an internal one between the supplied air and the methane and an external one between the methane and the surrounding air. The di!erence between B and C is clearly visible. For selected operating conditions, the temperature pro"les along the vertical axis of the central jet and the radial peak temperatures as a function of height were measured with a 0.008 Pt/10%Rh-Rh, Type S, thermocouple (Omega Engineering) insulated with an alumina sheath tube (ID 1/32, OD 1/8). The last centimetre of the two thermowires up to the welded connection point was not insulated. The measured temperatures are corrected for radiation loss according to Collis and Williams (1959). All pro"les were measured without precursor. Measurements with precursor present results in the deposition of particles on the wire by thermophoresis and di!usion, which changes the emissivity of the wire. The combustion of the precursor has little e!ect on the #ame temperature due to the low precursor concentration in this study.

3. Theoretical

time-smoothed variables are the stress tensor s, the turbulent momentum #ux tensor sR, the velocity u, gas density o and pressure p. The vector g is the gravitational acceleration. The local value of the tensor sR is calculated by the di!erential Reynolds stress turbulence model: The individual components of sR, the Reynolds stresses, in cartesian coordinates are de"ned by qR "ou u ; qR "ou u ;2 (3) VV V V VW V W where u etc. are the #uctuating component of the instant V velocity and the bars denote time smoothing. The individual stresses are computed from transport equations derived from the momentum equations using the modelling assumptions of Launder, Reece & Rodi (1975) and Launder (1989). This yields six scalar partial di!erential equations to account for the six independent stresses of the symmetric tensor sR. The turbulent kinetic energy k is de"ned by (4) k"(u #u #u ). W X  V To account for the dissipation rate e of the turbulent kinetic energy k an additional equation is needed (Launder & Spalding, 1974):

A model for the process of particle generation, based on a detailed calculation of the "elds of velocity, temperature and composition within the #ame, is developed. The dynamics of particle growth is simulated by integrating di!erential equations for the evolution of primary particle size and number concentration of aggregates along trajectories following the particle #ow.

* k oe# ' (oue)" ' R(e) *t p C k e k u:u#(1!C ) ! R go #C C Ck R op F e !C o C k

3.1. Flame model

where k , the turbulent viscosity is given by R k "oC k/e (6) R I and the coe$cients p , p , C , C and C have been asC I C C I signed the following empirical values (Launder & Spalding, 1972):

The #ame model comprises the time-smoothed equations of change together with models for the turbulent #uxes. Turbulence must be included since eddy di!usion a!ects the radial transport rates in particular. The model is set up and solved for the whole gas-phase domain contained within the protective shield surrounding the #ame (cf. Figs. 1 and 3) using the computational #uid dynamics solver,Fluent (Fluent Inc., 1996). The principal governing equations, which form the basis for the calculations are shown below. For further details, readers are referred to Bird, Stewart and Lightfoot (1960) and Fletcher (1997) and the software manuals (Fluent Inc., 1996. The equation of motion becomes (Bird et al., 1960): * (ou)"![ ' ouu]! p![ ' s]![ ' sR]#og. *t (2) In this notation, scalars are lightface roman, vectors are boldface roman and tensors are boldface Greek. The







 (5)

C "1.44, C "1.92, C "0.09, C C I p "1.3, p "1.0. (7) C I As boundary conditions for the equation of motion (2) the spatial variation of u and all Reynolds stresses are speci"ed at all inlet boundaries. Along the outlet boundary surfaces all normal gradients vanish, and an even distribution of pressure is assumed. There is no-slip #ow conditions along the walls, and in the near wall region turbulence is controlled by nonequilibrium wall functions (Launder & Spalding, 1974). The continuity equation for component i yields * (om )# ' om u"! ' (j #jR)#S G G G G G *t

(8)

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

181

where m is the mass fraction of species i, S the net rate of G G production of species i per unit volume by chemical reactions. The molecular di!usive mass #ux j vector G becomes

¹ is a reference temperature and c is the speci"c heat  NG capacity of the species i. The molecular and turbulent di!usion #uxes of heat are given by

j "!oD m (9) G GK G where D is the molecular di!usion coe$cient for speGK cies i in the gas mixture. The turbulent mass #ux vector jR of species i in terms of the turbulent viscosity is given by G jR"!(k /p )m (10) G R K G with the turbulent Schmidt number p "0.7. The rate of K methane combustion is calculated as the limiting value of a regular Arrhenius expression for the single step methane combustion (Bartok & Saro"m, 1991) and the rate calculated by the Magnussen and Hjertager (1976) turbulent combustion model. The list of chemical species implied by Eq. (8), in addition to CH , O , CO , H O, and N , includes a 6th      pseudo-component PM, which describes the transport of particulate matter within the #ame. In the CFD simulations, the PM component does not describe changes in number concentration of particles as a result of coagulation but only the mass fraction of particulate matter at any position in the computational domain. The in#uence of this dilution/dispersion on the coagulation rate is discussed below. Due to their high dilution in the gas, the in#uence of the particles on the #ow and temperature "elds can be neglected. Within the #ame, particles are transported solely by gas convection and turbulent di!usion, because for particles of all relevant sizes, the molecular, sizedependent di!usion coe$cient (Eq. (9)) is negligible compared to the eddy di!usion coe$cient (cf. Eq. (10)), which is independent of size. Molecular di!usion will only be signi"cant for particles consisting of relatively few monomers, which is only is a matter of microseconds. Inertial forces can also be neglected for the relevant particle sizes. Hence one can expediently and accurately simulate the transport of particles, irrespective of size, by that of a gas component PM with a negligibly small molecular di!usion coe$cient. The magnitude of the remaining physical properties of PM has virtually no in#uence on the #ame model, if they are chosen within the ranges of the other gas species. The equation of energy in terms of the speci"c enthalpy h becomes

(13)

* *p (oh)# ' ohu"! ' (q#qR)! ' h j # H H *t *t H # u ' p!s:u!oe#S F with



2 h" m h where h c d¹ G G G NG  2 H

kc q"!k ¹, qR"! R N¹ A p F

where k is the thermal conductivity and p "0.7 is the A F turbulent Prandtl number. In Eq. (11) the source term S denotes the rate of heat evolvement from chemical F reactions and radiant heat absorption. The second term on the right-hand side denotes the transport of enthalpy by species di!usion. The irreversible rate of increase of the enthalpy by viscous dissipation is partitioned into its laminar part, the scalar term s:u and its turbulent part oe. Radiant energy absorption and irradiation is incorporated in the enthalpy balance by the so-called P-1 model (Siegel & Howel, 1992), in which the radiation #ux is expressed by expanding the radiation intensity in series of orthogonal harmonics. The complete system of equations is solved using a "nite volume technique (Fluent Inc.), based on the partitioning of the physical domain into discrete control volumes using a computational grid. The governing equations are integrated on the individual control volumes and combined with the boundary conditions to construct algebraic "nite di!erence equations for the discrete unknowns: velocity, pressure, species mass fractions, etc., for each control volume. The set of discretized equations is solved by an algebraic equation solver. The burner and the shield surrounding it has complete rotational symmetry around the central jet axis which

(11)

(12)

Fig. 3. The computational grid for the #ame environment. Part of the grid above the burner mouth and at the outlet is omitted. There is complete rotational symmetry around the axis to the right. The grid is re"ned in areas with large gradients.

182

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

minimizes any swirling #ow. Therefore, an obvious choice is to carry out the simulations with a two-dimensional computational grid with a rotational symmetry axis (Fig. 3). The grid is structured with 65 nodes in the radial direction and 110 nodes axially. The solutions are invariant to a doubling of the grid density in axial and radial direction simultaneously, and the grid resolution is concluded to be beyond numerical in#uence. The computed discrete values of the dependent variables in the grid nodes are stored in arrays for the subsequent analysis of particle dynamics. The velocity components, temperature, etc. along given trajectories in the computational domain are computed by an e$cient piecewise polynomial interpolation function for two-dimensional gridded data (Boisvert, 1982). 3.2. The dynamics of particle growth The decomposition of the precursor is assumed to occur instantaneously when the gas temperature exceeds a given value ¹ , converting all precursor molecules to " free Al O `monomera molecules. The subsequent   formation and growth of particles occur by coagulation caused by Brownian collisions between particles with the monomers as starting point and leading to a steadily growing average particle size. The notion of a collision limited growth, even for the monomers, is based on the assumption that dimeric and larger clusters are stable * a plausible assumption due to the extreme supersaturation of the monomers at the time of their formation. Our model for the rate of change of the number concentration of particles is based on the model of Kruis et al. (1993), which assumes a monodisperse size distribution during the coagulation process. Although obviously approximate, this model re#ects the fact that the size distribution of a coagulating aerosol eventually attains a so-called self-preserving form (Friedlander, 1977) with a constant geometric variance (Matsoukas & Friedlander, 1991). By comparison with detailed sectional models of aggregate dynamics (Xiong & Pratsinis, 1993), the model by Kruis et al. (1993) was shown to predict the average primary particle size well. The model is modi"ed to include the e!ect of dilution caused by radial eddy dispersion in the #ame (Johannessen, 1999). The dilution factor f is calculated by " f "m /m " .+ .+

(14)

where m is the mass fraction of the PM pseudo gas .+ component and subscript 0 denotes the conditions at the start of the precursor jet. It is assumed that ensembles of particles move together along well-de"ned trajectories, i.e. it is assumed that dilution is the only e!ect of radial dispersion. Then the rate of change of the particle number concentration and the corresponding growth rate for the volume of one

particle, v, are given by (Johannessen, 1999): dN 1 Ndf " "! o bN# dt 2 E f dt " and




(15)



dv 1 1 1 " o bN v" o bNv dt N 2 E 2 E

(16)

where N is the number of particles per unit mass of gas, o is the gas density, b is the collision frequency function E and t is the elapsed time along a particle trajectory. The second term on the right-hand side of Eq. (15) denotes the rate of concentration change due to dilution of the central jet by the surrounding particle free gas. The use of a concentration based on mass of gas is expedient because it is invariant to the thermal expansion of the gas. The collision frequency b is calculated by Fuchs equation (Seinfeld, 1986; Fuchs, 1964):






r (2D A # . b"8pDr A 2r #(2g cr A A

(17)

The particle di!usion coe$cient, D, particle velocity, c, and transition parameter, g, are given by e.g. Seinfeld (1986) and Kruis et al. (1993). These relations are correct for collisions in the continuum regime, for Kn1, and in the molecular regime, for Kn1, and provide a semi-empirical interpolation in the intermediate regime. The collision radius, r , is equal A to the true radius for a solid sphere. However, when the particles are shaped like dendritic aggregates with fractal-like structure, their collision radius depends on the fractal dimension, primary particle size and number of particles in the aggregate. The collision radius for an aggregate consisting of n identical primary particles is N calculated by (Matsoukas & Friedlander, 1991):




v "D r "r "r (n )"D. (18) A N v N N N The fractal dimension is assumed to be D "1.8, a value D often observed for cluster}cluster coagulation (Kaplan & Gentry, 1988). A solid sphere has a fractal dimension equal to three. The net rate of change for the surface area of an aggregate, a, consists of two terms, the rate of area growth due to coagulation and the rate by which the surface area becomes diminished by surface tension forces, which seeks to minimize the surface area * the process usually called sintering or coalescence:






da 1 1 " o bN a!r (a, v, ¹). Q dt N 2 E

(19)

The exploration of sintering rates for ultra "ne particles is a quite new branch of science and the theory for the relationships imbedded in r (a, v, ¹) is as yet poorly Q

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

developed, although it is extremely important for the control of particle morphology, speci"c surface area, etc. For the "nal stages of sintering of a slightly non-spherical particle, initially consisting of two contacting spheres with diameter d , the following expression applies (Koch N & Friedlander, 1990; Wu & Friedlander, 1994): 1 r (a, v, ¹)" (a!a ) (20) Q Q q D where a is the surface area, v the volume of the particle and a is the surface area of a solid sphere with the same Q volume as the sintering particle. The characteristic coalescence time, q , depends on the surface tension, on a, v D and ¹ and on the sintering mechanism, which may be any combination of viscous #ow, lattice di!usion, surface di!usion or grain boundary di!usion (cf. e.g. Lehtinen, 1997). For each of these mechanisms Table 2 shows expressions for q , all yielding to the common form: D q +dK/C(¹). (21) D N Below the melting point of the particles, the three di!usion mechanisms of Table 2, i.e. lattice-, surface-, and grain boundary di!usion presumably are the most important mechanisms. They are thermally activated processes with the di!usion coe$cient following an Arrhenius expression (Kobata, Kusakabe & Morooka, 1991). Consequently, the general expression for q of Eqs. D (20) and (21) for the three di!usion mechanisms is







¹ E 1 1 q "k dK exp  ! (22) D  N¹ R ¹ ¹   where E is the activation energy of the di!usion coe$c ient and d is the initial diameter of the two contacting N Table 2 Expressions for the characteristic coalescence time q of Eq. (21) for D di!erent coalescence mechanisms (Lehtinen, 1997) Coalescence mechanism

C(¹)

m

Viscous #ow Evaporation condensation

3p 2k 6ppM (2no(R¹) N 64D pM J o R¹ N 225dD pM Q o R¹ N 192wD pM @ o R¹ N

1 2

Lattice di!usion Surface di!usion Grain boundary di!usion

3 4 4

p"surface tension, R"gas constant, o"particle density, k"viscosity, ¹"temperature, D"di!usion coe$cient, p"pressure, d"surface layer thickness (l"lattice, s"surface and b"grain boundary), M"molar mass, w"grain boundary width.

183

primary particles. The pre-exponential term, k dK, is the  N characteristic sintering time at the temperature ¹ .  By the generation of particles in a #ame, dendritic aggregates consisting of numerous, easily discernable primary particles often occur at an intermediate stage and/or as the end product. The sintering kinetics of Eqs. (20)}(22) in a strict sense is limited to fairly compact particles, i.e. particles only moderately perturbed from a spherical shape. Several attempts have been made to extend the kinetics for more complicated shapes. Either by a direct application of Eq. (20) to particles of any shape (Friedlander & Wu, 1994; Kruis et al., 1993), in the latter study with q based on the size of the primary D particles. A newer approach is to partition an aggregate particle into groups of primary particles which sinter independently (Lehtinen, Windeler & Friedlander, 1996). The predictions of these di!erent models di!er substantially (Johannessen, 1998). We propose a new model for the sintering kinetics of particles/aggregates, which elaborates the ideas of Ulrich and Subramanian (1977) and Lehtinen et al. (1996). When the particles are relatively compact, one can still use Eqs. (20)}(22). For aggregates of any size (more than two primary particle), we assume that they consist of equally sized, spherical primary particles. The local characteristic sintering time is computed from the size of the primary particles and the total sintering rate is related to the number of contact points. For two coagulated particles, the area reduction for a "nite period primarily occurs in the region close to the point of contact between the particles. The rate of this reduction for one contact point is assumed to follow Eqs. (20)}(22) with d equal to N the primary particle size. The net rate of reduction for an aggregate is obtained by multiplication with the number of contact points in the aggregate. The idea behind the points-of-contact concept for aggregate sintering, visualized on Fig. 4, implies that the coalescence within an aggregate not only increases the primary particle size but also reduces the number of contact points. The sintering rate expression becomes








(a!a ) 6(v/2)  Q ; dH" N q (dH) p D N

r " (n !1) Q N

for n )2, N

(2a !a ) N NQ q (d ) D N

0.41a N "(n !1) N q (d ) D N

(23)

for n '2. N

The diameter and number of primary particles in an aggregate, d and n , and a , the surface area of a primary N N N particle are computed from the particle volume v and surface area a by simple geometrical considerations,

184

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

sintering rate r of Eq. (23) and the coalescence time of Q Eq. (22) inserted in Eq. (19). The integration occurs along particle trajectories in the r,z-space. Assume that the trajectory is given by a parametric form as s(t)"[r(t),z(t)] Fig. 4. The POC (Point Of Contact) approach to the coalescence of aggregates. The local coalescence rate depends on the primary particle size according to Eq. (23). The overall coalescence rate is obtained by multiplying the local rate by the number of connecting points.

together with v the volume of one primary particle: N 6v p v d " , v " d, n " , a "pd. N N 6 N N v N N a N

(24)

In Eq. (23), q (d ), computed by, e.g. Eq. (22), is the D N characteristic coalescence time for two contacting spherical particles with diameter d . N The derivation of Eq. (23) follows directly from the model of Fig. 4: When n (2, the rate of sintering is N proportional to the excess surface, a!a , and follows the Q &compact particle' kinetics of Eqs. (20)}(22) with a sintering diameter dH, de"ned as the diameter of a sphere with N volume v/2. When n '2, the sintering rate at each point of contact N is proportional to an excess surface of 2a !a with N NQ "2a the surface area sintering diameter d and a N NQ N of a sphere with volume 2v . The number of contact N points is assumed equal to n !1, corresponding to N a completely branched aggregate without cyclic structures. In the points-of-contact (POC) model, the coalescence of an aggregate is completely described by the expression for the coalescence of two primary particles and the number of contact points in the aggregate. There are no restrictions on the form of the expression for the characteristic sintering time. The POC-model ensures a correct asymptotic behaviour for the "nal stage of coalescence when the whole particle is compact. The transition from aggregate sintering to compact particle sintering occurs without discontinuity in the rate of sintering. Based on the limited experimental data below, we cannot with signi"cance discriminate between the three alternative models for r , i.e. the POC-model and the Q models of Kruis et al. (1993) and Lehtinen et al. (1996). Presently, we do prefer the POC-model for two intuitively attractive qualities, i.e. the sintering time constant q refers to the primary particle size, and the relative rate D of sintering r /a is independent on n, the number of Q primary particles in an aggregate, when n is large. 3.3. Particle trajectories The three di!erential Eqs. (15), (16) and (19), with N,v and a as dependent variables are integrated with the

(25)

where t is the time elapsed for a particle moving along the trajectory. This trajectory actually de"nes an envelope surface surrounding the #ame axis due to the cylindrical symmetry of the #ame. Any given trajectory is de"ned such that it envelops a constant fraction a of the total mass #ow of particles w (kg/s) or for any z: .+



2p

PY o m u r dr"aw for z"z E .+ X .+ 

(26)

where a is constant for a given trajectory and m is the .+ mass fraction of particulate matter and u is the axial X component of the mean velocity. The total particle mass #ow rate is known from the precursor #ow of the jet and is given by



w "2p .+

P  

o m u r dr. E .+ X

(27)

The lines de"ned by Eq. (26) are trajectories only in the sense that they depict the most likely subsequent paths of particles positioned anywhere at a trajectory. Due to the random nature of turbulent di!usion, individual particles may deviate from this path, which hence is an approximation to the real particle paths, strictly correct only for purely laminar #ow. The axis of time t in the three di!erential Eqs. (15), (16) and (19) is established by the requirement that the tangential particle velocity along the trajectory computed from Eq. (25) as ds/dt is equal to the tangential component of the mean gas velocity u in any point of the trajectory. For any value of t, the gas temperature and composition, required for the right-hand sides of the three di!erential equations, is calculated via Eq. (25) by interpolation between the grid values of the solution matrices. The initial conditions for N,v and a are given by the assumption of an instantaneous combustion of the precursor at the decomposition temperature ¹ , at which " point N and v , the concentration and volume, respec  tively, of the initially formed single monomer molecules are calculated from the precursor content in the gas. The initial condition for a is 6v /d , where d is the diameter  N N of the monomer molecule. In both cases the cup-mix average values of N,v and a #owing through a horizontal section of the #ame is given by integration over all trajectories at constant z: P No u r dr E X NM "  P o u r dr  E X

(28)

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

with similar expressions for a and v . The advantage of the trajectory method is its very fast execution time on the computer. It allows rapid iterations on the sintering kinetics, even for several di!erent #ames simultaneously, because the #uid #ow and temperature "elds need not be recalculated between iteration steps when the in#uence of the particles on these "elds is negligible. The error introduced by including only the e!ect of dilution * not of particle mixing * as a result of turbulence, is insigni"cant as long as the variation of particle characteristics between trajectories is only as moderate as for the present study. The implementation of the particle dynamics equations directly into the CFD code avoids this error, but is rendered di$cult due to numerical di$culties from the highly non-linear source terms (coagulation and sintering), and yields a much slower code. Results from work with this type of model will be presented in an ensuing article.

4. Results and discussion Fig. 5 shows the calculated #ame temperature distribution for four runs from Table 1, selected to illustrate the variation of the #ame structure during the experiments. The #ame photographs inserted on Fig. 5 show that the calculations qualitatively conform with the visual ap-

185

pearance of the #ame with respect to size and shape including the di!erence between the single di!usion #ame C and the double di!usion #ame B, for which the cold core region is broader. A high #ow rate of methane (B1 and C1) obviously yield a larger #ame region than a low rate (B3 and C3). The double di!usion #ames (B1 and B3) have methane surrounded by air on both sides and hence yield a more e$cient combustion than the single di!usion #ame (C1 and C3), and they have a somewhat more contracted #ame region and a slightly higher peak temperature. For B3 the temperature peak actually occurs away from the axis and removed from the precursor jet. Measured and calculated #ame temperatures are compared on Figs. 6 and 7. Fig. 6 shows the central axis temperature and Fig. 7 the maximum temperature at given heights above the burner mouth. Due to the di$culty of measuring precisely with a thermocouple in the extreme environment of the #ame some deviation between measurement and calculation is to be expected. This is particularly true for the centreline temperature, where the target point is surrounded by steep temperature gradients during the initial heating up of the jet. The measurement of the maximum temperature, which is made by small radial displacements of the thermocouple to obtain the maximum reading, is more robust, which is also evident from the agreement with the theoretical

Fig. 5. Photographs of the Al O -#ames compared with temperature contour plots from the CFD-simulations. The structural di!erence between   B-#ames (double di!usion #ame) and C-#ames (single di!usion #ame) is clearly seen on both the simulation and the photographs. The terms B1-3 and C1-3 refer to the operating conditions in Table 1. Two particle trajectories, corresponding to a"0.1 (trajectory 1) and a"0.9 (trajectory 5) respectively (cf. Eq. (26)), are superimposed on each contour plots. Each trajectory has a unique time}temperature}dilution history.

186

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

Fig. 6. Measured and simulated centre axis temperature pro"les for #ames B4 and C4. Operating conditions as in Table 1.

Fig. 7. Measured and simulated maximum temperatures as a function of height above burner mouth for #ames B4 and C4. The plateau at high temperature for #ame C4 is evidently well reproduced by the simulations. The curve for #ame B4 has a single peak. Operating conditions as in Table 1.

curves for the two types of measurements. We consider the agreement between measurements and calculation in Figs. 6 and 7 a satisfactory validation of the simulated temperature "eld. For the simulation of particle dynamics "ve di!erent trajectories are selected. They are de"ned by Eq. (26) with a"0.1,0.3,0.5,0.7, and 0.9. The innermost and outermost trajectories are shown for each of the #ame "elds of Fig. 5. Both convection and eddy di!usion contributes to the radial spreading of the particle jet, however the latter mechanism is by far the most important. Radial convection of particles is caused by the radial component of the mean gas velocity u, and if this was the only mechanism the particles would follow streamlines s (cf. Eq. (25)) de"ned by ds/dt"u.

(29)

A streamline thus de"ned with the same starting point as the outermost trajectory of Figs. 5 and 6 is removed 0.5}1 mm from the axis at the top of the computational domain as opposed to 10}14 mm for the outermost trajectories de"ned by Eq. (26), thus demonstrating the importance of radial eddy dispersion. The in#uence of temperature especially a!ects the sintering rate, while the dilution of the particles a!ects the coagulation rate. Hence, the variation of temperature and dilution factor f with time along the trajectories are " the dominating factors for the eventual shape and size of the aerosol particles. Theses relations depend on the #ame con"guration and may vary from one trajectory to another. Figs. 8 and 9 depict the dilution factor and temperature for selected operation conditions. The rate of dilution of the particles in the central jet is evidently fastest along the outermost trajectories which creates a signi"cant variation of f between di!erent trajectories " for intermediate values of t, until the particle concentration eventually attains an almost even radial distribution at the outlet of the computational domain (Fig. 8). The rate of dilution is faster for #ame con"guration B than for C due to the in#uence of the high air #ow rate adjacent to the central jet for the former con"guration (Fig. 8). Fig. 9 reveals the important di!erence between the double di!usion #ame B and the single di!usion #ame C. The peak temperature for #ame B varies considerably, i.e. from 1250 to 1950 K for the core trajectory when the methane #ow rate is changed from 0.38 to 0.145 l/min. The same change in methane #ow rate a!ects the peak temperature for #ame C far less, merely shifting the position of the peak towards lower values of t. Furthermore, the di!erence between the ¹,t curves for di!erent trajectories is more pronounced for the B than for the C #ames. The only unknown parameters in the model for the particle dynamics are the sintering rate constant k and  activation energy E of Eq. (22). These two values are 

Fig. 8. The dilution factor f of Eq. (14) for the alumina #ames B3 " and C1. The turbulence created by the variation of the linear gas velocities at the burner mouths a!ects the radial dispersion rate and hence f . "

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

187

Fig. 10. Measured values (points) and simulated values (curves) of the speci"c surface area of the Al O -particles. The operating conditions   are listed in Table 1.

Fig. 9. The variation of temperature along "ve trajectories for four alumina #ame con"gurations. The "ve trajectories correspond to a"0.1,0.3,0.5,0.7, and 0.9 respectively.

estimated by "tting the model to the measured particle surface areas with a least squares routine minimizing SSQ" (ln(S )!ln(S ))      G G

(30)

where i denotes the twelve measurements of Table 1. The value of S is calculated as the mean value of a/(o v)   N at the outlet of the #ame shield. The mean value is obtained from Eq. (28) after integration of the particle dynamics model (Eqs. (15), (16) and (19)) along "ve trajectories with a"0.1,0.3,0.5,0.7, and 0.9, respectively, for each of the twelve measurements. The model is solved iteratively under the control of a modi"ed Marquard routine (Fletcher, 1971). The derivatives of SSQ with respect to the estimated parameters, required in each iteration step by the Marquard code, are computed by numerical perturbation of the particle dynamics model. The initial conditions for the particle volume, surface area and number concentration are calculated from the assumption that all the solid oxide is instantly released as monomers, considered to be single, spherical Al O mol  ecules at a given decomposition temperature ¹ "470 " K, approximately 20 K above the inlet temperature. The value of ¹ a!ects the "nal particle characteristics insig" ni"cantly as discussed below. However the chosen value minimizes SSQ of Eq. (30). Di!erent values of m in Eq. (22) were attempted. A value of m"4 yields a slightly better model "t than m"3, which points to surface or grain boundary di!usion as the most likely coalescence mechanisms. However, the di!erence is hardly signi"cant enough to discriminate with certainty between the two values of m. Thus the parameters in Eq. (22) were estimated as: k "(1.3$0.6);10 s/m and E "68$24  

kJ/mol. These results are based on the following values of the constants in Eq. (22): m"4 and ¹ "1400 K. The  estimated two-p con"dence limits are shown together with the parameter estimates. The estimated standard error in predicting S by the model is 10%.  Fig. 10 compares the measured and computed speci"c surface areas. There is a rather moderate variation of the surface area for the twelve alumina measurements, i.e. from 48 to 94 m/g. This is in accordance with the small sintering activation energy estimated for alumina. Windeler, Friedlander & Lehtinen (1997) studied coagulation-coalescence kinetics for the #ame generation of oxide particles, including TiO , Nb O and Al O . The      ultimate primary particle size varies inversely with the bulk solid di!usion coe$cient indicating that bulk solid di!usion is important for the coalescence mechanism. However, they conclude that Al O probably coalesces   by a di!erent mechanism than the other materials. Our results supports this, since the activation energy of the bulk di!usion coe$cient of alumina (Windeler et al., 1997) far exceeds the coalescence rate activation energy of 68 kJ/mol found in our study. We can further validate the model by comparison with the measurements by transmission electron microscopy (TEM) and di!erential mobility analyser (DMA). The TEM photos (Fig. 11) provide a clear depiction of the dendritic structure of the aggregate particles. The aggregate size cannot be determined by TEM, because, "rstly, the individual aggregates cannot be identi"ed due to overlap with other aggregates and, secondly, the deposited particles are not necessarily representative for those in the gas with respect to aggregate size. However, the primary particles are clearly discernible and their size appears uncorrelated with the aggregate size. The primary particle size is evidently more homogeneous for the single di!usion #ame C, than for the double di!usion #ame B, in accordance with the more homogeneous C-trajectory curves (cf. Figs. 8 and 9). The average primary particle size and its standard deviation for 100 randomly selected primary particles on the TEM photos

188

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

Fig. 11. TEM images of alumina particles from di!erent #ame operating conditions. Large dendritic aggregates of primary particles are seen for all conditions. The size distribution of the primary particle size seems quite narrow for #ames C1 and C3, whereas #ames B1 and especially B3 have a wide range of primary particle sizes. For three of the four images, the di!raction images have been inserted. The spots and formation of rings show that at least some of the particles are crystalline.

was measured for all C-#ame alumina series and are shown in Table 3 together with the computed by Eq. (1). For the B-#ames the primary particle size cannot be measured reliably because they are too inhomogeneous. The &BET'-equivalent primary particle size in Table 3 varies with the methane and nitrogen #ow in qualitatively the same way as the TEM-based size, but is constantly, approximately 30% larger. This discrepancy is attributed to the necking phenomenon, i.e. the vanishing of surface area around the point of contact between two partly fused particles and therefore is a natural di!erence between the two di!erent de"nitions of the primary particle size. The logarithmic size distributions of aggregates measured by the di!erential mobility analyser (DMA) are

normally mono-modal, quasi-Gaussian shaped curves (Johannessen, 1999). The geometric standard deviation for the distributions varies in the range 1.8}2.4. The size measured by the DMA is the mobility diameter, de"ned as the diameter of a sphere with the same mobility as the particle. The average mobility diameter for a given distribution is closely related to the collision diameter 2r A * the model variable computed by Eq. (18), although the exact relation between the two di!erently de"ned diameters are not yet fully understood for dendritic aggregates. The measured mean mobility diameter and the

 The mobility is de"ned as the steady state velocity of a particle a!ected by a unit force.

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191

189

Table 3 The primary Al O -particle size for the C-#ames measured from TEM images and computed from BET surface area data   Flame

C1

C2

C3

C4

C5

C6

d (TEM) N d (BET) N

18.1$1.2 29.5

16.9$1.8 24.7

15.8$0.9 22.1

20.9$2.5 32.7

20.2$2.3 27.5

16.5$2.5 23.6

All numbers are in units of nm, and the TEM-results are shown with the estimated standard deviation.

Fig. 13. The variation of temperature along the "ve trajectories for #ame B3. Operating conditions of Table 1. Fig. 12. Measured and simulated Al O -aggregate sizes. The experi  mental values are the mean mobility diameters of the size distributions measured by the DMA. The simulated values are based on the aggregate collision diameter in the simulation model. Operating conditions of Table 1.

computed collision diameter for the twelve alumina series are compared in Fig. 12. The model, in accordance with the measurements, predicts a fairly small variation of aggregate size for the entire series, and there is an interestingly "ne agreement between collision diameter and measured mobility diameter, in particular for the B#ames, where the detailed variation of size with the imposed experimental conditions are followed very neatly by the simulations. The agreement between measurements and simulations in Figs. 10 and 12 and Table 3 is quite gratifying, and indeed proves that the model in great detail re#ects the experimental observations. The ratio between the rates of coalescence and coagulation varies signi"cantly along a trajectory, and the model simulation reveals intricate details about the resulting variation of particle morphology and size (Figs. 13}16). The number of primary particles per aggregate, n , increases when this ratio is large, which is the case for N the &cold' period during the initial heating up and also for the "nal stage of the #ame. During the intermediate period, where the temperature is high, the rate of coalescence dominates and n decreases. Accordingly, n typiN N cally shows both a maximum and a minimum at intermediate positions (Fig. 14), which causes the colli-

sion diameter to exhibit a similar behaviour (Fig. 16). The size of the primary particles grows monotonically with a growth rate which has a maximum at some intermediate position and becomes negligible below approximately 1000 K (Fig. 15). Figs. 13}16 show the importance of including radial dispersion of particulate matter in the simulations. In systems with large radial gradients, one cannot assume that all particles will experience identical time}temperature}dilution histories. The homogeneity of the process outlet depends highly on #ame con"guration and material properties. The value of the decomposition temperature ¹ and " the fractal dimension D are both subject to some uncerD tainty. Therefore, the sensitivity of the model to changes in these parameters has been studied. The sensitivity of a calculated model variable towards changes in a parameter m is de"ned as Sensitivity"* ln /* ln m.

(31)

Hence, a sensitivity equal to one implies that a fractional change in m leads to an identical fractional change in . The sensitivity of the speci"c surface area and the aggregate size to ¹ and D are summarized in Table 4. Cha" D nges in D primarily a!ect the aggregate size but is D virtually without in#uence on the speci"c surface area. Hence, the establishment of a proper value of D can be linked with measurements of the aggregate D size, which for the present investigation show a satisfactory agreement with the model (cf. Fig. 12). The

190

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191 Table 4 Sensitivity analysis: The e!ect of the decomposition temperature, T , " and fractal dimension, D , on the outlet speci"c surface area and the D aggregate size Speci"c surface area(S ) 

Flame Flame Flame Flame Fig. 14. The variation of the number of primary particles in one aggregate along the "ve trajectories for #ame B3. Operating conditions of Table 1.

C1 C3 B1 B3

Aggregate collision diameter (d ) A

¹ "

D D

¹ "

D D

0.06 0.02 0.64 0.36

0.02 0.01 !0.07 0.05

!0.15 0.03 !0.12 !0.04

!1.77 !1.59 !1.26 !2.16

The individual elements are computed by Eq. (31).

5. Summary

Fig. 15. The variation of the primary particle diameter along the "ve trajectories for #ame B3. Operating conditions of Table 1.

Fig. 16. The variation of the aggregate collision diameter along the "ve trajectories for #ame B3. Operating conditions of Table 1.

The combination of a simple batch model for coagulation and coalescence of aerosol particles with a computational #uid dynamics model of temperature, velocity and gas composition in a di!usion #ame has provided useful insight into the fundamental processes during #ame synthesis of alumina particles. The model requires the values of the two parameters for particle sintering. These are not generally available, but by "tting the model to experimental data, an empirical sintering rate for alumina, including an e!ective activation energy was obtained. The resulting model is capable of predicting the speci"c surface area, the primary particle size and the average mobility diameter of the generated particles with satisfactory accuracy for a range of conditions. The simulation model therefore should be useful as a tool for the optimization of #ame processes for particle production. In this study, the thermophoretic forces on particles and the in#uence of particles on radiative heat transfer in the gas phase has been neglected. These phenomena may have some e!ects on the particle dynamics and needs to be further studied. There is a need for further study of the coalescence and recrystallization kinetics of ultra "ne particles. The simple kinetics used here probably only apply for a limited range of temperature and particle size.

Acknowledgements This work was supported by a grant from the Danish Council for Technological Research.

decomposition temperature ¹ only a!ects the speci"c " surface area of the double di!usion #ames and only to a very moderate degree. It has negligible e!ect on the surface area for the single di!usion #ames and on the aggregate size. Evidently, the model is not signi"cantly a!ected by the uncertainty of ¹ and D . " D

References Bartok, W., & Saro"m, A. F. (1991). Fossil fuel combustion. New York: Wiley.

T. Johannessen et al. / Chemical Engineering Science 55 (2000) 177}191 Bird, R. B., Stewart, W. E., & Lightfoot, E. N. (1960). Transport phenomena. New York: Wiley. Boisvert, R. (1982). Fortran Subroutine B2INK. Scienti"c Computing Division, National Bureau of Standards, Washington, DC 20234. Collis, D., & Williams, M. (1959). Two-dimensional convection from heated wires at low Reynolds numbers. Journal of Fluid Mechanics, 6, 357. Fletcher, C. A. J. (1997). Computational techniques for yuid dynamics (2nd ed.). New York: Springer. Fletcher, R. (1971). A modixed Marquard subroutine for non-linear least squares. Harwell Report AERE.6799. A computer code, VA07, is included in the Harwell subroutine library. Fluent, Inc. (1996). Fluent: User's Guide, 1}4, Release 4.4, Fluent Incorporated, Centerra Resource Park, 10 Cavendish Court, Lebanon, NH 03766. Friedlander, S. K. (1977). Smoke, dust and haze. New York: Wiley. Friedlander, S. K., & Wu, M. K. (1994). Linear rate law for the decay of excess surface area of a coalescing solid particle. Physics Reviews B, 49, 347. Fuchs, N. A. (1964). Mechanics of aerosols. New York: Pergamon Gurav, A., Kodas, T., Pluym, T., & Xiong, Y. (1993). Aerosol processing of materials. Aerosol Science and Technology, 19, 411}452. Haanappel, H. D., van Corbach, H. D., Fransen, T., & Gellings, P. J. (1994). The Pyrrolytic decomposition of aluminum-tri-sec-butoxide during chemical vapour deposition of thin alumina "lms. Thermochimica Acta, 240, 66}77. Johannessen, T. (1999). Synthesis of nano-particles in yames. Ph.D Thesis, Department of Chemical Engineering, Technical University of Denmark, Denmark. Kaplan, C. R., & Gentry, J. M. (1988). Agglomeration of chain-like combustion aerosols due to Brownian motion. Aerosol Science and Technololgy, 8, 11}28. Kobata, A., Kusakabe, K., & Morooka, S. (1991). Growth and Transformation of TiO crystallites in aerosol reactor.. A.I.Ch.E. Journal,  37, 347}359. Koch, W., & Friedlander, S. K. (1990). The e!ect of particle coalescence on the surface area of a coagulating aerosol. Journal of Colloid and Interface Science, 140, 419. Kruis, F. E., Kusters, K. A., Pratsinis, S. E., & Scarlett, B. (1993). A simple model for the characteristics of aggregate particles undergoing coagulation and sintering. Aerosol Science and Technology, 19, 514. Launder, B. E., & Spalding, D. B. (1972). Lectures in mathematical models of turbulence, London: Academic Press.

191

Launder, B. E., & Spalding, D. B. (1974). The numerical computation of turbulent #ows. Computer Methods in Applied Mechanics and Engineering, 3, 269}389. Launder, B. E., Reece, G. J., & Rodi, W. (1975). Progress in the development of a Reynolds-stress turbulent closure. Journal of Fluid Mechanics, 6, 537}566. Launder, B. E. (1989). Second-moment closure: Present 2 and future?. International Journal of Heat Fluid Flows, 10, 282}300. Lehtinen, K. E. J., Windeler, R., & Friedlander, S. K. (1996). A note on the growth of primary particles in agglomerate atructures by coalescence. Journal of Collied and Interface Science, 182, 606}608. Lehtinen, K. E. (1997). Theoretical studies on aerosol agglomeration processes. Ph.D thesis, VTT Technical Research Centre of Finland. Magnussen, B. F., & Hjertager, B. H. (1976). On mathematical model of turbulent combustion with special emphasis on soot formation and combustion. 16th Symposium (International) in combustion, Cambridge, MA, August 15}20. Matsoukas, T., & Friedlander, S. K. (1991). Dynamics of aerosol agglomerate formation. Journal of Collied and Interface Science, 146(2), 495}506. Pratsinis, S. E. (1998). Flame Aerosol Synthesis of Ceramic Powders. Progress in Energy and Combustion Science, 24, 197}219. Pratsinis, S. E., Zhu, W., & Vemury, S. (1996). The role of gas mixing in #ame synthesis of titania powders. Powder Technology, 86, 87. Seinfeld, J. H. (1986). Air polution: New York. Siegel, R., & Howel, J. R. (1992). Thermal radiation heat transfer. Washington D.C: Hemisphere. Tranto, J. (1997). Supported ceramic membranes for high temperature gas separation. Ph.D thesis, Dept. of Material Sci., Tech. University of Denmark. Ulrich, G. D., & Subramanian, N. S. (1977). Particle growth in #ames. Combust. Sci. ¹echnol., 17, 119 Vemury, S., & Pratsinis, S. E. (1995). Dopants in #ame synthesis of titania. Journal of the American Ceramic Society, 78, 2984}2992. Windeler, R. S., Friedlander, S. K., & Lehtinen, K. E. J. (1997). Production of nanometer-sized metal oxide particles in a free jet I: Experimental system and results. Aerosol Science and Technology, 27(2), 174}190. Wu, M. K., Windeler, S. E., Steiner, C. K. R., Bors, T., & Friedlander, S. K. (1993). Controlled synthesis of nanosized particles by aerosol processes. Aerosol Science and Technology, 19, 527. Xiong, Y., & Pratsinis, S. E. (1993). Formation of agglomerate particles by coagulation and sintering * Part I. A two.dimensional solution to the population balance equation. Journal of Aerosol Science, 24, 283}300.