Modelling of particulate processes

We review models and numerical methods used in flame synthesis of organic and inorganic nanoparticles. We discuss a general model in which particles form in the gas phase and grow through mass-adding surface reactions, condensation, and coagulation. They shrink or reshape by sintering and mass-abstracting surface reactions. The model is formulated in terms of a population balance which can incorporate a range of levels of detail, i.e. a varying number of internal coordinates. These coordinates can not only describe the geometry of a particle but also its chemical composition or age. In the simplest version a particle is modelled as a sphere whereas in the most complicated form a particle is modelled as an agglomerate of smaller or primary particles where the geometrical shape is known exactly. For these population balance models a number of different numerical approaches exist. We review the method of moments, sectional, finite element, and Monte Carlo methods and give examples of their applications in flame synthesis. Different strategies for coupling a population balance to laminar and turbulent flows are reviewed. For turbulent flows the closure problems arising from chemical reactions and the population balance are briefly discussed. We then summarize the literature on nanoparticle modelling from laboratory to industrial scale and highlight important areas for future research.


Introduction
In this article we review models of combustion synthesis of nanoparticles and the numerical methods used with which to solve them.Nanoparticles can be found in everyday life, for example as pigments, reinforcement material, or even in sunscreen.They are produced on a million ton per year scale and the industrial processes are well established.Also, based on recent research efforts, it is believed that nanoparticles hold great promise for future applications.However, there are also major concerns about the environmental and specifically human health effects of nanoparticle emissions, for example from diesel engines.To improve yields or product quality and to avoid particle emissions, numerical modelling has been extensively used.In this paper we shall discuss models and numerical methods for two types of nanoparticles.
Organic nanoparticles come in two forms namely carbon black and soot.Carbon black is a name for a well-defined industrially manufactured product which is a produced at a rate of several million tons per year.Amongst a variety of applications it is used as filler and reinforcement material with the properties of the individual particles determining their application.There are several production routes which all involve combustion of gases and liquids.The most important is the furnace process in which liquid hydrocarbons are sprayed into the exhaust of a turbulent lean natural gas flame.Soot on the other hand is an unwanted product of combustion, mainly in diesel engines.New emissions regulations make it necessary to reduce the amount of soot formed during combustion drastically.
The variety of inorganic nanoparticles is much greater.Two main routes of synthesis are in use: wet-phase chemistry synthesis and flame synthesis.In many cases flame technology has been identified to be superior as it is a direct, continuous process with little waste and by-product generation [110].For this reason we only examine models which have been developed for flame synthesis.Important inorganic nanoparticles are titania, fumed silica, and alumina powders.However there are many others and recently the number of ceramic powders has grown as a result of new applications, for example as catalysts and surface coatings.These powders are mainly produced similar to carbon black by spraying metal oxide precursors into a high temperature reaction zone.
Two very important aims which one hopes to achieve with the knowledge gained through modelling are the improvement of yield and quality of products and the scale-up of new processes from a laboratory or bench scale level to an industrial level.In both cases the model needs to capture the important physics and chemistry of the formation of nanoparticles, their interactions with each other and with the surrounding gas phase, and their transport through surrounding media.This requires detailed models of the chemical reactions, the population of particles and for the mostly turbulent transport.Despite recent progress in some of these areas the predictive power of the current models remains poor.
In the research community it has been recognized that it is necessary to study sys-tems which are simpler than a spray injected into a reactor of complex geometry and in which the flow is turbulent.For this reason a number of laboratory experiments have been developed, the simplest being a shock tube; other experimental configurations include plug flow reactors, jet stirred reactors, premixed flames, counterflow and axisymmetric co-flow diffusion flames.Whereas for industrial applications only simplified models are available there exist detailed models for the laboratory scale experiments.The findings obtained from these experiments and attempts to model them are summarized in a number of text books and review articles.These are excellent sources of further information on nanoparticles and their applications and a selection is discussed below.The books by Seinfeld and Pandis [105] and Williams and Loyalka [125] are comprehensive treatises on the transport of aerosols and gas phase chemistry in the atmosphere.Although these books do not treat the industrial application of nanoparticles directly they contain a lot of material which is relevant to their modelling.The book by Friedlander [37] is also a classic text and is more relevant to nanoparticles as it contains details on population balance modelling, in particular on agglomerate formation and restructuring.A more general text on population balances is the book of Ramkrishna [95] which describes several numerical techniques for population balance problems each of them relevant for the modelling of nanoparticle population dynamics.Important review articles in the area of soot include the articles by Kennedy [48] and Bockhorn [19] which contain comprehensive collections of references but also describe a number of empirical and detailed soot models.In addition, the book on soot formation edited by Bockhorn [18] is recommended.It contains a collection of papers which address all the important issues concerning the modelling and measurement of soot in a variety of reactors, from lab-scale up to diesel engines.In the area of inorganic nanoparticles, the papers by Pratsinis and co-workers [94,92,66] include a number of interesting references for experimental and modelling work.For educational purposes the article by Rosner [97], the excellent brochure on carbon black from Degussa [74] and the chapter on carbon black in Ullmann's Encyclopedia of Industrial Chemistry [120] are recommended.
The aim of this paper is to provide an up-to-date review of the literature which is concerned with modelling of organic and inorganic nanoparticle population dynamics considering both models and numerical methods.A variety of numerical approaches are discussed.A distinction is made between detailed models describing particle inception, their growth, and their change in structure and models which describe particle transport in laminar and turbulent flows.
The paper is structured as follows.First, a general model for a population of nanoparticles is presented which assumes no spatial gradients in any of the physical quantities.Then different numerical approaches for population balances are reviewed.The next section presents different strategies for coupling the population balance equation with the fluid flow equation in laminar and turbulent flows.A number of applications are presented in the next section.Finally, some new areas for future research are identified.

Outline of a General Model
We first lay out a general model for the population dynamics of nanoparticles as displayed in Figure 1.We assume that all precursors and particles are in a perfectly mixed control volume which means that there are no gradients in any of the physical quantities contained in this volume.The model is defined as follows: • Individual nanoparticles may be completely described by elements of some type space E on which addition corresponding to coagulation is defined.
• The nanoparticle population is described at time t by the number n (t, x) per unit volume of particles of type x ∈ E.
• n evolves according to the discrete Smoluchowski coagulation equation: In these definitions we make the implicit assumption that E is countable so that summations are meaningful.This assumption is common in the literature but not essential and can be removed by replacing the sums with integrals.The (time dependent and non-linear) coagulation operator K is then defined by ( The first sum represents coagulation to form particles of type x and the second loss of particles of type x due to coagulation.K t (x, y) defines a map from the concentrations of particles of types x and y to their coagulation rate at time t given by K t (x, y) n (t, x) n (t, y).K is known as the coagulation kernel.
Surface reactions, condensation, and sintering which only involve one physical particle at a time are described by the linear operator S defined by where l ∈ U. U is an index set for a process or a type of event which is either one element of a set of surface reactions, sintering steps or condensation steps.β (l) t (x) is the rate at which a particle of type x undergoes the change of index l at time t.g (l) (x) is the result of a particle of type x undergoing an event of index l.If a surface reaction removes a particle from the population the function will take the special value 0 in E. We allow for g (l) to be a random function.Note in the deterministic case the probability P in (3) reduces to an indicator function.For example, in the case of x describing the number of monomers in a particle then the addition of a monomer is given by g(y) = y + 1. Models for sintering and surface reactions will be discussed below.
The Particle inception I (t, x) is the rate at which particles of type x enter the system at time t.
As boundary conditions general we will use n (0, x) = n 0 (x) for all x ∈ E .

Representations of nanoparticles
Nanoparticles can be modelled using different degrees of detail which is specified by the set E. In the following, three examples for the set E are given.The models discussed in this paper will either be identical to or are combinations of the following examples.
The Primary Particle Model is the most detailed particle model discussed in the literature.There, a nanoparticle is described as an unordered finite sequence of 'primary particles' and their locations e.g.displacement from the first particle in the list.Primary particles can be characterized as spheres of constant density described by, say, mass but can also contain additional information as for example chemical composition of the surface of a primary particle.The number of internal coordinates is the length of the sequence times the number of particle attributes like distance from the first member in the sequence.Therefore a member of the population of nanoparticles can have several thousand internal coordinates.
The simplest and most widely used model is the Coalescent Sphere Model in which all particles are assumed to be spheres of common density.It is then convenient to describe a particle in terms of the number of monomers.Addition is defined just as for the natural numbers: the result of the coagulation of two spheres is a new sphere with volume equal to the sum of the volumes of the initial spheres, i.e. coagulation is completely coalescent.The internal coordinate of the nanoparticle is its volume or number of contianied monomers.
A slightly more detailed model with two internal coordinates is the Surface and Volume Model.Here we say that a particle is described by volume and surface area.We define addition componentwise, i.e. volume and surface area of the new particle are the sums of the volumes and the surface areas of the interacting particles respectively and hence coagulation is modelled without any coalescence.It is an intermediate stage between the Primary Particle Model and the Coalescent Sphere Model.Processes like sintering, condensation, or surface reactions can transform the volume to surface ratio of an aggregate.The definition of g (l) is non-trivial even with a clear picture of the underlying chemical processes.

Gas phase chemistry
The gas phase chemistry leading to organic soot and Polycyclic Aromatic Hydrocarbons (PAH) has been subject of research for a long time because the combustion of hydrocarbons is so important for many technical and chemical applications.Gas phase chemistry models used for modelling soot and carbon black can be found in [121,3,96].
The gas phase chemistry leading to inorganic nanoparticles is far less well understood and it is very clear that this is an important field of future research.Kinetic data for SiO 2 particles formed from SiH 4 can be found in [60].Reference [113] contains data on the formation of titania particles (TiO 2 ) made from tetraisopropoxide and reference [93] contains data if TiCl 4 is used as precursor.Simple models for the gas phase chemistry describing the formation of Fe 2 O 3 particles from Fe(CO) 5 can be found in [38,77].

Particle inception
Particle inception is probably the most difficult part when modelling nanoparticle synthesis.It is fair to say that it is still not really understood.One approach to model particle inception for inorganic nanoparticles is to apply the theory of homogeneous nucleation of Vollmer-Becker-Döring-Zeldovich which dates back to the early thirties of the last century.In their model supersaturated vapor containing monomers and sub-critical clusters of molecules collide to form stable clusters of molecules which exceed a critical radius.This critical radius depends on the surface energy, the volume of the formed structure, the temperature and the degree of supersaturation.The surface energy for clusters of such small size is not known and is usually estimated from bulk properties.These estimates can deviate from the true values by orders of magnitude.Unfortunately surface energies can only be measured for clusters which are significantly larger than the critical size and therefore cannot provide sufficient accuracy.However molecular dynamics simulation can be used to calculate the properties of a critical cluster, for instance surface tension.These simulations are computationally very expensive and contain significant approximations.
For the simulation of SiO 2 and TiO 2 the critical cluster size has been determined by homogeneous nucleation theory using bulk properties to be not more than one molecule in [60].However in a more recent paper Artelt et al. [6] studied the effect of varying surface energy on the critical particle diameter and proposed that the size of critical clusters can exceed one.
The nucleation of organic nanoparticles is believed to be different but also still unknown.The inception species for the first soot particles are polycyclic aromatic hydrocarbons (PAH).The model in [3], also called ABF model, which is widely used in the combustion community, assumes that the first particles form when two pyrene molecules collide.In this model every collision leads to successful formation of a soot particle.Despite the simplicity of this assumption the prediction for number density and soot volume fractions in premixed laminar flames are within an order of magnitude of the experimental observations.A more refined model is that of Frenklach and Wang [36] where a population balance of PAHs of different sizes is considered which form a soot particle on collision.
More recently, quantum mechanic calculations have been used to shed light on the problem of particle inception.Appel et al. [4] argue that dimers should be stable under flame conditions if the interaction energy in the dimer is larger than the internal energy of the two molecules.The result of their quantum mechanical calculations indicate that PAH dimers larger than the pyrene dimer are sufficiently stable to form a three dimensional structure at temperatures above 1500K.Schuetz and Frenklach [104] make use of molecular dynamics to study the collision of two pyrene molecules and conclude that aromatic dimers of species as small as pyrene can survive long enough to evolve into soot nuclei.Violi [118] combines kinetic Monte Carlo simulations on a mesoscopic level with molecular dynamics simulations to study the formation of the first soot particles from PAHs.She explains successfully the difference in H/C ratio, particle sphericity, and depolarization ratio as a function of the fuel properties.Despite first success of molecular dynamics modelling the time scales on which these calculations are carried out are in the order of nanoseconds which is too short to bridge the time gap necessary to establish the statistics of collision.Fortunately new experimental data on particle size distribution down to very small particle sizes are available for validation [130] and will help in gaining a better understanding of particle inception.

Coagulation
The model for particle coagulation depends on the choice of the particle model and on the physical conditions like temperature and pressure of the environment the particles are in.The coagulation kernel K determines the rate of collision of two particles.For nanoparticles even in turbulent flows Brownian coagulation is the mechanism with which particles coagulate.In [14] it is shown that for typical flame conditions the Brownian collision frequency is dominating over the Saffman Turner collision frequency which describes the collision frequency of particles due to turbulence.Brownian coagulation is modelled by the Fuchs kernel, e.g.[125,37,86] which covers different regimes classified by the Knudsen number, the ratio of twice the gas mean free path to the particle diameter.The free mean path of course depends on temperature and pressure of the particle surrounding gas.The different regimes are the continuum, slip-flow, transition and free molecular regime where at one end of the spectrum, in the continuum regime, particles perform a Brownian motion and on the other end of the spectrum, in the free molecular regime, particles perform a free flow followed by a collision type of movement.In addition to these kinetically derived collision frequencies an enhancement factor due to van der Waals forces between the particles needs to be accounted for [43].These forces lead to coagulation of particles even if they do not collide but get "close enough".
The collision diameter which appears in the coagulation kernel is unknown except for the very simple Coalescent Sphere Model where all particles are spheres.The collision diameter depends on the fractal dimension of the aggregate which is formed.The fractal dimension can be obtained from Monte Carlo simulations.In the continuum regime the diffusion limited cluster-cluster aggregation model (DLCA) gives a fractal dimension of 1.80.For the free molecular regime, the prevailing in atmospheric flames, the ballistic cluster-cluster aggregation model (BCA) leads to a slightly higher fractal dimension of 1.94.
However, sintering, surface growth, or condensation can change the fractal dimension significantly.Mitchell and Frenklach [75,76] developed a Monte Carlo method including coagulation and surface growth.Using the BCA model, i.e. assuming a particle sticks rigidly at the first point of contact to a target particle, they performed surface growth on each primary particle in the target particle by integration, losing the discrete nature of the surface events.The collision diameter and fractal dimension was then obtained from the target particle.
In [103] a Monte Carlo simulation of aggregation and sintering is performed.Similar to the approach above a target particle is allowed to coagulate with other particles.The authors then employ a 'shrinking back backbone' approach to change the surface to volume ratio and therefore the particles' fractal dimension.
Very recently Balthasar et al. [10] combine the approach developed in [76] with a Monte Carlo algorithm developed in [40] to simulation for the first time a full population of nanoparticles using the Primary Particle Model for a nonpremixed laminar flame.
The most important result of the above-mentioned simulations is that the fractal dimension of nanoparticles changes in time depending on surface growth, coagulation, and sintering.The detailed simulations can then be used to extract a functional relationship for the change of fractal dimension and therefore collision diameter with the aim to include this information in a simpler Surface and Volume Model.

Surface reaction
In order to model surface reactions accurately it is necessary to understand the gas-phase chemistry and the chemical composition and area of the surface of the nanoparticles.In many applications surface reactions play a very important role.Along with particle sintering which will be discussed later surface reaction can significantly change the fractal dimension of a particle.Unfortunately the importance on the one hand is reflected by the lack of knowledge about these processes on the other hand.
In the case of carbon nanoparticles the number of active sites has been observed to reduce as a function of the age of a soot particle in a premixed laminar flame.In [3] a measure for surface reactivity was introduced and fitted to a number of experiments.In [107] the surface reactivity has been associated to the age of the soot particles which accounts for the change of the surface chemistry of each particle.However this very simplistic model needs to be replaced by a more detailed one in which the chemical composition of the particle's surface is known.
In a series of papers [31,32,35] Frenklach and co-workers used a kinetic Monte Carlo technique to study this deactivation of surface reactivity on soot particles.The fluctuations of the spatial structure of PAHs on the surface of the particles are responsible for the deactivation of active surface sites.However these findings have not been implemented in a population balance model.
The number of papers which present models for surface growth of inorganic nanoparticles is quite small ( e.g.see [92] and references therein).Here, not only the chemical nature of the active sites needs to be known but also an accurate model for the change of the surface through sintering needs to be in place.

Condensation
Particle condensation is a process in which a chemical species from the gasphase sticks to the surface of a particle not through chemical reactions but through the weaker van der Walls force.A typical example for condensation is the attachment of PAHs on a soot particle's surface.In the case of soot or carbon black this can be a significant growth process.Typically species that play a role in the particle inception mechanism, like PAHs, are likely to condense on a particle's surface.For laminar flames where the particle inception zone is very thin condensation is less likely to play an important role in particle growth.

Sintering
The driving force for the sintering of nanoparticles is the reduction of potential energy due to the decrease in surface area.A lot of research has been carried out to understand sintering of metallic particles.However there much fewer papers on nanoparticles that simultaneously perform sintering and coagulation.
One of the most widely used sintering models in population balances has been introduced by Koch and Friedlander [49,128].The model is based on the Particle Volume Surface model where a particle with a surface to volume ratio equivalent to an agglomerate of primary particles changes its surface area at a characteristic sintering time to the surface area of a perfect sphere while keeping its density constant.For some nanoparticles, e.g.Si-particles the driving force for this physical mechanism is grain boundary diffusion.This simple model describes only the later stages of the sintering process sufficiently well where the particles are close to spherical.Structural changes are not taken into account.
To obtain deeper understanding into the structural changes, and therefore the changes in fractal dimension over time, Monte Carlo simulations have been performed by a number of research groups.Pratsinis and co-workers [1] carried out one of the first Monte Carlo simulation considering coagulation with simultaneous restructuring due to sintering.They examined two-dimensional clusters with finite binding energy using random particle walks on the surface of these clusters.This model is able to reproduce at least qualitatively the influence of residence time and temperature on nanoparticles formed in a reactor.Much more recently Schmid et al. [103] present a very detailed Monte Carlo simulation of agglomeration and sintering where the particles are shrunk along a common backbone and the overlapping volume is then distributed over the accessible surface.Depending on the sintering time the fractal dimensions vary from 1.86 to 2.99.They establish a functional relationship for the change of the fractal dimension over time and use this relationship in a simpler population balance model [5].

Numerical Methods
The numerical solution of equation ( 1) is a very challenging problem for several reasons.The two most important ones are the large size of the system of ordinary differential equations and the nonlinearity of the coagulation operator (2).For example in the case of the coalescent sphere model the difference between the smallest and largest particle sizes and therefore the size of (1) can easily be several orders of magnitude.Surface reactions and particle inception can also cause severe numerical challenges.In this section we shall discuss different strategies to obtain meaningful numerical solutions to (1) despite these difficulties.

Method of Moments
The method of moments (MOM) is computationally the most efficient approach to obtain a numerical approximation to the moments of a population balance.For this reason this method is often used when simulating problems where transport of particles in a flow with complex geometry is essential.In the area of nanoparticle modelling two techniques have been used so far.The first technique used by Frenklach and co-workers is the method of moments with interpolative closure (MOMIC) [34].The second approach is based on the quadrature method of moments which is a more recent technique based on the work of McGraw [73].
In the first instance MOMIC has been developed to describe the formation and oxidation of soot particles.In its early form the method is based on a univariate description of spherical soot particles in the free molecular regime, for instance in [36] MOMIC is used to simulate the formation of soot in a burner stabilized premixed ethylene flame where the soot moments are calculated in a post processing step.In this work two sets of moments of the size distribution of PAH and the size distribution of soot particles are solved simultaneously.Some theoretical remarks on the validity of the MOMIC approach can be found in [33].Numerical tests of the accuracy of the interpolative closure are contained in [12] for inorganic and in [41] for organic nanoparticles.For unimodal particle size distributions where no surface reactions that reduce particle size are present MOMIC produces excellent results being accurate and fast.The method has been extended to include coagulation in the transition and the continuum regime and taking into account the aggregate structure of soot particles by introducing a preset size for a primary particle in an agglomerate [47].In the latest development in MOMIC a shape descriptor, which is related to the surface fractal dimension, is used to model the aggregate structure of soot particles using particle volume and particle surface area [9].
An alternative approach for obtaining the moments of the PSD is the quadrature method of moments (QMOM) [73].There the moments are calculated assuming the PSD can be represented as weighted multi-dimensional Dirac delta function.The weights and the nodes are then chosen to satisfy the transport equations for the moments of the PSD.The advantage of this approach is that due to the choice of delta functions there exists no closure problem.Rosner and co-workers use this method to simulate uni-variate and bivariate population balances to model alumina nanoaggregate evolution in counterflow diffusion flames.The bivariate model includes coagulation and sintering of the particles.In [98] Rosner uses a Monte Carlo method to validate this approach and discusses the the role of mixed moments in QMOM.Despite the success of this method Marchisio and Fox [68] point out two main problems with QMOM.First, the solution of the transport equations for the moments of the PSD makes it difficult to treat systems where the dispersed-phase velocity strongly depends on the internal coordinates and more importantly for multivariate systems the QMOM becomes numerically very challenging.The authors present an improved approache called the 'direct quadrature method of moments' (DQMOM) which is computationally more attractive and allows the extension to more internal variables.Despite all the computational advantages there are some shortcomings associated with the method of moments.The most significant shortcoming is the non-uniqueness of the reconstruction of the particle size distribution function, i.e. the PSDF is not available.As a consequence, reactions that take place on the particles' surface which lead to decomposition of particles back to the gas phase can only incorporated with additional model assumptions.

Sectional Methods
Sectional methods (SM) are widely used to solve population balance equations.There the size spectrum is divided into a set of size classes.One distinguishes between zero-order and higher-order methods.Higher order methods use low order polynomials to represent the particles within each section and can be regarded as a simple form of finite element methods which will be discussed in the next section.They can suffer from stability problems and artificial dispersion whereas zero-order methods are more robust.There the computational domain is divided into rather small intervals in which the solution is approximated by step functions.For each interval one obtains an ordinary differential equation which is coupled to the neighbors depending on the discretization scheme used.Batterham et al. [16] divide the size domain in a geometric series of 2 and derive a numerical scheme which is mass conserving.Hounslow, et al. [44] extend this method using a correction to conserve particle number and particle mass.Litster et al. [62] modified this method by introducing an adjustable geometrical size discretization so that higher moments and self-preserving shape of the PSD is correctly predicted.However, there are large numerical errors in the presence of discontinuities and surface growth.To alleviate the problem with numerical diffusion in the presence of surface growth Kumar and Ramkrishna [54,55,56] introduced a pivot technique combined with a moving grid as well as the method of characteristics.An additional set of equations is solved to ensure that a chosen set of properties is conserved.All these methods are compared by Vanni [115] for a comprehensive set of test cases.Pope and Howard [89] couple a sectional model to a detailed gas phase mechanism to calculate the soot particle size distribution in a CSTR.Sectional methods have also been developed for bivariate population balance equations.Based on a model by Koch and Friedlander [49] Xiong and Pratsinis [128] present a two-dimensional particle size distribution including particle volume and particle surface as internal variables.The numerical method is based on a sectional approach which conserves volume.Nakaso et al. [79] also use a two-dimensional discrete-sectional representation of the size distribution solving the aerosol general dynamic equation for chemical reaction, agglomeration, and sintering.However, the use of these methods has been limited by their high computational cost.Hence in [78] a number of approximations to the detailed twodimensional model in [128] are made.The most important being the introduction of constant average coagulation coefficients and employing look up tables.Also a volume correction term was introduced that ensured the conservation of both number and mass.An alternative to bivariate population balances with volume and surface area as internal coordinates is to average the surface area information for each volume section turning a bivariate population balance into two univariate population balance.Different strategies are used to do this.Tsantilis and Pratsinis [114] proposed a sectional method in which they assumed that particle surface area in a given section decreases monotonically due to sintering and solved two sets of equations for particle volume and surface area using the sectional method developed in [44].This method is applied in [113] to simulate titania particles and compared with the bivariate model.Recently, Jeong and Choi [45] used a single surface fractal dimension to correlate particle volume and area in a given interval to simplify the bi-variate model of Xiong and Pratsinis [128] to two sets of one-dimensional equations for particle volume and surface area.A slightly different approach has been developed by Park and Rogak [87] where the agglomerate volume and the number of primary particles within a aggregate is tracked.This model considers the effect of surface growth due to ultrafine particle deposition on the primary particle size.Wen et al. [124] combined this with moving sections as well as the addition of new sections for incepted particles and applied it to model soot formation in a plug flow reactor.In each section the number of soot clusters and the number of primary particles is stored.

Finite Element Methods
An alternative to sectional methods are the more sophisticated finite element methods.In the finite element approach the solution of the population balance is expanded in a series of polynomials.For the coefficients of this expansion a set of equations has to be solved which is obtained by inserting the expansion into the population balance equation.Various methods can be derived by a different choice of basis functions, nodes, and time stepping schemes.The mathematical discipline of functional analysis provides the theoretical framework with which errors can be estimated.This is of course a very attractive feature of finite element methods.
In the early 1990s, Deufelhard, Bornemann, and Wulkow at the ZIB in Berlin [25,24,20] developed a discrete Galerkin h-p method using for the simulation of molecular weight distributions in polymerization reactions.Wulkow commercialized this idea and successfully implemented an algorithm in his software product PREDICI [126].
The structural similarities between the population balances in polymerization and nanoparticle dynamics then led to the development of a new software package PAR-SIVAL, a new dynamic flow sheet simulator especially suited for process units which are modelled by population balance equations.This new approach is based on a fully adaptive Galerkin h-p method with an automatic error control for time discretization as well as for the discretization of the property coordinate.The minimization of the required number of degrees of freedom is combined with an adaptivity in time and internal coordinate.Details of the numerical procedure are given in [127].Recently, PARCIVAL has been extended to solve bivariate population balances but no examples involving the simulation of nanoparticles are available in the open literature.
PARSIVAL and PREDICI are employed not only in industry but also by several research groups.Appel and Bockhorn used the PARSIVAL to simulate soot particle size distributions in laminar premixed flames [4].Artelt, Schmidt, and Peukert [5,6] have been using PARSEVAL for modelling Si production.They incorporated the mean value of fractal dimension into their model as well as the agglomerate collision diameter and the number of primary particles.Not using PARSEVAL or PREDICI but based on Wulkow's earlier work, Vlasov, Warnatz and co-workers employed a finite element code to simulate the formation of soot in a shock tube [119,80].
Independent of the work that originated in Berlin, also in the early 1990s, Sabelfeld and co-workers from Novosibirsk in Russia [101,52] developed a finite element approach using a local set of either B-spline or Hermitian spline basis functions.For the time stepping they employ an explicit scheme with enlarged stability region.This method has been applied to aerosol dynamics describing in SiH 4 decomposition in a quartz reactor [85].Nicmanis and Hounslow in Cambridge developed a finite element method to calculate the steady state of a general population balance equation [81,82].In [83] some error estimates are derived for this method.Although this method has not been applied to nanoparticle dynamics the authors make some comparison to sectional methods finding the finite element approach numerically more accurate.
More recently, Liu and Cameron [63] presented a wavelet based approach which has been validated against analytical solutions.Wavelets are a special choice of basis function and provide through their multiresolution properties a number of advantages.They seem particularly suited for modelling bimodal particle size distributions.However, up to now this technique has not been used for modelling nanoparticle population dynamics.

Monte Carlo Methods
An alternative to sectional methods for solving the population balance equation are Monte Carlo (MC) methods.They are easy to implement, can account for fluctuations, and can easily incorporate several internal coordinates.In the case of nanoparticle modelling the number of particles is so large that the fluctuations in particle numbers can be neglected.
Monte Carlo methods used to obtain approximations to the solution of equation ( 1) are all generalization of the classic Marcus-Lushnikov process [70,65].The theory of probabilistic methods related to coagulation is reviewed in [2].One of the first application of Monte Carlo methods for the solution of the Smoluchowski equation was presented [39].In this very influential work Gillespie used a Direct Simulation Monte Carlo (DSMC) method to simulate the dynamics of aerosol particles.Around the same time Ramkrishna and co-workers presented a Monte Carlo method for population balance equations [106,95].
The DSMC algorithm works as follows.First, a normalization volume is chosen then the solution to the population balance equation is approximated by N particles.The normalization volume plays the role of a numerical approximation parameter and determines the error between the approximation and the solution to (1).This particle system evolves in time mimicking the physical nanoparticle system.For a waiting time τ nothing happens.Then one of the following events, particle inception, coagulation/coalescence, surface growth, or sintering takes place depending on their probability, i.e. on the rate of the corresponding processes.The particle system is modified accordingly and the new probabilities are calculated.Then a new waiting time is determined.
Over the last decade the DSMC has been improved by several research groups in the engineering community.Most of these improvements can be classified as variance reduction techniques.Matsoukas and co-workers developed a constant number algorithm for coagulation [108], coagulation and fragmentation [58], and [59] contains an extension to other source terms and discusses a variation to the constant number approach namely a constant mass approach which is found to be superior.Maisels et al. [67] discuss a particle doubling strategy to reduce variance for simultaneous nucleation coagulation and surface growth.Independently of those efforts MC methods have been developed by mathematicians and physicists mainly to simulate problems arising in nuclear weapon development and space programmes.In Russia probability theory was applied to aerosol transport quite early.In the mid-eighties Sabelfeld and co-workers in Novosibirsk developed a weighted MC method that made use of an acceptance rejection method also called fictitious jumps to speed up the algorithm.The idea is based on the acceptance rejection technique introduced by John von Neumann where the coagulation kernel is replaced by a majorant.This leads to more events/jumps but some of them are rejected to correct the overall number.This method is computationally attractive if the majorant reduces the numerical effort when choosing the collision pairs.In [101] Sabelfeld and co-workers publish in English a variety of techniques that significantly improve the standard DSMC techniques.They not only use a constant majorant kernel but also discuss a more sophisticated version of the particle doubling technique, introduce a weighted particle method for variance reduction and discuss the Nanbu type stochastic algorithm.These methods are nothing else but a mathematically rigorous way of describing what is known in the engineering community as time driven, event driven, constant number and constant volume methods.
Another breakthrough has been achieved at WIAS in Berlin by Wagner and his student Andreas Eibeck.Motivated by Wagner's work on the Boltzmann equation they refined the method of fictitious jumps by introducing more effective majorants [28] for the coagulation Kernel and, more importantly, managed to prove rigorously the convergence [27,57] of their stochastic simulation procedure for (1).This algorithm will be called direct simulation algorithm (DSA) in this paper.Other interesting convergence results and algorithms are published [84,22,42,57].One choice of weighted MC methods deserves more attention.Based on work by Babovsky [7] Eibeck and Wagner developed a mass flow algorithm (MFA) [29] which is in many cases more efficient than DSA [17] and can be generalized to include particle inception, surface growths [23,77], and sintering [122].Kolodko and Sabelfeld, now also at WIAS, discuss MFA as a special case of a general class of stochastically weighted algorithms and perform some numerical experiments to compare the different methods [51].
Monte Carlo methods can be easily extended to multiple internal coordinates and for this reason they have been employed to simulate various systems of nanoparticles.Akhtar et al. [1] included two internal coordiantes to describe the restructuring of the particles and showed the time evolution of the mass fractal dimension.Tandon and Rosner studying a very similar system using the technique of fictitious jumps with a constant majorant for the coagulation sintering equation [112].Methods introduced by Wagner have been applied to inorganic nanoparicle dynamics in [41,40,123] and organic nanoparticles in [12,107].
Monte Carlo methods have also been used to model coating of nanoparticles.For example in [26] the enclosure of Fe 2 O 3 by SiO 2 is studied.There a population balance within a population balance is modelled.In [53] a DSMC algorithm is applied to model the coating of particles by aggregation.However, in the paper not much detail on the technique is given.
Monte Carlo methods have also been used to simulate the growth of single aggregates and study the effect of surface growth or sintering on the fractal dimension and other properties of the aggregate.For example, Mitchel and Frenklach [75] developed a Monte Carlo method to study soot agglomeration and surface growth of a single particle.A similar study but for inorganic aggregates which change through sintering has been performed in [103].So far a full simulation of aggregation of thousands of aggregates has to not been feasible due to excessive computational demands.However, for the first time, in [13] efficient majorants [40] and a ballistic MC algorithm including surface growth [76]have been combined to simulate a full population balance of soot particles in a premixed laminar flame.Figure 2 shows a two dimensional projection of a soot particle from a premixed laminar ethylene flame simulation.

Transport of Nanoparticles
So far we have only discussed models based on the assumption that the physical quantities do not vary in space.In practice this is of course never the case and therefore the models need to include the transport of particles due to convection and diffusion.Simulating reactive flows is a difficult task even if there are no particles present.We distinguish between laminar flow and turbulent flow.

Laminar flows
Transport of nanoparticles in laminar premixed flames has been modelled using the method of moments and sectional methods.Historically the effect of diffusion has been neglected.For stationary premixed laminar flames the space coordinate can be transformed into a time coordinate so that the zero-dimensional population balance equation ( 1) can be solved.If transport by diffusion and thermophoretic effects are included then additional equations that either account for the sections or the moments of the population balance have to be solved.The first coupling of a MOMIC approach to flow equations simulating a laminar counterflow diffusion flame was published by Mauss et al. [72] and was subsequently used also for premixed laminar flames in [71].Giesen et.al. [38] coupled the equations for the moments of the particle sizes and area to the commercial CFD package Fluent to simulate the formation of iron particles in a heated reactor.Marchisio et al. implemented QMOM [69] into Fluent taking simultaneous aggregation and breakage into account.Also the sectional approach has been coupled to flow equations by Smooke et al. [109] to simulate the formation of soot in a co-flow ethylene diffusion flame.The authors are using a sectional approach which is coupled to a detailed chemical model along with a velocity-vorticity formulation including buoyancy.Sun et al. [111] present a fully coupled two-dimensional sectional method where they model the formation of sodium coated carbon particles in a sodium halide diffusion flame.These simulations consume a lot of CPU time even on fast machines, in particular, if the gas phase chemical mechanism is very large, as it is the case for combustion of hydrocarbons leading to soot.A strategy to alleviate this problem is to decouple the fast processes, like the fast chemical reactions in the gas phase, from the slower ones.This was achieved for a laminar co-flow diffusion flame using a flamelet library for the chemical species and the source term of soot volume fraction.In the simulation of the flame only equations for soot volume fraction and mixture fraction had to be solved along with velocity and temperature [11].The authors report good agreement with measurements and low CPU times.

Turbulent flows
As most industrial reactors are turbulent it is important to incorporate a model for nanoparticle transport into a model which describes turbulent flows.Modelling of reactive turbulent flow is a very challenging problem due to its nonlinear, chaotic nature and the magnitude of length scales present in the flow.The books and review articles by Fox [30] and Pope [90,91] give an excellent treatment of the physics of turbulent flow and models which are in use.Although it is possible to formulate detailed models for the synthesis of nanoparticles in turbulent flames it is necessary to simplify these models to make them computationally tractable.The most detailed statistical approach is based on the probability density function (PDF) transport equation from which the moments of all physical quantities can be obtained.The PDF transport equation including an empirical model for soot formation and a global chemistry model has been applied to turbulent flames by Kollmann et al. [50].Borghi et al. uses a variant of the PDF approach [102] also with global chemistry and an empirical expression for soot volume fraction.More recently, a PDF approach with a reduced chemical mechanism with 15 species and 144 reaction, Curl mixing model and a radiation model has been used by Lindstedt and co-workers [61].The particulate phase is modelled using the method of moments including detailed surface reaction chemistry.
Almost all PDF based models use Monte Carlo techniques in conjunction with a CFD code and therefore the computational cost is very high.An alternative and much less computationally expensive are stochastic reactor models (SRMs) [17] which are based on the assumption that there are no spatial gradients in the PDF of the physical quantities in the computational domain.Because of this gross simplification, detailed reaction chemistry as well as the moments of a particle population can be included in SRMs making them a useful and computationally attractive tool to simulate engines [129] and industrial carbon black reactors [14].The alternative to SRMs, where the information on the flow field is sacrificed in favor of the chemistry, is to use a very simple chemistry model and solve Reynolds or Favre averaged transport equations using a closure model for the mean reaction rates.For example Lockwood [64] used a global chemistry step combined with an empirical model for soot mass fraction and soot number density employing a presumed beta PDF closure.
Soot formation in turbulent flames, gas turbines and diesel engines has also been modelled using a CFD flamelet approach combined with the method of moments.Two methods are used to decouple the chemistry from the flow field.Mauss and co-workers [8,46,15] use flamelet libraries which contain the species information as well as the source terms for surface growth, particle inception, and condensation.This approach is now implemented in the commercial code STAR-CD.A computationally more demanding but more accurate approach is to transport representative interactive flamelets (RIF) through the computational domain, average over their properties and couple these properties back to the CFD code.This method has been successfully used to simulate diesel engines by Pitsch et al. [88].
Nanoparticles in turbulent reactive flow have only been modelled using the method of moments incorporated into a CFD code.However for the transport of aerosols [100] and the transport of bubbles in an extraction column [116] a population balance has been coupled to a CFD code using Monte Carlo methods.These methods have good potential if it is necessary to couple population balances with several internal coordinates to turbulent reacting flows.

Applications
In this section we attempt to collect a list of references which report on numerical models for nanoparticle synthesis.Table 1 contains a list of applications where the numerical methods described above have been used to simulate nanoparticle population dynamics.As throughout this paper, both organic and inorganic nanoparticles have been treated together.In [48] a more detailed but older list on soot modelling can be found.The entries of the table are ordered with respect to the mathematical complexity of the model.For example, entries in the left column range from a shock tube (ST) to real industrial furnaces, gas turbines (GT), and internal combustion engines (ICE).It is important to note that this table is not complete and cannot be given the vast amount of literature.However, examples for the most important models and numerical methods are included.One important field of future research will be the chemistry of precursors to inorganic nanoparticles.Homogeneous and heterogeneous chemical reactions need to be understood to model particle formation and the change in particle size due to chemical reactions on the particles' surface.Currently, data on the thermodynamic properties of relevant species and reaction mechanisms are either missing or have high uncertainties associated with them.It will be necessary to conduct more high quality experiments in order to determine the relevant data.However, the role of computational chemistry will become more important as computers become more powerful so that thermodynamic properties and reaction kinetics can be obtained.
Another area which will certainly develop further is the coupling of population balances to CFD codes.Here the method of moments seems to be very promising but it is necessary to extend the method so that more internal variables can be tracked.The direct quadrature method of moments seems to be the method of choice at the moment.However, with the increase of computational power coupling Monte Carlo methods to CFD seems also viable.A more fundamental problem is the modelling of turbulent two phase flow which will remain a challenge for a long time to come.Apart from coupling of the population balance to the continuous phase, closure models for chemical reactions and for the turbulent transport terms need to be found.The use of LES in conjunction with flamelets and the PDF transport equations could be a possible but computationally expensive way forward.
An area which is very important is the development of mathematical techniques for estimating parameters in population balance models from experimental data.Recently, a weighted Monte Carlo method for the calculation of parametric derivatives has been developed [117].This method has a large potential to determine physical parameters form experimental data in the area of nanoparticles.
In the scientific community only the modelling of homogeneous gas phase nanoparticle synthesis has been studied.However in many applications the precursors of the particles are sprayed into a flame where they then form particles under the appropriate conditions [66].This means that the dynamics of the droplets including the precipitation within the droplets needs to be included in the population balance model.Models which have been developed in the area of spray drying are a suitable starting point to achieve this.
There exist a number of applications where nanoparticles are coated or are composites of different types of particles.Models for these processes have not been treated in this paper but will be an important field in the future.First steps have been taken in [111] where the primary size of a particle has been controlled by condensing a sodium fluoride layer on the particles and thus avoiding the formation of hard agglomerates.The authors of that paper used a sectional method with two internal coordinates to model this process.Monte Carlo methods have been used in similar applications [26].

Conclusion
In this paper we reviewed the literature of nanoparticles from flame synthesis and discussed the models and the numerical methods used in this field.We started by describing a general population balance model for nanoparticles in the absence of spatial gradients of the physical quantities.The model is comprised of several submodels which are: a detailed homogeneous gas phase reaction mechanism, a particle inception model, a model for particle agglomeration, surface reactions, condensation, and sintering.References for these submodels were given and their accuracy were discussed.Then the numerical treatment of population balances was reviewed.The literature on the method of moments, sectional, finite element, and Monte Carlo methods were discussed.Special attention was given to the latest developments in the area of Monte Carlo methods.They are very promising for solving population balance models with a very large number of particle properties.Then the coupling of a population balance to laminar and turbulent flow was examined.The method of moments has been identified as a suitable method to efficiently simulate laminar flames.The method of moments is also a good choice to couple a population balance to turbulent flow where several closure problems have to be addressed.The separation of time scales is an important technique when modelling particle synthesis in industrial type combustion chambers.A number of examples were summarized in a table to make the vast amount of literature more accessible.Finally, future areas of research were discussed.
New materials and processes as well as the ever growing computational power will turn modelling of particular processes from specialized software used in research laboratories to a development tool for engineers in industry.This development needs to be accounted for in the education of engineers.

Figure 1 :
Figure 1: Schematic of processes in nanoparticle flame synthesis

Figure 2 :
Figure 2: Two-dimensional TEM-style projections of a sample particle taken from premixed laminar ethylene flame simulation (Provided by Neal Morgan).

Table 1 :
List of applications.