Stochastic modelling of silicon
nanoparticle synthesis
William Jefferson Menz
Sidney Sussex College
October 15, 2013
A dissertation submitted for the degree of Doctor of Philosophy
Stochastic modelling of silicon nanoparticle synthesis
William J. Menz
This thesis presents new methods to study the aerosol synthesis of nano
particles and a new model to simulate the formation of silicon nanoparticles.
Population balance modelling is used to model nanoparticle synthesis and a
stochastic numerical method is used to solve the governing equations. The
population balance models are coupled to chemical kinetic models and offer
insight into the fundamental physiochemical processes leading to particle
formation.
The first method developed in this work is a new mathematical ex
pression for calculating the rate of Brownian coagulation with stochastic
weighted algorithms (SWAs). The new expression permits the solution
of the population balance equations with SWAs using a computationally
efficient technique of majorant rates and fictitious jumps. Convergence
properties and efficiency of the expression are evaluated using a detailed
silica particle model.
A sequentialmodular algorithm is subsequently presented which solves
networks of perfectly stirred reactors with a population balance model using
the stochastic method. The algorithm is tested in some simple network
configurations, which are used to identify methods through which error in
the stochastic solution may be reduced. It is observed that SWAs are useful
in preventing accumulation of error in reactor networks.
A new model for silicon nanoparticle synthesis is developed. The model
includes gasphase reactions describing silane decomposition, and a detailed
multivariate particle model which tracks particle structure and composition.
Systematic parameter estimation is used to fit the model to experimental
cases. Results indicated that the key challenge in modelling silicon systems
is obtaining a correct description of the particle nucleation process.
Finally, the silicon model is used in conjunction with the reactor net
work algorithm to simulate the startup of a plugflow reactor. The power
of stochastic methods in resolving characteristics of a particle ensemble is
highlighted by investigating the number, size, degree of sintering and poly
dispersity along the length of the reactor.
Preface
This dissertation is the result of my own work and includes nothing which
is the outcome of work done in collaboration, except where specifically in
dicated in the Acknowledgements. The work was performed at the Depart
ment of Chemical Engineering and Biotechnology at the University of Cam
bridge between October 2010 and June 2013. No other part of this thesis
has been submitted for a degree.
The thesis contains 42 figures and 15 tables and approximately 30,000
words. Parts of this thesis have been published in the following articles:
• W.J. Menz, S. Shekar, G.P.E. Brownbridge, S. Mosbach, R. Körmer,
W. Peukert and M. Kraft ‘Synthesis of silicon nanoparticles with a nar
row size distribution: A theoretical study’ Journal of Aerosol Science,
44:46–61, 2012.
• W.J. Menz and M. Kraft ‘A new model for silicon nanoparticle synthe
sis’ Combustion and Flame, 160:947–958, 2013.
• W.J. Menz and M. Kraft ‘The suitability of particle models in capturing
aggregate structure and polydispersity’ Aerosol Science and Technology,
47:734–745, 2013.
• W.J. Menz, R.I.A. Patterson, W. Wagner and M. Kraft ‘Application of
stochastic weighted algorithms to a multidimensional silica particle
model’ Journal of Computational Physics, 248:221–234, 2013.
• W.J. Menz, J. Akroyd and M. Kraft ‘Stochastic solution of popula
tion balance equations for reactor networks’ Journal of Computational
Physics, 256:615–629 2014.
William J. Menz
October 15, 2013
iv
Acknowledgements
I would like to acknowledge the members of the Computational Modelling
Group for their assistance: in particular, Prof. Markus Kraft (supervisor)
for his guidance and Alastair Smith for tirelessly maintaining the group’s
computer resources.
The initial work towards this thesis was supported by collaboration with
Dr. Richard Körmer and Prof. Wolfgang Peukert at the FredrichAlexander
Universität ErlangenNürnberg. Their input on early silicon modelling at
tempts was invaluable.
The numerical aspects of the population balance modelling was devel
oped with insight from Dr. Robert Patterson and Prof. Wolfgang Wagner at
the Weierstrass Institute of Applied Analysis and Stochastics, Berlin. I would
like to acknowledge them both for many useful discussions and some coding
help.
Most importantly, I thank the Cambridge Australia Trust for providing
the financial support, without whom my coming to Cambridge would not
be possible.
v
Contents
Summary iii
Preface iv
Acknowledgements v
List of Figures ix
List of Tables xi
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Novel aspects of the thesis . . . . . . . . . . . . . . . . . . . . 3
1.3 Structure of the thesis . . . . . . . . . . . . . . . . . . . . . . 4
2 Background 5
2.1 Gasphase chemical kinetics . . . . . . . . . . . . . . . . . . . 5
2.1.1 Rate expressions . . . . . . . . . . . . . . . . . . . . . 6
2.1.2 Threebody reactions . . . . . . . . . . . . . . . . . . . 7
2.1.3 Falloff reactions . . . . . . . . . . . . . . . . . . . . . . 7
2.1.4 Thermochemistry . . . . . . . . . . . . . . . . . . . . . 8
2.1.5 Gasphase library . . . . . . . . . . . . . . . . . . . . . 8
2.2 Population balance modelling . . . . . . . . . . . . . . . . . . 8
2.2.1 The collision kernel . . . . . . . . . . . . . . . . . . . . 11
2.2.2 Stochastic solution of the population balance equations 12
2.2.3 Other solution methodologies . . . . . . . . . . . . . . 15
2.3 Reactor models . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.1 Gasparticle coupling . . . . . . . . . . . . . . . . . . . 17
2.4 Error measurements . . . . . . . . . . . . . . . . . . . . . . . 17
vi
2.5 Silicon nanoparticles . . . . . . . . . . . . . . . . . . . . . . . 19
2.5.1 Manufacture . . . . . . . . . . . . . . . . . . . . . . . 19
2.5.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . 21
3 Choosing a particle model 23
3.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.1 Typespace . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.2 Derived properties . . . . . . . . . . . . . . . . . . . . 26
3.1.3 Particle processes . . . . . . . . . . . . . . . . . . . . . 28
3.1.4 Average properties . . . . . . . . . . . . . . . . . . . . 30
3.2 Case studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.1 Coagulation and sintering . . . . . . . . . . . . . . . . 31
3.2.2 Effect of polydispersity . . . . . . . . . . . . . . . . . . 35
3.2.3 Computational resources . . . . . . . . . . . . . . . . . 37
3.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4 A new coagulation kernel for stochastic weighted algorithms 40
4.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.1.1 Typespace . . . . . . . . . . . . . . . . . . . . . . . . 42
4.1.2 Particle processes . . . . . . . . . . . . . . . . . . . . . 42
4.2 Development of stochastic weighted algorithms . . . . . . . . 45
4.2.1 Linear particle processes . . . . . . . . . . . . . . . . . 46
4.2.2 Coagulation process . . . . . . . . . . . . . . . . . . . 46
4.3 Numerical studies . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.4.1 Highprecision solution . . . . . . . . . . . . . . . . . 56
4.4.2 Convergence of algorithms . . . . . . . . . . . . . . . 58
4.4.3 Computational efficiency . . . . . . . . . . . . . . . . . 61
4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5 Solution of reactor networks 66
5.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.1.1 Stochastic method . . . . . . . . . . . . . . . . . . . . 69
5.1.2 Discrete Model . . . . . . . . . . . . . . . . . . . . . . 72
5.1.3 Solution of reactor networks . . . . . . . . . . . . . . . 73
5.2 Numerical studies . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.2.1 Comparison to Discrete Model . . . . . . . . . . . . . . 74
vii
5.2.2 Linear networks . . . . . . . . . . . . . . . . . . . . . 76
5.2.3 Networks with recycle loops . . . . . . . . . . . . . . . 80
5.2.4 Reactors with multiple inlets . . . . . . . . . . . . . . 82
5.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.3.1 Approximation of a PFR . . . . . . . . . . . . . . . . . 84
5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
6 A detailed model for silicon nanoparticle synthesis 88
6.1 Model development . . . . . . . . . . . . . . . . . . . . . . . 91
6.1.1 Gasphase kinetic model . . . . . . . . . . . . . . . . . 91
6.1.2 Particle model . . . . . . . . . . . . . . . . . . . . . . . 92
6.1.3 Particle processes . . . . . . . . . . . . . . . . . . . . . 94
6.1.4 Coupling algorithm . . . . . . . . . . . . . . . . . . . . 98
6.2 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . 99
6.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 100
6.3.1 Numerical results . . . . . . . . . . . . . . . . . . . . . 100
6.3.2 Case studies . . . . . . . . . . . . . . . . . . . . . . . . 102
6.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
7 Putting it all together 114
7.1 Particle model analysis . . . . . . . . . . . . . . . . . . . . . . 114
7.2 Unsteadystate response . . . . . . . . . . . . . . . . . . . . . 118
7.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
8 Conclusions 122
8.1 Findings of this work . . . . . . . . . . . . . . . . . . . . . . . 122
8.2 Suggestions for future work . . . . . . . . . . . . . . . . . . . 124
A Appendix to Chapter 3 126
A.1 Collision of particles in the nearspherical zone . . . . . . . . 126
A.2 Comparison of maximum sinterable surface area . . . . . . . 127
B Appendix to Chapter 6 129
B.1 Conversion of SURFACE CHEMKIN parameters . . . . . . . . . 129
Nomenclature 131
viii
List of Figures
2.1 Illustration of population balance solution methodologies . . 12
2.2 Flowchart of the stochastic particle algorithm . . . . . . . . . 13
2.3 Schematic of the architecture of MOPS . . . . . . . . . . . . . 17
2.4 Temporal evolution of a macroscopic quantity in a particle
system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.5 SEM image of highlyspherical silicon nanoparticles . . . . . . 20
3.1 An illustration of different particle models’ description of co
agulation events . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.2 Comparison of normalised mean collision diameter as a func
tion of coagulation and sintering times . . . . . . . . . . . . . 33
3.3 Contours of difference between the binarytree and surface
volume models . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.4 Influence of polydispersity on particle morphology . . . . . . 36
3.5 CPU time required for different particle models . . . . . . . . 38
4.1 Comparison of the coagulation process implemented in direct
simulation algorithm and weighted particle methods. . . . . . 53
4.2 PSDs for the silica test case . . . . . . . . . . . . . . . . . . . 56
4.3 Temporal evolution of functionals describing the silica test case 57
4.4 Final values and confidence intervals of the key metrics for
the silica model . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.5 Relative statistical error in the final values of the key metrics
for the silica model . . . . . . . . . . . . . . . . . . . . . . . . 60
4.6 Computational time required for the coagulation algorithms . 62
4.7 Temporal evolution of number of coagulation events . . . . . 63
5.1 Algorithm used to solve the reactor network with Strang split
ting. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
ix
5.2 Schematic of the test systems used in the present study. . . . . 75
5.3 Transient evolution of moments for the threereactor system . 77
5.4 Convergence of metrics for a linear network . . . . . . . . . . 78
5.5 Comparison of statistical error in process metrics . . . . . . . 79
5.6 Convergence of metrics for fR = 0.5 . . . . . . . . . . . . . . 81
5.7 Statistical error as a function of recycle fraction . . . . . . . . 82
5.8 Statistical error accumulation through a reactor network . . . 83
5.9 Convergence of network solutions to PFR solutions . . . . . . 85
5.10 Comparison of network and PFR solutions . . . . . . . . . . . 86
6.1 Unimolecular rate constants for silicon decomposition . . . . 103
6.2 Experimental temperature profiles used in this work . . . . . 104
6.3 Comparison of experimental and theoretical PSDs for the case
of Körmer et al. [75] . . . . . . . . . . . . . . . . . . . . . . . 106
6.4 Comparison of experimental and theoretical PSDs for the case
of Frenklach et al. [49] . . . . . . . . . . . . . . . . . . . . . . 107
6.5 Comparison of experimental and theoretical PSDs for the case
of Wu et al. [163] . . . . . . . . . . . . . . . . . . . . . . . . . 108
6.6 Particle diameter dependence on silane for the system of Nguyen
and Flagan [103] . . . . . . . . . . . . . . . . . . . . . . . . . 109
6.7 Comparison of experimental and theoretical PSDs for the case
of Nguyen and Flagan [103] . . . . . . . . . . . . . . . . . . . 110
6.8 Temporal evolution of diameter for the case of Onischuk et al.
[113] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.9 Comparison of experimental and computergenerated TEM
images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.10 Temporal evolution of particle diameter in a laserdriven flame 112
7.1 Evolution of process metrics for the silicon model . . . . . . . 116
7.2 Trajectory of the case of Wu et al. [163] through the zones of
Figure 3.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
7.3 Transient evolution of polydispersity and sintering level . . . 119
7.4 Transient evolution of collision diameter and number concen
tration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
A.1.1Ratio of collision diameters, as predicted by different particle
models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
x
List of Tables
3.1 Comparison of derived properties of the particle models . . . 26
4.1 Particle selection properties for the weighted transition kernel 51
4.2 Definition of the coagulation weight transfer functions. . . . . 51
4.3 Numerical and model parameters used in the silica test case . 54
4.4 Description of functionals used to investigate the silica model 55
5.1 Summary of process metrics studied in the present work. . . . 75
6.1 Gasphase mechanism used for silane decomposition . . . . . 93
6.2 Heterogeneous reaction processes for silicon nanoparticles
synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.3 Numerical and model parameters used in the silicon model . . 98
6.4 Experimental data used to fit the silicon model . . . . . . . . 100
6.5 Results of the parameterfitting procedure . . . . . . . . . . . 101
6.6 Comparison of rate expressions for silane decomposition . . . 102
6.7 Process conditions of experimental test cases . . . . . . . . . . 104
B.1 Surface reaction mechanism of Ho et al. [63] . . . . . . . . . 130
B.2 Estimation of surfacereaction preexponential for silicon . . . 131
xi
Chapter 1
Introduction
It was not until the late 80’s that the term nanotechnology entered the scien
tific vernacular [35]. Since then, the field has rapidly subdivided and mul
tiplied in size, like doomsdaybringing ‘grey goo’. One of the biggest, and
possibly most wellknown areas of nanotechnology research is that which
focuses on nanoparticles, defined as particles of diameter 1–100 nm.
One of the earliest identified uses of nanoparticles was as part of a dec
orative glaze applied to ceramic earthenware [39]. Gold, silver and copper
‘nanoclusters’–formed during firing of the object–gave the ceramic its char
acteristic lustre. The nanoparticle glazes were made as early as the 10th
century AD. Today, nanoparticles are both an engineered product for indus
trial or research purposes, and an undesirable byproduct of human activity.
While it is common to think nanoparticles are an invention of humans,
there are many which occur naturally. Fires and volcanic eruptions produce
soot, which is composed of many different nanosized moieties [105]. Phys
ical and chemical processes in the soil may form inorganic nanoparticles [9],
and organic structures such as viruses are often considered to be a type of
nanoparticle [43].
Regardless of their origin, nanoparticles are of intense interest due to
their unique properties. The small size and porous or aggregate structure of
particles yield a large surface area to volume ratio [141], one of the factors
to which these interesting properties are attributed. Quantum effects can
come into effect at the ‘nano’ lengthscale, in which strange but potentially
useful characteristics manifest [80, 97].
The work in this thesis is directed towards understanding how nanopar
ticles can be made. In particular, population balance modelling will be used
1
to simulate the physiochemical processes leading to formation of nanoparti
cles in the gasphase. New models and methods will be developed to eluci
date the pathways to synthesising silicon nanoparticles, a material of consid
erable interest to industrial and research communities worldwide [88]. The
motivation behind this will be discussed in Section 1.1, and the structure of
the thesis will be presented in Section 1.3.
1.1 Motivation
The synthesis of nanoparticles in the gasphase occurs over very short time
scales, on the order of microseconds. Therefore, it is difficult to characterise
the processes responsible for particle formation with standard experimental
techniques [128]. Modelling particle synthesis may be used as an alter
native to experimentally characterising the formation processes. However,
the challenge in obtaining experimental data can leave uncertainties in the
kinetic expressions of the model. Such models typically utilise parameter
estimation to arrive at approximate rates of particle processes [144].
As the properties of nanoparticles are strongly dependent on their size
and morphology [2], it has become necessary to develop more detailed mod
els which betteraccount for the physiochemical processes which yield the
nanoparticle product. These models can potentially enhance process control
[56], reduce waste [14] and improve product properties.
Models for nanoparticle synthesis could also be used to investigate where
particles are not formed. Such knowledge could be used to suppress particle
formation in conditions where inhalation of particles is a hazard [28, 106],
or in conditions where particles are regarded as a contaminant in the man
ufacturing process [104].
The principal motivation of this work is to develop new methods and
models for simulating silicon nanoparticle synthesis. Traditionally consid
ered as an unwanted contaminant in semiconductor manufacture [41], they
are now being evaluated as printable inks [7], medical imaging tracers [86]
and in many other applications (see review by Mangolini [88]). The im
provements to methods for modelling nanoparticle systems presented in this
thesis will expand the flexibility of existing tools for solving the population
balance equations. These methods are not restricted to the particular case
of silicon.
2
1.2 Novel aspects of the thesis
In the first part of this thesis, new tools and methods for understanding
nanoparticle synthesis with population balance models will be developed.
Following that, a new model for silicon nanoparticle synthesis will be pre
sented, and the new modelling tools will be used in conjunction with the
silicon model to explore properties of an experimental system. Novel as
pects this thesis include:
• A comparison of particle models commonly used in population bal
ance modelling of nanoparticle formation. It is demonstrated that sub
stantial errors can be incurred when using a particle model inappro
priate for a particular modelling application.
• New insight into the structure and morphology of aggregate particles
formed due to coagulation and sintering.
• A new coagulation kernel for stochastic weighted algorithms (SWAs).
A detailed study of the numerical behaviour of the new kernel with a
multidimensional silica particle model is conducted.
• A sequential modular approach to solve a generic network of reactors
with a population balance model using a stochastic numerical method.
The convergence of the stochastic particle algorithm in test networks
is evaluated as a function of network size, feedback and numerical
parameters.
• New methods for accounting for particle transport to and from
reactors with stochastic particle methods. The optimal transport and
coagulation algorithms for some reactor systems are identified, yield
ing more efficient simulation of reactor networks.
• A detailed multivariate silicon model for the simulation of silicon
nanoparticle synthesis from silane decomposition. The model is tested
against a six different experimental configurations, with excellent fit
observed for the majority of cases.
3
1.3 Structure of the thesis
Chapter 2 gives a review of the theory behind modelling gasphase chem
ical kinetics, population balance modelling (particularly with reference to
stochastic methods), and coupling of the two systems. The chapter is con
cluded with an overview of the methods for manufacture and potential fu
ture applications of silicon nanoparticles.
A mathematical description of particle models used in population bal
ance modelling is provided in Chapter 3. Three models commonly used in
the literature are compared under different conditions of coagulation and
sintering. The results of this chapter highlight the importance of using mul
tivariate particle models to capture the structure and morphology of aggre
gate particles.
In Chapter 4, a new mathematical expression is presented for calculating
the rate of coagulation when using SWAs to solve the population balance
equation. It is tested using a multivariate silica particle model, and the
convergence and efficiency of the new expression are evaluated.
An algorithm for solving reactor networks with a population balance
model using stochastic methods is presented in Chapter 5. The algorithm
is tested for some simple reactor network configurations. These are used
to identify methods through which systematic and statistical error in the
stochastic solution may be reduced.
The formation of silicon nanoparticles is investigated in Chapter 6, where
a new model for silicon nanoparticle synthesis is presented. Parameter es
timation is used to fit the model to experimental cases. The new model
is evaluated against a variety of experimental configurations at different
process conditions in order to examine the goodness of fit. The physical
significance of the model is also assessed.
The methods and models developed in this thesis will be tiedtogether
in Chapter 7. First, the new silicon mechanism is analysed using the particle
models discussed in Chapter 3. Secondly, the new coagulation kernel of
Chapter 4 and reactor network algorithm of Chapter 5 are used with the
new silicon model to simulate the startup of a plugflow reactor.
The findings of the thesis are concluded and suggestions for future work
are presented in Chapter 8.
4
Chapter 2
Background
This chapter provides a theoretical background for the modelling work con
ducted in this thesis. In Section 2.1, the principles of kinetic modelling are
outlined, which form the basis of the silicon chemical mechanism devel
oped in Chapter 6. The population balance equations and the stochastic
methodology to solve them are introduced in Section 2.2, which are used
throughout this work for modelling nanoparticle formation.
The equations which compose the reactor models, and the theory of cou
pling the gasphase and particle solvers are given in Section 2.3. Numerical
metrics to measure the accuracy and precision of the stochastic population
balance solver, particularly relevant to Chapters 4 and 5, are discussed in
Section 2.4. Finally, a literature review of the manufacture and applications
for silicon nanoparticles is given in Section 2.5, which will be built upon in
Chapters 6 and 7.
2.1 Gasphase chemical kinetics
Kinetic modelling is used in order to understand the rate at which a chemi
cal species converts into other chemical species. The collection of gasphase
chemical species and reactions will be referred to throughout this text as
the ‘gasphase mechanism’ or the gasphase statespace, G. Combining a
mechanism with a reactor model (discussed in Section 2.3) and solving the
resultant system of ordinary differential equations (ODEs) yields the tempo
ral evolution of chemical species concentrations.
5
2.1.1 Rate expressions
Chemical reaction i in Nrxn reactions may be represented in the general form
[71]:
NGP∑
k=1
ν′ki χk
NGP∑
k=1
ν′′ki χk (2.1)
where ν′k and ν
′′
k are the forward and reverse stoichiometric coefficients of
species k for NGP chemical species with chemical symbols χk. k will be
used here to index the gasphase species. The production rate of species
k is given by the sum of the rate of all forward and backwards reactions
involving species k:
ω˙k =
Nrxn∑
i=1
νkiqi (2.2)
where νki is the net production rate:
νki = ν
′′
ki – ν
′
ki (2.3)
and qi is the rate of progress variable for the ith reaction. The form of the
rate of progress variable depends on the type of reaction under considera
tion. Generally, it is given by:
qi = kfi
NGP∏
k=1
c(χk)
ν′ki – kri
NGP∏
k=1
c(χk)
ν′′ki (2.4)
where c(χk) is the molar concentration (mol/m
3) of chemical species k,
kfi is the forward rate constant and kri is the reverse rate constant. The
bestknown form for the forward rate constant is the modifiedArrhenius
equation:
kf = A T
β exp
(
–
EA
R T
)
(2.5)
where A, β and EA are constants specific to the reaction. If reaction param
eters are not specified for the reverse reaction, the reverse rate coefficient is
given by:
kr = kf/Kc (2.6)
where Kc is the equilibrium constant for the reaction, which can be obtained
from the thermodynamic properties of the reacting species [71].
6
2.1.2 Threebody reactions
A third body, or collision partner, is often necessary for some reactions to
proceed [63, 71]. These bodies are typically denoted by a generic collision
species with chemical symbol M. An example of how such reactions can be
written is:
χ1 + χ1 + M → χ2 + M (2.7)
where χ1 and χ2 are reactant and product chemical species, respectively.
In this work, it is assumed that all chemical species in a mixture con
tribute equally as third bodies. This results in the following formulation of
the rate of progress variable:
qi =
NGP∑
k=1
c(χk)
kfi NGP∏
k=1
c(χk)
ν′ki – kri
NGP∏
k=1
c(χk)
ν′′ki
(2.8)
The summation term is often represented as the concentration of species M:
c(M) =
NGP∑
k=1
c(χk)
(2.9)
The calculation of the forward and reverse rate coefficients remains un
changed.
2.1.3 Falloff reactions
Equations (2.4) and (2.5) for standard Arrhenius reactions are typically used
to describe reaction rates at the highpressure limit, where the rate is depen
dent on temperature. However, threebody reactions imply a dependence on
mixture pressure due to the c(M) term in the rate of progress variable. This
condition is met when the system is near the lowpressure limit [71].
Expressions for the rate constant which span the low and highpressure
limits are termed falloff reactions, as they include a blending function to
account for the overlap between regions [71]. A falloff reaction is denoted
by a (+M) term in a chemical reaction:
χ1 + χ1 (+M) → χ2 (+M) (2.10)
There are a variety of expressions to account for the reaction rate in
the falloff region [87, 148, 151]. The equation of Lindemann et al. [87]
7
is the simplest and is used in this thesis. Two sets of modifiedArrhenius
parameters are first defined: those for the highpressure limit (k∞) and the
lowpressure limit (k0) according to Equation (2.5). The rate constant at
any pressure is then given by:
k = k∞
(
Pr
1 + Pr
)
(2.11)
where Pr is the reduced pressure, given by:
Pr =
k0c(M)
k∞
(2.12)
2.1.4 Thermochemistry
In order to calculate equilibrium coefficients, heat capacities, and other ther
modynamic properties of chemical species, reliable thermochemical data
is required. Such data is readily available online [150] or in other peer
reviewed literature (e.g. [149]). Alternatively, it is possible to derive this
data through use of quantummechanical calculations and statistical me
chanics (e.g. [126, 145]). In this work, thermochemical data in the form of
NASA polynomials is used [20].
2.1.5 Gasphase library
It is necessary to solve a system of many ODEs when investigating the chem
ical behaviour system with many reactions. Before this can be done, the
mechanism must be interpreted by a gasphase chemical kinetics library,
which supplies the rate of progress and mixture properties to an ODE solver.
This work uses an inhouse built library named Sprog to determine chem
ical reaction rates and gasphase mixture properties. Sprog is written in
C++ and has been validated against other popular libraries such as CHEMKIN
[71, 132] and Cantera [24].
2.2 Population balance modelling
The nucleation and growth of particles may be simulated with a population
balance model. In order to formulate the population balance equations,
which govern the behaviour of such models, it is necessary to mathemat
ically describe a particle. This representation is called the typespace of a
8
particle, denoted as P. The collection of all particles and the information
used to describe them is called the particle statespace, Z.
The evolution of the number density of particles n(t, Pq) at time t and of
type Pq ∈ Z is described by a population balance equation:
d
dt
n(Pq) = I(c, Pq) +K(n, Pq) +R(c, n, Pq) + S(n, Pq) (2.13)
where I(c, Pq), K(n, Pq), R(c, n, Pq) and S(n, Pq) are the inception, coag
ulation, surface growth and sintering operators of particles of type P, re
spectively. The index q will be used to index particles throughout this work.
The variable c is the vector of gasphase species concentrations and n is
the vector of particle concentrations. Each of these processes represents a
reaction which alters typespace of a particle, however the actual transfor
mation of the particle depends on the particle model used. These reactions
and operators are explained in a general form below.
Inception
Particles are created in the population balance by an inception process.
The terms inception and nucleation will be used interchangeably in
this thesis. This is represented by the following reaction:
χi + χj → Pq (2.14)
where χi and χj are nucleating gasphase species. The rate of nucle
ation is described by the inception operator I(c, Pq):
I(c, Pq) = 12
∑
χi,χj∈G;
χi+χj=Pq
K(χi,χj)c(χi, t)c(χj, t) (2.15)
where K is the collision kernel, a function which defines the rate of
collision, for the process. Equation (2.14) defines the + operator in
the summation term of Equation (2.15). This representation will be
used for the subsequent particle processes.
Coagulation
Particles coagulate and stick together according to the following reac
tion:
Pi + Pj → Pq (2.16)
9
The rate of coagulation depends on the choice of coagulation kernel
K, which defines the probability of collision. Particles of type Pi co
agulate with particles of type Pj at a rate K(Pi, Pj) according to the
Smoluchowski equation, such that:
K(n, Pq) = 12
∑
Pi,Pj∈Z;
Pi+Pj=Pq
K(Pi, Pj)n(Pi, t)n(Pj, t)
– n(Pq, t)
∑
Pi∈Z
K(Pq, Pi)n(Pi, t)
(2.17)
The first term in Equation (2.17) represents the formation of particles
due to the collision of smaller particles. The second term represents
the depletion of particles due to coagulation with larger particles.
Surface growth
Surface growth refers to the addition of mass to particles from a dif
ferent phase. This is typically described as a collisional process of a
particle with a gasphase species (condensation) or a heterogeneous
reaction of a gasphase species with the particle surface (surface reac
tion). Such events occur in the following manner:
Pi + χj → Pq + χk (2.18)
The surface growth operator is given by:
R(c, n, Pq) =
∑
Pi∈Z;
χj∈G;
Pj+χi=Pq
RSG(χj, Pi)n(Pi, t)c(χj, t)
–
∑
χj∈G
RSG(χj, Pq)n(Pq, t)c(χj, t)
(2.19)
where RSG is the surface growth rate:
RSG(χj, Pq) =
{
K(χj, Pq) for condensation
RSR(χj, Pq) for surface reaction
RSR is the surface reaction rate, given by a Arrhenius expression, e.g. Equa
tion (2.5).
10
Sintering
Sintering refers to the reduction in particle surface area due to internal
rearrangement of the component atoms or molecules. The form of a
sintering reaction is strongly dependent on the particle model chosen.
In general, they take the form:
S(Pi)
ς(Pq)−−−→ S(Pq) (2.20)
where S is the surface area of the particle and S(Pq) ≤ S(Pi) through
transformation ς. The sintering operator is given by:
S(n, Pq) =
∑
Pi∈Z;
ς(Pi)=Pq
Rsint(Pi)n(Pi, t)
– Rsint(Pq)n(Pq, t) (2.21)
where Rsint is the rate of sintering for the particle.
2.2.1 The collision kernel
Equations (2.15), (2.17) and (2.19) utilise a ‘coagulation kernel’ to esti
mate the collision probability of two entities. The transition kernel Ktr, a
computationallyefficient approximation of Brownian coagulation [70, 121],
will be used in this thesis. It is given by:
Ktr(Xi, Xj) =
(
1
Kfm(Xi, Xj)
+
1
Ksf(Xi, Xj)
)–1
(2.22)
where X ∈ G, Z; either a gasphase or particle species. Kfm is the free
molecular kernel and Ksf is the slipflow kernel. The dominant kernel is
dependent on the Knudsen number (Kn). The freemolecular kernel is dom
inant for Kn 1 [70] and takes the form:
Kfm(Xi, Xj) = kfm
√
pikB T
2
(
1
mi
+
1
mj
) 1
2 (
di + dj
)2 (2.23)
where kB is Boltzmann’s constant, T is the temperature and mi and di are the
mass and collision diameter of species i, respectively. kfm is the free molec
ular enhancement factor, applied for the enhancement of the coagulation
rate due to van der Waals forces [60]. After substitution of the Cunningham
11
smpV
Real particle system Stochastic system
1 3 4
5 6 8
2
7
V
Figure 2.1: Illustration of the stochastic solution methodology for the pop
ulation balance equations.
slip correction factor and the Knudsen number, the slip flow kernel is given
by [70]:
Ksf(Xi, Xj) =
2 kB T
3µ
(
1
di
+
1
dj
+ 1.257(2σ)
[
1
d2i
+
1
d2j
])
(di + dj) (2.24)
where µ is the gas viscosity and σ is the meanfreepath of gas molecules.
2.2.2 Stochastic solution of the population balance equations
This thesis is primarily concerned with using stochastic methods to solve the
population balance equations, however other methods are briefly discussed
in Section 2.2.3. Stochastic methods approximate the real particles with a
collection of computational particles in a sample volume Vsmp.
In the illustration of Figure 2.1, an ensemble of eight computational
particles is used to approximate the real system. Here, N(t) will be used to
denote the number of computational particles in an ensemble, with the max
imum capacity given the symbol Nmax. Based on the properties of these par
ticles, the stochastic particle algorithm calculates the particle process rates
described by Equations (2.15), (2.17), (2.19) and (2.21).
A process is selected to be performed, based on the rates. A particle
is then selected as the object of the process and is modified according to
the physiochemical reaction of the process, defined by Equations (2.14),
(2.16), (2.18) and (2.20). The process is repeated until the final time (tf) is
reached. This algorithm is illustrated in Figure 2.2.
The stochastic population balance solver is written in C++ and named
Sweep. In previous work, the basic stochastic algorithm has been modified
12
Start t ← 0
Calculate exponentiallydistributed
waiting time δt
t ← t + δt
Select process with probability
rate of process
sum of all process rates
Select particles(s) to adjust
Perform selected process
Update the particle ensemble statistics
t < tfFinish
FALSE TRUE
Figure 2.2: Stochastic algorithm for solving the population balance equa
tion.
13
to increase its computational efficiency [37, 53, 120]. These improvements
are described in the following subsections.
Binary tree for ensemble statistics
The stochastic algorithm requires calculation of the total process rate in ev
ery timestep δt (Figure 2.2). This requires looping over each possible colli
sion pair for coagulation, resulting in an order N(t)2 computational expense
[53]. This problem is avoided by storing information about particles in a bi
nary tree, which reduces the order of expense to Nmax log(Nmax) [54]. Each
node of the tree stores a partial sum of the properties of a particle from the
nodes below it. Thus, when a particle is changed, only the ‘branch’ of the
tree on which the particle lies must be updated. It also enables rapid selec
tion of a particle based on individual structural properties, e.g. surface area
[54].
Majorant rates and fictitious jumps
The expressions typically used for the collision kernel K in the coagulation
process contain crossterms, that is, they cannot be written as linear sums of
particle properties (e.g.
√
a + b). However, as the binarytree precalculates
certain sums, the rate calculation could be accelerated if the coagulation
kernel was rewritten to take advantage of these sums. Eibeck and Wagner
[37] proposed majorant kernels are used to resolve this issue. A majorant
kernel Kˆ must be greater than or equal to K for all P ∈ Z and should be fast
and efficient to calculate. A majorant kernel for the freemolecular regime
was presented by Goodson and Kraft [53], and a majorant formulation for
the transition kernel was given by Patterson et al. [121]. A new majorant
kernel is presented in this work in Chapter 4.
To correct for the overestimation of the coagulation rate from use of a
majorant kernel, the concept of fictitious jumps is introduced. When two
particles are selected for coagulation, the true rate R (that calculated with
the full kernel K) and majorant rate Rˆ are calculated for the particle pair.
Fictitious jumps, those where the process is not actually performed for the
particle pair, will occur with probability 1 – R/Rˆ.
14
Linear process deferment
Linear processes are those which involve a singleparticle interaction, e.g. sur
face growth or sintering. It is often the case that the rate of these processes is
much greater than that of nonlinear processes such as coagulation [8, 92],
resulting in more computational effort being devoted to linear interactions.
Patterson et al. [120] presented the Linear Process Deferment Algorithm
(LPDA) to address this, in which the stochastic algorithm (Figure 2.2) does
not include linear processes in the main loop. The deferred processes are
performed in bulk on each particle at the end of an arbitrary timestep ∆t.
2.2.3 Other solution methodologies
As this thesis focuses on stochastic solution of the population balance equa
tions, the two most popular alternative methods for solution are mentioned
briefly here.
Moment methods
The method of moments [47, 48] describes the evolution of a population
of particles by the moments of the distribution. It is typically fast, however
the particle size distribution (PSD) can not be resolved from the results of
such calculations. Despite this, their speed often makes them the choice
for coupling with detailed computational fluid dynamics (CFD) calculations
[3, 89].
Sectional methods
The sectional methods split the distribution of particle properties into sec
tions [65], allowing resolution of PSD. The choice of sections is nontrivial,
and often they must be dynamically chosen [156]. These methods are of
ten more computationally expensive [165] than the moments method and
can only be extended to multiple internal dimensions with great difficulties
[61]. A type of sectional method will briefly be introduced in more detail in
Chapter 5.
15
2.3 Reactor models
A reactor model defines how the temperature, pressure and volume of a re
actor change with time. These factors are dependent on the assumptions
made about the reactor and the source terms provided by the gasphase
and particle solvers. The choice of reactor model depends on the system
which is to be modelled: two of the most common target modelling systems
for nanoparticle synthesis are flames and plugflow reactors (PFRs). As a
stochastic approach for solving the population balance equations considers
the transient evolution of a sample volume, it is common to use reactor mod
els which are compatible with such restrictions. A common simplification is
to take a Lagrangian view of particles (and gas) evolving along a streamline
of a reactor or flame [99]. In this approximation, the equations govern
ing the reactor are identical to a zero (spatial) dimensional batch reactor
[26, 143].
The equations governing such systems are now presented. The rate of
change in concentration of species k is given by:
dck
dt
= ω˙k(c, T) + g˙k(c, n, T) – γ(c, n, T)ck (2.25)
where g˙k(c, n, T) is the molar production rates of species k at temperature
T due to particle reactions. The production rate due to gasphase reactions
ω˙k(c, T) is given by Equation (2.2). The gasphase expansion factor γ for
the reactor is:
γ(c, n, T) =
1
ρ
NGP∑
k=1
[
ω˙k + g˙k
]+ 1
T
dT
dt
(2.26)
where ρ is the molar density of the mixture. Ignoring energy changes due to
the particle population balance, the adiabatic energy balance of a constant
pressure reactor is:
ρCP
dT
dt
=
NGP∑
k=1
–Hˆkω˙k (2.27)
where CP is the bulk heat capacity and Hˆk is the specific molar enthalpy
of species k. The energy balance is analogous for a constantvolume reac
tor, however the specific molar internal energy is used in place of enthalpy.
These equations will be expanded in Chapter 5 for modelling networks of
reactors where transport of material between reactors can not be neglected.
16
stochastic population
balance solver
gas−particle
coupling
reaction rates
chemical
calculates reactor model
provides
ODE solver
Sweep
Sprog MOPS
Figure 2.3: Schematic of the architecture of MOPS.
2.3.1 Gasparticle coupling
Equations 2.15, 2.19, 2.25 and 2.26 imply a twoway relationship between
the particle and gasphase. Thus, it is appropriate to adopt a method of
coupling to link the two. Here, the technique of operator splitting is em
ployed to couple the gas and particlephases. In this work, Strang splitting
is the chosen form of operator splitting. The algorithm by which this is done
has been described in detail by Celnik et al. [25] and Shekar et al. [143].
The stochastic population balance solver (Sweep) is coupled with the
gasphase chemistry interface (Sprog and a separate ODE solver) in a code
called MOPS. The flow of information through these libraries is illustrated
in Figure 2.3. Implementation of the reactor models and network infras
tructure developed for the present work (Chapter 5) is also administered by
MOPS.
2.4 Error measurements
Chapters 4 and 5 investigate the convergence of the stochastic particle method’s
solution as a function of numerical parameters. The stochastic particle
method has three parameters which affect its numerical error, namely the
maximum number of computational particles (Nmax), the number of inde
pendent runs (L) and the splitting timestep size (∆ts) [25]. This numerical
error is composed of two parts: the systematic and statistical error (stat).
17
The former is related to the accuracy of the stochastic solution. It is
associated with choice of numerical parameters, for example too few com
putational particles. The statistical error is representative of the fluctuations
due to the psuedorandom choices made in selecting particles and processes
in the stochastic algorithm.
The systematic and statistical error can be assessed by generating L in
dependent estimates of the particle system, and comparing the macroscopic
quantities of the system (X(t)) for a given set of numerical parameters. The
temporal evolution of these functionals is averaged over the number of in
dependent runs:
ξ(Nmax,L)(t) =
1
L
L∑
l=1
X(Nmax)l (t) (2.28)
The confidence interval IP for ξ(Nmax,L)(t) is given by:
I(Nmax,L)P (t) = aP
√
Var(Nmax,L)(t)
L
(2.29)
where P is the confidence level. In this work, P = 0.999 is used, for which
a0.999 = 3.29 from normal distribution tables [143]. The variance is ap
proximated by:
Var(Nmax,L)(t) =
1
L – 1
L∑
l–1
X(Nmax)l (t)
2
– [ξ(Nmax,L)(t)]2 (2.30)
The total error is then estimated by calculating the normalised average ab
solute error () over M time intervals [25]:
(Nmax,L) =
1
M
M∑
m=1
∣∣∣ξ(Nmax,L)(tm) – X*(tm)∣∣∣
X*(tm)
(2.31)
where X*(t) is the ‘true solution’, which could be that found analytically.
Where analytical solutions are unavailable, a highprecision stochastic so
lution is used [143]. The normalised statistical error can be estimated by
taking the ratio of the confidence interval and the mean:
(Nmax,L)stat (t) =
I(Nmax,L)p (t)
ξ(Nmax,L)(t)
(2.32)
The normalised average statistical error may then be obtained:
(Nmax,L)stat =
1
M
M∑
m=1
(Nmax,L)stat (tm) (2.33)
18
X(
t)
time
syst. err.
stat. err.
ξ(t)
ξ(t)±I0.999(t)
X*(t)
Figure 2.4: Temporal evolution of macroscopic quantity X(t), showing sta
tistical and systematic error.
These measurements will provide a basis on which to investigate the nu
merical properties of the stochastic algorithms developed in this work. The
influence of these errors on a macroscopic property of the particle ensemble
X(t) (e.g. average diameter) is illustrated in Figure 2.4.
2.5 Silicon nanoparticles
Knowledge of the properties of silicon has underpinned the revolution of in
formation technology. A great deal of study has therefore been undertaken
into its synthesis, purification and resultant properties. Silicon nanoparti
cles were first made in the late 1970’s [101]. Since then, they have been
the subject of considerable work towards improving their manufacture and
identifying potential applications.
2.5.1 Manufacture
Gasphase, laser and plasma synthesis of particles are the most common
methods with which silicon nanoparticles are manufactured [88]. In gen
eral, these processes begin with silane (SiH4), which may be decomposed
19
Figure 2.5: Highlyspherical silicon nanoparticles made in a hotwall reac
tor. Reprinted from Journal of Aerosol Science, 41/11, Körmer,
R. et al., “Aerosol synthesis of silicon nanoparticles with narrow
size distribution–Part 1: Experimental investigations”, pp. 998–
1007, Copyright (2010), with permission from Elsevier.
by thermal, laser or microwave radiation sources [22, 72, 103]. The decom
position of silane forms reactive silicon hydrides, which combine with each
other to nucleate into silicon nanoparticles [149].
A tubular hotwall reactor is frequently used to study silicon nanoparti
cles in a research environment [75, 158]. These reactors allow control of
the temperature profile through the length of the reactor [161], and par
ticles can be directly sampled from within the reactor [75]. The particle
size, distribution and morphology is strongly dependent on the process con
ditions of the reactor, however large losses of silicon to wallgrowth is often
reported [75, 113]. Despite this, the product properties can be tightly con
trolled [56], and highly spherical particles, such as those in Figure 2.5, can
be obtained.
Production of silicon nanoparticles through laser pyrolysis typically uses
a CO2 laser to promote silane decomposition [22, 44]. Laser heating causes
rapid localised heating of the silane, leading to fast particle production. The
short timescales over which particles are produced makes controlling the
particle size distribution challenging. Laser ablation of solid silicon has also
been identified as a means to produce nanoparticles [157].
20
In plasma synthesis, silane is fed into a parallelplate reactor, where ra
diofrequency discharge ionises the precursor and initiates decay. A great
deal of study has focused on characterising the mechanisms of silane anion
formation, particle nucleation, and silicon deposition [50, 66]. Microwave
reactors have also been used to generate silane plasmas [72]. The major ad
vantage in plasma synthesis is that the particles become negatively charged,
which suppresses coagulation and confines them in the plasma [88].
Once particles are produced, they can be modified in many ways. Bring
ing silicon particles into contact with air will cause oxidation of the particle
surface, removing the SiH bonds at the surface. Particles can be made sta
ble in water by functionalising the surface with acrylic acid [86], or be made
stable in air by replacing the hydride surface layer with hydroxide [88]. Fi
nally, particle sizes can be reduced by etching with hydrogen fluoride and
nitric acid [85].
2.5.2 Applications
There are two key properties of silicon nanoparticles which make them in
teresting to researchers. The first is photoluminescence, the emission of light
due to the absorption of a photon. Early studies in silicon nanostructures
observed photoluminescence of visible wavelength after excitation with UV
light. It has since been confirmed to also occur through electrical excitation
[88]. As such, silicon nanoparticles could make promising candidates for
lightemission components. Cheng et al. [29] suggested they could be intro
duced into organic lightemitting diode (OLED) screens, where they could
be used in future touchscreen technology.
Photoluminescent silicon nanoparticles have also been considered for
medical imaging applications. Erogbogbo et al. [40] developed biocompati
ble luminescent silicon nanoparticles and successfully tracked them through
pancreatic cancer cells. Park et al. [115] reported that silicon quantum dots
exist in the body for a short time in comparison to other similar materials
(e.g. gold nanoparticles or carbon nanotubes), with no evidence of toxicity.
The other property of silicon nanoparticles is its use as a semiconduc
tor. Much study into applications of silicon nanoparticles in photovoltaics
has been conducted. As the quantum confinement effects of a particle are
strongly dependent on its size and structure, it may be possible to tune
a photovoltaic cell to better match solar emission than in bulk material
21
[88]. Finally, printable semiconductor and photoluminescent inks of silicon
nanoparticles have been proposed as a straightforward way of preparing
nextgeneration circuits or displays [7].
22
Chapter 3
Choosing a particle model
Three particle models commonlyused in the literature are com
pared in a detailed numerical study. The conditions under which
these models differ are identified, and the effect of particle polydis
persity upon each of the models is investigated. The study concludes
with an analysis of the computational time required for the models.
Material from this chapter has been published in [91].
The population balance equation (Equation (2.13)) describes the for
mation and growth of nanoparticles. When solving this equation, different
assumptions can be made about the nature of the particles, and it is this
level of detail which is defined as the ‘particle model’. Mathematically, the
particle model is represented by the typespace.
The aerosol system to be modelled and the choice of population balance
equation solution methodology define which particle model is used. The
sectional methods and moment methods discussed in Sections 2.2.3 are dif
ficult to extend to multiple internal dimensions [61] and will therefore use
relatively simple particle models [47]. The stochastic particle method used
in this thesis is ideal for comparing particle models, as the particle model is
completely independent of the stochastic algorithm, and a typespace of any
level of detail can be used.
Some of the most commonlyused particle models will now be reviewed.
In early modelling attempts, particles were treated as spheres which would
coalesce into larger spheres upon collision with other particles; for example,
see Frenklach and Harris [48]. This corresponds to the spherical particle
23
model; which is a onedimensional model tracking only the volume or some
number of chemical elements in a particle.
However, not all nanoparticles are spherical [117, 139, 152]. It is very
common to obtain aggregates (chemically bound collections of particles)
and agglomerates (physically bound particles) [36], especially when parti
cle growth occurs by coagulation or sintering. The surfacevolume (surfvol)
particle model extends the basic spherical model by tracking the volume and
total surface area of a particle [78, 165]. Another twodimensional model
providing similar information tracked the number of aggregates and the
number of primary particles [116]. Such knowledge provides a basic un
derstanding of aggregate structure of particles. The surfacevolume model
was also extended in the surfacevolumehydrogen model of Blanquart and
Pitsch [15], which considered the number of active sites as an additional
independent variable.
The surfacevolume model is limited in the respect that it assumes that
primary particles are monodisperse. Further, it has been identified that typ
ical expressions used to calculate the primary particle size are not ideal
[155]. Heine and Pratsinis [61] developed a sectional model to address
this issue, where the polydispersity of primary particles was incorporated
in a sectional model. Balancing the population of primary particles over
all aggregate particles enabled the simultaneous effects of coagulation and
sintering to be captured.
Even more detailed particle models have been developed. Sander et al.
[136] presented a binarytree particle model as part of a stochastic popu
lation balance model. It incorporated the full aggregate structure of soot
particles, allowing for resolution of individual primary particles and their
connections to other primary particles in an aggregate. This model has been
extended to silica [143, 144] and silicon [92] in later studies.
Despite these advances in modelling the structure and morphology, many
studies still use basic spherical or surfacevolume particle models [4, 56,
76]. In this chapter, the conditions under which use of these models is
appropriate will be investigated, and the conditions where their use could
represent a significant oversimplification will be identified. This is accom
plished using the generic stochastic population balance solver presented in
Section 2.2.2.
The structure of this chapter is as follows. Section 3.1 presents the math
24
ematical description of the particle models and explains how their derived
properties (e.g. particle collision diameter) are calculated. Case study are
presented in Section 3.2 and the implications of these studies discussed in
Section 3.3.
3.1 Model
3.1.1 Typespace
This work uses a generic multicomponent description of the primary parti
cles’ chemical composition, that is:
p = p (η1, η2, . . . , ηn) (3.1)
where p will contain ηi chemical components (e.g. number of carbon atoms).
A particle Pq of index q in Nmax computational particles will always contain
at least one primary particle p. Throughout this work, the term ‘particle’ will
be used to refer to the computational particle P, whereas ‘primary’ will refer
to a primary particle owned by a computational particle. For the spherical
particle model, the typespace is very simple:
Pq = Pq
(
p(q)
)
(3.2)
Here the superscript is used to denote ownership of the primary particle
p(q) by computational particle Pq. The surfacevolume model tracks the
surface area S(q) in addition to the chemical composition of primaries:
Pq = Pq
(
p(q), S(q)
)
(3.3)
The surfacevolume model, however, assumes that primaries in a particle
Pq are uniformly distributed. The ability to include polydispersity of primary
particles was proposed by Sander et al. [136] in their binarytree model.
This model represents a particle as a binary tree of primary particles with
connectivity stored in the nonleaf nodes:
Pq = Pq
(
p1, p2, . . . , pnq , C
(q)
)
(3.4)
where Pq contains nq primary particles and C is a lowerdiagonal matrix
representing the common surface area between two neighbouring primary
particles. The form of this matrix is discussed in [136, 143, 144].
25
Table 3.1: Comparison of the derived properties of the particle models.
Property Spherical Surfvol Binary tree
typespace Pq
(
p(q)
)
Pq
(
p(q), S(q)
)
Pq
(
p1, p2, . . . , pnq , C(q)
)
Vq
n∑
i=1
ηiMi
NAρi
n∑
i=1
ηiMi
NAρi
nq∑
j=1
V(pj)
dsph,q
3
√
6
pi
Vq
3
√
6
pi
Vq
3
√
6
pi
Vq
Ssph,q pid
2
sph,q pid
2
sph,q pid
2
sph,q
Sq Ssph,q S
(q) Ssph,q
sq(1 – n
– 13
q ) + n
– 13
q
dpri,q dsph,q
6Vq
Sq

dpriq  
1
nq
nq∑
j=1
dpri(pj)
nq 1
S3q
36piV2q
nq
dcol,q dsph,q
1
2
(
dsph,q +
√
Sq
pi
)
dpriq
[
S3q
36pi V2q
] 1
Df
3.1.2 Derived properties
The typespace of the particles allows for all other properties of the particles
to be determined [143]. The volume of a primary is based on the number
of chemical units ηi and bulk densities ρi:
v(p) =
n∑
i=1
ηiMi
NAρi
(3.5)
where Mi is the molecular weight of element i and NA is Avogadro’s
number. For spherical or surfacevolume particles, the volume of the particle
Vq is equal to that of the primary. In the binary tree model, the volume of a
particle is the sum of the volume of its constituent primaries, that is:
Vq = V(Pq) =
nq∑
j=1
v(pj) (3.6)
The derived properties of the surfacevolume model have been published
in [78]. Essentially, the volume and surface area Sq of a particle are used
26
to estimate the size of primary particles, the number of primaries and the
collision diameter. These equations are summarised in Table 3.1. For the
binary tree model, the derived properties are analogous to those used in
previous studies [90, 136, 143]. Of key importance is the surface area:
Sq =
Ssph,q
sq(1 – n
– 13
q ) + n
– 13
q
(3.7)
where Ssph,q is the equivalent spherical surface of a sphere with the same
volume as the particle Pq, and sq is the average sintering level between the
primaries of the particle. The sintering level (0 ≤ s ≤ 1) of two primary
particles pj and pk is given by [136]:
s(pj, pk) =
Ssph(pj,pk)
Cjk
– 2–
1
3
1 – 2–
1
3
(3.8)
where Cjk represents the element of the commonsurface matrix C describ
ing the common surface area of the two primaries pj and pk. For two primary
particles in point contact, this initially represents the sum of their surface
area. As sintering occurs, the common surface decreases until it is equal
to the surface area of the sphere with the same volume as the sum of the
primaries’ volume, Ssph(pj, pk). Thus, we can always write:
Ssph(pj, pk) < Cjk ≤ Ssph(pj) + Ssph(pk) (3.9)
The definition of sintering level implies that a spherical particle has a
sintering level of 1, while two primaries in point contact (with no sinter
ing) have a sintering level of 0. Note that the primary particle diameter
described by the surfacevolume model is an estimate basedon the surface
area and volume. The primary particle diameter for the binarytree model
is an average within the particle, given by:
dpriq =
1
nq
nq∑
j=1
dpri(pj) (3.10)
The collision diameter of a particle in the binary tree model is given by
[143]:
dcol,q = dpriq
[
S3q
36pi V2q
] 1
Df
(3.11)
27
where dpriq is the average primary diameter (as calculated by Equation
(3.10)) and Df is the fractal dimension, assumed to be 1.8 for the present
work [138]. These are summarised and compared with the derived proper
ties of other particle models in Table 3.1.
3.1.3 Particle processes
The three particle models interact with several particle processes. Exam
ples of such processes are inception (or nucleation), heterogeneous growth,
coagulation and sintering. These processes are documented and mathemat
ically described in detail in [90]. In the present work, coagulation and sin
tering are the two which differ most between models so are outlined here.
The rate of coagulation is calculated using the transition regime coag
ulation kernel (Equation (2.22)), which is a computationally efficient ap
proximation to true Brownian coagulation [70]. In the spherical particle
model, a coagulation event will form a sphere with volume equal to the sum
of the coagulating spheres’ volume. This is represented by the following
transformation:
Pq
(
p(q)
)
+ Pr
(
p(r)
)
→ Ps
(
p(q) + p(r)
)
(3.12)
The addition of primary particles is equivalent to taking the sum of their
component vectors:
p(q)(η(q)1 , . . ., η
(q)
n ) + p
(r)(η(r)1 , . . . , η
(r)
n )
= p(s)(η(q)1 + η
(r)
1 , . . . , η
(q)
n + η
(r)
n )
(3.13)
A similar transformation takes place for a coagulation event in the surface
volume model, except the total surface area is also conserved:
Pq(p(q), S(q)) + Pr(p(r), S(r))→ Ps(p(q) + p(r), S(q) + S(r)) (3.14)
A coagulation event in the binary tree model will retain the primary
particle information and connectivity of both coagulating particles:
Pq
(
p1, . . . , pnq , C
(q)
)
+ Pr
(
p1, . . . , pnr , C
(r)
)
→ Ps
(
p1, . . . , pnq , pnq+1, . . . , pnq+nr , C
(s)
) (3.15)
The commonsurface matrix C is changed to reflect the structure of old
particles Pq and Pr, and their new connection. The process through which
28
p(1) + p(2)
p(1) + p(2)
P(p(1) + p(2)) P(p
(1) + p(2),S(1) + S(2))
p1 p2
C12
P(p1, p2,C)
spherical surfvol binarytree
Figure 3.1: An illustration of different particle models’ description of coag
ulation events
this is done is described in detail by Shekar et al. [143]. The physical inter
pretation of Equations (3.12)–(3.15) is illustrated in Figure 3.1.
The three models also include different descriptions of sintering. For
the spherical particle model, sintering effectively occurs infinitely fast, as
particles coalesce upon coagulation. For the surfacevolume and binary tree
models, the exponential excess surface area decay formula–as popularised
by Koch and Friedlander [73]–is employed. In surfacevolume model, the
surface area of the whole particle is reduced:
∆Sq
∆t
= –
1
τS(p(q))
(
Sq – Ssph,q
)
(3.16)
where τS is the characteristic sintering time. It is an empirical function
of temperature T and diameter d which is related to the time required for
two neighbouring primaries to coalesce. For the surfacevolume model, it is
expressed in the form [155]:
τS(p
(q)) = AS d
nS
pri,q T exp
(
ES
R T
)
(3.17)
where AS, nS and ES are empirical parameters. These parameters are
usually specific to a particular type of material and are not trivial to deter
mine [19, 166]. Tracking of the connectivity of adjacent primary particles
allows the binary tree model to sinter each connection individually. For a
29
single connection, the rate of sintering is given by [136]:
∆Cjk
∆t
= –
1
τS(pj, pk)
(
Cjk – Ssph(pj, pk)
)
(3.18)
where Ssph(pj, pk) is the equivalent spherical surface area of a particle with
volume Vj + Vk. The characteristic sintering time also takes a slightly differ
ent form to account for two primaries of a different size touching:
τS(pj, pk) = AS min
[
dpri(pj), dpri(pk)
]nS T exp( ES
R T
)
. (3.19)
3.1.4 Average properties
In order to study the behaviour of a system of particles, it is useful to de
fine some average quantities which describe the size and morphology of the
particles and ensemble. In this thesis, the arithmetic and geometric mean
of the collision and primary diameter is used. These are denoted as µa and
µg respectively. For example, the geometric mean of the primary diameter
is given by:
µg(dpri) =
N(t)∏
q=1
dpri,q
1
N(t)
(3.20)
where N(t) computational particles are in the ensemble at time t. Note that
for the surfacevolume model, µg(dpri) represents the mean of the estimated
primary diameter of each particle; while in the binary tree model it repre
sents the mean of average primary diameter of each particle.
The (numberbased) geometric standard deviation is also particularly
useful for classifying polydispersity of particles [61]. Here, it is calculated
in the context of the binary tree model for the average primary diameter as:
σg(dpri) = exp
√√√√√ 1
N(t)
N(t)∑
q=1
ln
dpri,q
µg(dpriq)
2
(3.21)
For the binary tree model, two estimates of σg(dpri) can be made. The
first involves using Equation (3.21) to calculate the geometric standard de
viation of the average primary diameter. The second takes the arithmetic
mean of the geometric standard deviation of the primary diameter in each
30
particle:
µa
[
σg(dpri)
]
=
1
N(t)
N(t)∑
q=1
exp
√√√√√ 1
nq
nq∑
j=1
ln
[
dpri,j
µg(dpri)
]2 (3.22)
These two estimates provide different descriptions of the geometric stan
dard deviation of primary diameter. σg(dpri) refers to the average geomet
ric standard deviation of primary diameter across all particles. In contrast,
µa
[
σg(dpri)
]
yields information about the geometric standard deviation of
primaries contained within a particle.
It should be noted that the average properties presented here are dis
tinct from those in Section 2.4, as the former represent averages over the
computational particles, and the latter over independent runs.
3.2 Case studies
In order to study the structure of particles simultaneously undergoing coag
ulation and sintering, several properties are first defined. The characteristic
coagulation time τC is given by [103, 168]:
τC =
1
KtrM0
(3.23)
where Ktr is the transition regime coagulation kernel evaluated for the par
ticle ensemble under study (m3/s) and M0 is the zeroth moment, or particle
number density (1/m3). For the cases in Sections 3.2.1 and 3.2.2, the ini
tial characteristic coagulation time τC,i is used. The characteristic times are
normalised by the residence time τ , yielding a dimensionless characteristic
time.
In terms of numerical parameters, 16,384 computational particles (Nmax)
and 8 independent runs (L) were used to evaluate the following studies
[143]. All results presented in the following sections are runaverages (Equa
tion 2.28), with the ξ symbol neglected for clarity.
3.2.1 Coagulation and sintering
The three models presented in the present work use very different repre
sentations of coagulation and sintering processes. To compare the models,
31
a sample ensemble of monodisperse spherical particles of diameter di was
generated. The initial characteristic coagulation time τC,i of these particles
was estimated basedon their size and number density; and the sintering
time τS was set to a constant value. The space of normalised coagulation
time and normalised sintering time was scanned, and the resulting struc
tures of particles analysed. Each point in this space represents a calculation
for the same initialised ensemble of particles where only coagulation and
sintering could operate at a specified rate.
The results of this analysis are presented in Figure 3.2, where the colli
sion diameter (normalised by initial particle diameter) is plotted against the
dimensionless sintering and initial coagulation time. The collision diameter
was chosen as it is a strong function of typespace and is a good descriptor
of overall aggregate size.
It is unsurprising that the models converge for high values of τ/τS.
This corresponds to instantaneous sintering (or coalescence) of particles
upon collision. It is interesting to observe that for some combinations of
(τ/τS, τ/τC,i) the particles obtained through the binary tree model are larger;
whereas for other combinations they appear smaller than the surfacevolume
model. Indeed, an intersection between these surfaces defines a line at
which the particles obtained through both models are identical. This is
shown in more detail in Figure 3.3, where the contours corresponding to
the intersection of the binary tree and surfacevolume models, and ±5%
relative difference between the two are plotted. Specifically, the 5% error
shading refers to the conditions under which:
0.95 ≤ µg(dcol)
bintree – µg(dcol)
surfvol
µg(dcol)bintree
≤ 1.05 (3.24)
Four zones in this figure are identified, corresponding to different types
of particles. Particles obtained in the spherical zone sinter sufficiently quickly
to coalesce into spherical particles. In this zone, the all particle models effec
tively reduce to the spherical particle model; and there is no benefit gained
from tracking the additional structural detail of particles.
In the nearspherical zone, aggregates containing few primary particles
are formed, either due to a fast sintering rate or a slow coagulation rate.
Figure 3.2 illustrates that the particles predicted by the binary tree model
have a larger collision diameter, corresponding to smaller primary particles.
The reason for this can be understood by considering the derived properties
32
Figure 3.2: Comparison of the normalised geometric mean collision diam
eter as a function of dimensionless coagulation and sintering
time for each of the three particle models. The thick line depicts
the intersection of the binary tree and surfacevolume model
predictions.
33
τ/τ
C,
i,

τ/τS, 
intersection
5 % error
102
101
100
101
102
105 104 103 102 101 100 101 102 103
Partialsintering zone
Nearspherical zone
Spherical
zone
Aggregate
zone
Figure 3.3: A contour plot depicting where the binary tree model is equiv
alent to the surfacevolume model. The shading corresponds to
within 5% deviation (with respect to collision diameter) of the
surfacevolume model from the binary tree model.
34
of the particle models, in particular the expressions for collision diameter. It
is shown in Appendix A.1 that, for coagulation of similarlysized particles,
the surfacevolume model will underestimate the predictions of the binary
tree model to a maximum of 10% difference. While this deviation is small,
it causes potentially large changes in the rate of coagulation (rate ∝ d2col in
the freemolecular regime), thus increasing the maximum error.
In the partialsintering zone, particles coagulate to form partiallysintered
aggregates with more than 2–3 primary particles. Here, the surfacevolume
model overpredicts the results of the binary tree model. The reason for
this can be understood by considering the sinterable surface area of an ag
gregate. As each connection node in the binary tree model is individually
sintered (according to Equation (3.18)) there must always be more sinter
able surface area in an aggregate described by the binary tree model than
one described by the surfacevolume model. This is mathematically demon
strated in Appendix A.2. If sintering occurs within the same characteristic
time, the model with more sinterable surface area will sinter the particles
faster (Equations (3.16) and (3.18)). In this case, this causes the binary tree
model to predict the formation of smaller aggregates with larger primaries.
The intersection of the partialsintering and nearspherical zones corre
sponds to the point where primaries are uniformly sized within an aggre
gate. As τ/τS → 0 the binary tree model becomes identical to the surface
volume model: aggregates are composed of unsintered monodisperse pri
mary particles (the aggregate zone).
3.2.2 Effect of polydispersity
The study in the previous section highlights that under certain conditions,
the surfacevolume model yields aggregate particles with approximately the
same size as the binary tree model. This, however, was determined for
an initially monodisperse ensemble of particles (initial geometric standard
deviation σg(di) = 1.0). The purpose of this case study is to investigate
how polydispersity can affect model performance, under conditions where
the surfacevolume and binarytree models are equivalent for monodisperse
particles. To investigate this, a point on the intersection line of Figure 3.3
was chosen (τ/τS = 1.4, τ/τC,i = 5.0) and ensembles initialised with spher
ical particles of increasing polydispersity (up to σg(di) = 3.0). The results
are given in Figure 3.4.
35
0
5
10
15
20
25
30
1 1.5 2 2.5 3
µ g
(d c
o
l) /
µ g
(d i
), 
σg(di), 
bintree
surfvol
spherical
2
4
6
8
10
12
14
16
18
20
1 1.5 2 2.5 3
µ g
(d p
ri)
/ µ
g(d
i),

σg(di), 
1
1.2
1.4
1.6
1.8
2
1 1.5 2 2.5 3
σ
g(d
co
l),

σg(di), 
1
1.2
1.4
1.6
1.8
2
1 1.5 2 2.5 3
σ
g(d
pr
i),

σg(di), 
bintree µa[σg(dpri)]
Figure 3.4: Effect of initial polydispersity, σg(di), on morphology of parti
cles. Shown are the normalised geometric mean µg and stan
dard deviation σg of the collision and primary diameter. These
were generated for τ/τS = 1.4, τ/τC,i = 5.0, on the intersection
line of Figure 3.3.
36
It is evident that the the surfacevolume and binarytree models perform
comparably with σg(di) < 1.5. Under these conditions, there appears to be
good agreement with the geometric mean collision and primary diameter
and a small difference in geometric standard deviation of both. Beyond
σg(di) = 1.5, the binarytree model appears to predict larger aggregates
with smaller primaries than the surfacevolume model. This is attributed to
the phenomenon discussed in Section 3.2.1 for the prediction of collision
diameter in the nearspherical zone.
It is also interesting to observe that for the binary tree model, the stan
dard deviation of the average primary diameter is significantly smaller than
the average standard deviation of the primary diameter. This suggests that
the local polydispersity within a particle can often be quite different to the
polydispersity of the ensemble. The right panels of Figure 3.4 show that the
spherical particle model yields ensembles with σg(d) ≈ 1.4 for σg(di) > 1.5.
This indicates that the system has reached the selfpreserving particle size
distribution [153].
In summary, this case study demonstrates that the prediction of the prop
erties of an ensemble of polydisperse particles undergoing coagulation and
sintering is strongly dependent on choice of particle model. Highly polydis
perse σg(di) > 1.5 systems are poorly described by the spherical or surface
volume models where coagulation and sintering occur on similar timescales.
3.2.3 Computational resources
While the binary tree model is capable of providing more information on
particle structure and morphology, this comes at the cost of increased com
putational resources required. In order to capture all particle structure
zones, a ‘slice’ across Figure 3.3 at τ/τC,i = 10 is taken, and the com
putational time required per run plotted in Figure 3.5. Calculations were
conducted on a SGI Altrix Cluster with nodes composed of two 3.00 GHz
quadcore Intel Harpertown processors and 8 GB of RAM.
In the stochastic formulation of the population balance, it is evident that
the surfacevolume and spherical particle models require almost the same
computational time, and this time is independent of the sintering rate. This
is unsurprising, given that the additional effort required to adjust a particle’s
surface area after a sintering event is minimal: it merely requires calculation
of Equation (3.16).
37
100
101
102
103
104
105 104 103 102 101 100 101 102
t C
PU
/L
, s
τ/τS, 
Partialsintering zone
N
ea
rs
ph
er
ica
l
zo
n
e
Sp
he
ric
al
z
on
e
Aggregate zone
bintree
surfvol
spherical
Figure 3.5: Comparison of computational time (tCPU) required per run for
τ/τC,i = 10 to calculate τ = 0.5 s of process time.
The three models require the same computational time when spherical
particles are obtained. Only when coagulation becomes important (relative
to sintering) does the binarytree model require more time, and in this case
it requires orders of magnitude more. In the partialsintering zone, par
ticles containing more than one primary experience ‘merge’ events, where
two distinct primaries are combined into one due to sintering [143]. These
events alter the binary tree structure and as such, for large trees (i.e. par
ticles with many primaries), more time is required to regenerate the tree
structure after a node is lost. This phenomenon is responsible for the peak
in computational time in the partialsintering zone.
3.3 Conclusions
This chapter has presented the mathematical formulation of particle models
commonly used in aerosol dynamics of nanoparticles. A detailed numerical
study was conducted in order to understand under which conditions these
models differ. Starting from a set of identical conditions, the predictions
38
of the resulting average properties of the particle ensemble were compared
for each of the particle models. It was identified that under certain circum
stances, all models were equivalent to a spherical particle model.
Outside of this range, the performance of the popular and computation
ally inexpensive surfacevolume model with respect to a model which tracks
full aggregate structure (binarytree model) was evaluated. The two mod
els gave similar results (to within a small error margin) where sintering
was slow. However, where coagulation and sintering occurred on a similar
timescale, the models predicted substantially different results. Further, the
spherical and surfacevolume particle models were shown to incur a large
error margin where the coagulating ensemble was highly polydisperse.
There is still considerable scope for further investigation of the behaviour
of particle models. Heterogeneous growth processes are welldocumented
contributors to rounding of particles [83, 144] and have not been explicitly
addressed in the present work. Further, only the most popular version of the
surfacevolume model has been used, when alternative twodimensional for
mulations [61, 116] and different expressions for derived properties [155]
exist.
As many population balance models use some degree of fitting or pa
rameter estimation [27, 92, 142], the error in use of a ‘simple’ model will
be reflected in the parameters obtained through the fitting procedure. The
present work highlights that the choice of particle model does matter, and
that a target modelling system should be wellcharacterised experimentally
before proceeding to modelling.
39
Chapter 4
A new coagulation kernel for
stochastic weighted algorithms
The transition regime coagulation kernel is adapted for use with
stochastic weighted algorithms (SWAs). A numerical study into use
of the new kernel with a multidimensional silica particle model is
conducted, and the convergence properties of SWAs are evaluated.
Material from this chapter has been published in [93].
The key challenge in solving the population balance equations with a
stochastic particle method is accounting for the coagulation mechanism.
The Direct Simulation Algorithm (DSA)–introduced in [37]–was developed
to solve the Smoluchowski coagulation equation through the technique of
fictitious jumps and majorant kernels.
In direct simulation, every computational particle in the ensemble rep
resents the same number of real particles. In basic implementations of DSA,
coagulation events deplete the ensemble until there is only one computa
tional particle remaining in the ensemble. To avoid this, the ensemble may
be duplicated when it is below 50 % capacity [81]. Even when extended
to avoid computational particle depletion, DSA produces statistically noisy
estimates of the concentrations of physically rare particles and other asso
ciated quantities [34, 122]. This is usually overcome by computationally
expensive changes in the numerical parameters, for example, by increasing
the number of computational particles in the ensemble [122].
To avoid these issues, stochastic weighted algorithms can be used [38,
40
74, 122, 134, 164]. These methods have been presented in various forms
such as the Mass Flow Algorithm (MFA) [38, 54, 98, 99, 155] and the algo
rithms developed in [59, 170]. Stochastic weighted algorithms differ from
direct simulation and constantnumber methods in that coagulation events
do not reduce the number of computational particles. Instead, simulated co
agulations in the SWAs adjust the statistical weight that is attached to each
computational particle. This statistical weight is proportional to the number
of physical particles represented by the computational particle.
A recent paper gave a comprehensive review of present adaptations of
SWAs [122]. It also gave the general formulation of weighted particle meth
ods. This framework allowed for the development of general weighting rules
which describe how the statistical weights are adjusted during coagulation.
Further, it assumed a general form of the coagulation kernel and formulated
this in terms of majorant techniques: an important method of reducing com
putational costs [37, 121].
The majority of previous studies of SWAs presented either used a simple
spherical particle model [33, 34, 133] or a simple coagulation kernel for
which an analytical solution to the population balance can be found [122,
170]. In one case, it was suggested that discrepancies with experimental
results were a product of not using a coagulation kernel valid across the full
range of Knudsen numbers [98].
In the previous chapter, the importance of using detailed particle models–
especially where coagulation and sintering operate on similar timescales–
was highlighted. Further, it has been previously identified that SWAs are
particularly useful when extended to spatially inhomogeneous systems [57,
69, 118]. Thus, if stochastic methods are to be used to model systems with
detailed typespace particle models, the numerical behaviour of these meth
ods must be investigated. This will be achieved by extending the numerical
studies on a silica particle model, conducted by Shekar et al. [143].
The structure of this chapter is as follows. In Section 4.1, the stochastic
particle method used to solve the population balance equation is outlined.
Section 4.2 shows how to implement the stochastic weighted algorithm for
the transition regime coagulation kernel, and Section 4.3 gives the numer
ical studies used to test this system. The efficiency and convergence of the
SWAs is discussed in Section 4.4 and the chapter is concluded with a critical
analysis of these methods.
41
4.1 Model
A fullycoupled gasphase and particle model is used to simulate the for
mation of silica from tetraethoxysilane. The kinetic model describing the
gasphase decomposition is described in [126, 142] and consists of 58 re
versible gasphase reactions among 27 chemical species.
4.1.1 Typespace
The binarytree model presented in Section 3.1 is adapted here to the par
ticular case of silica. A detailed mathematical formulation of this model’s
typespace is given in [143], however a brief overview of its salient features
is presented here. Primary particles p are represented as:
px = px(ηSi, ηO, ηOH) (4.1)
where ηSi is the number of silicon atoms, ηO is the number of oxygen (non
OH) atoms and ηOH is the number of OH units. Particles are represented
as:
Pq = Pq(p1, . . . , pnq , C
(q)) (4.2)
where particle Pq contains nq primary particles. As the particle model tracks
the number of chemical units of each primary and the aggregate structure,
it is able to provide detailed information about the nature of a single silica
particle.
4.1.2 Particle processes
In this model, particles may be created and changed through several differ
ent processes. The implementation of these processes in the framework of
the stochastic weighted algorithm is described in Section 4.2.2.
Inception
The collision of two gasphase Si(OH)4 monomers yields a particle
with an initial state given by p1(ηSi = 2, ηO = 1, ηOH = 6). This
occurs through the reaction:
Si(OH)4 + Si(OH)4
→ Pq
(
p1(ηSi = 2, ηO = 1, ηOH = 6), C
(q)
)
+ H2O
(4.3)
42
where the newlyincepted particle is denoted by P. Their rate of for
mation is based on the transition coagulation kernel:
Rincep =
1
2
Ktr
(
Si(OH)4, Si(OH)4
)
N2A c
2
Si(OH)4 (4.4)
where Ktr
(
Si(OH)4, Si(OH)4
)
is the transition kernel as applied to
Si(OH)4 monomers, NA is Avogadro’s constant and cSi(OH)4 is the gas
phase concentration of Si(OH)4.
Coagulation
The collision and sticking of two silica particles is referred to as coag
ulation. The coagulation of particles Pq and Pr is represented by the
following change to the particles’ state:
Pq
(
p1, . . . , pnq , C
(q)
)
+ Pr
(
p1, . . . , pnj , C
(r)
)
→ Ps
(
p1, . . . , pnq , pnq+1, . . . , pnq+nr , C
(s)
) (4.5)
The matrix C is transformed to include the connectivity of the old
particles and that of the new pointconnection [143]. The rate of co
agulation is calculated using the transition coagulation kernel.
Surface reaction
Surface reactions may occur when a Si(OH)4 monomer reacts with an
OH site on the surface of a primary particle in a fashion analogous to
inception. This alters the composition of the particle in the following
manner:
Pq
(
p1, . . . , px, . . . , pnq , C
(q)
)
+ Si(OH)4
→ Pr
(
p1, . . . , p
′
x, . . . , pnr , C
(r)
)
+ H2O
(4.6)
where nq = nr. A primary px is uniformly selected and transformed
according to:
px(ηSi, ηO, ηOH)→ p′x(ηSi + 1, ηO + 1, ηOH + 2) (4.7)
The common surface matrix (C) is changed by a surface reaction event
(C(q) → C(r)) on px according to the following equation:
C′xy =
0 px, py are nonneighbouring
Cxy +
2∆vσSR
dpri(px)
px, py are neighbours
(4.8)
43
where ∆v is the volume element resulting from addition of mass and
σSR is the surfacerounding factor, a smoothing factor such that 0 ≤
σ ≤ 2 [137]. Cxy represents the element of the commonsurface ma
trix C which describes the common surface area of the two primaries
px and py. The rate of surface reaction is proportional to the con
centration of Si(OH)4 precursor and the number of OH sites on the
particle surface:
Rsurf(Pq) = Asurf exp
(
–
EA
RT
)
ηOH(Pq) NA cSi(OH)4 (4.9)
where Asurf and EA are the gasphase Arrhenius constants.
Sintering
Neighbouring primary particles may sinter with each other according
to the excess surface area decay formula popularised by Koch & Fried
lander [73]:
∆Cxy
∆t
= –
1
τS(px, py)
(
Cxy – Ssph(px, py)
)
(4.10)
where Ssph(px, py) is equivalent surface area of primaries px and py.
Here, a viscousflow model is used to describe the characteristic sin
tering time τS as a function of diameter and temperature:
τS(px,py) =
AS min
[
dpri,x, dpri,y
]
exp
(
ES
T
[
1 –
dp,min
min
[
dpri,x, dpri,y
]]) (4.11)
where AS, ES and dp,min are empirical constants obtained from model
fitting [142]. The degree of sintering is measured by the sintering
level parameter [136]:
s(px, py) =
Ssph(px,py)
Cxy – 2
– 13
1 – 2–
1
3
. (4.12)
where Ssph(px, py) is the surface area of the sphere with volume equal
to the sum of the primaries’ volume.
When the sinter process is called, the sintering level, s, of the neigh
bouring primaries is first checked against the coalescence criterion. If
44
the sintering level is greater than 0.95, the particles are merged; oth
erwise, the particles are sintered according equation (4.10) for time
∆tsinter with
∆tsinter = t – tlast update (4.13)
where tlast update was the last time the sintering characteristics of the
particle were updated. Sintering also results in the release of water,
the amount of which, like surface reaction, is proportional to the par
ticle’s statistical weight.
Intraparticle reaction
The OH units on the surface of a primary particle may also react
with each other to release water. A particle is transformed by a intra
particle reaction as:
Pq
(
p1, . . . , px, . . . , pnq , C
(q)
)
→ Pr
(
p1, . . . , p
′
x, . . . , pnr , C
(r)
)
+ H2O
(4.14)
where nq = nr. C is altered in the same manner as for surface reaction.
A primary px is uniformly selected and adjusted by:
px(ηSi, ηO, ηOH)→ p′x(ηSi, ηO + 1, ηOH – 2) (4.15)
The rate of intraparticle reaction is calculated for each particle such
that the Si:O ratio tends to 1:2 as t → ∞. It is mathematically de
scribed by the difference between the surface reaction rate and the
wholeparticle sintering rate:
Rint(Pq) = Rsurf(Pq) –
ηOH(Pq)
2 S(Pq)
nq∑
x,y=1
Cxy – Ssph(px, py)
τS(px, py)
(4.16)
The reader is referred to [143] for the derivation and explanation of
this formula.
4.2 Development of stochastic weighted algorithms
In stochastic weighted algorithms, computational particles are described by
a type Pi and a statistical weight wi. The physical concentration of particles
of type Pi is given by wi/Vsmp.
45
4.2.1 Linear particle processes
While the rate expressions of linear particle processes such as inception and
surface reaction remain unchanged in this weighted particle method [122],
the interaction of these with the gasphase must be carefully considered. A
process which removes one molecule of χ1 from the gasphase and release
one molecule of χ2 into the gasphase can be written as:
(Pi, wi) + χ1 → (P′i, wi) + χ2 (4.17)
The amount of χ1 removed (∆cχ1) and χ2 added is proportional to the
statistical weight of the particle wi and the scaling factor Vsmp, that is:
∆cχ1 = –
wi
NA Vsmp
(4.18)
where Vsmp refers to the computational sample volume, and:
∆cχ2 =
wi
NA Vsmp
(4.19)
Note that the statistical weight wi remains unchanged throughout this pro
cess. Particles are created in the population balance with weight wi = 1.0
where inception occurs.
4.2.2 Coagulation process
Coagulation events in the SWA are asymmetric [38]. An ordered pair of
particles (Pi, wi) and (Pj, wj) coagulate at rate K(Pi, Pj)wj such that:
(Pi, wi), (Pj, wj)→ (Pi + Pj, γ(Pi, wi, Pj, wj)), (Pj, wj) (4.20)
The computational particle (Pj, wj) is unaffected by the simulated coagula
tion event, but the new particle with index i is constructed using (Pj, wj).
This applies to the new type Pi + Pj and the new statistical weight, which is
calculated with the weight transfer function γ [122].
In the stochastic simulation of a coagulation process, the total jump rate
R for coagulation using the direct simulation algorithm is given by [121]:
R =
1
2
N(t)∑
i6=j
K(Pi, Pj)
Vsmp
(4.21)
where N(t) is the number of stochastic particles in the ensemble at time
t. When a weighted particle method is used, the jump rate is given by
46
K(Pi, Pj)wj and coagulation is no longer symmetric [122], so the total rate
expression is modified to:
R =
N(t)∑
i 6=j
K(Pi, Pj)wj
Vsmp
(4.22)
The function K(Pi, Pj) represents the physics of the coagulation process
and is dependent on characteristics of the particle pair and the chemical
conditions of the surroundings. The transition kernel (Equation (2.22)) is
a popular choice as it approximates the more detailed Fuchs Interpolation
Formula and is valid across a wide range of Knudsen numbers [70, 127].
The implementation of the transition kernel within direct simulation is
explained in detail in [121]. Of central importance is the use of majorants
as they provide computationally efficient approximations of the true kernel
such that K ≤ Kˆ. In the context of SWAs, the majorant kernel Kˆ is written:
Kˆ
(
Pi, wi, Pj, wj
)
= Kˆ
(
Pi, Pj
)
wj (4.23)
To adapt the transition kernel to weighted methods, the majorant tech
niques of [122] are used. Provided the coagulation kernel satisfies:
Kˆ
(
Pi, wi, Pj, wj
) ≤ h(1)1 (Pi, wi) h(1)2 (Pj, wj)
+ h(2)1 (Pi, wi) h
(2)
2 (Pj, wj)
+ . . .+ h(Nh)1 (Pi, wi) h
(Nh)
2 (Pj, wj)
(4.24)
where h(j)i is part of an upper bound for Kˆ and Nh is the number of terms
in the factorised expression, the majorant rate can be factorised as:
N(t)∑
i,j=1
Kˆ ≤ λ(1)1 λ(1)2 + λ(2)1 λ(2)2 + . . .+ λ
(Nh)
1 λ
(Nh)
2 (4.25)
where:
λ(k)l =
N(t)∑
i=1
h(k)l (Pi, wi), k, l = 1, 2, . . . , Nh (4.26)
The quantities h(k)l (Pi, wi) are stored in a binary tree which enables rapid
calculation of the sums λ(k)l [54, 122]. They are also used to select the pair
of particles for a coagulation event: more information on such processes is
given in [54]. To apply this formulation to the transition kernel, it needs to
be considered separately as its two components.
47
Freemolecular kernel
The freemolecular kernel is given in Equation (2.23). In this chapter, it will
be written as:
Kfm(Pi, Pj) = α
(
1
mi
+
1
mj
) 1
2 (
di + dj
)2 (4.27)
where the constant α is given by:
α = 2.2
√
pikB T
2
(4.28)
This is decomposed into its majorant rate form Kˆfm(Pi, Pj) to allow for
efficient jump rate calculation:
Kˆfm(Pi, Pj) = α kmaj
(
1√
mi
+
1√mj
)(
d2i + d
2
j
)
(4.29)
where kmaj is the majorant rate scaling factor (kmaj = 2.0), necessary to
satisfy the inequality Kfm ≤ Kˆfm [53]. To calculate the total majorant rate,∑
i 6=j Kˆfm(Pi, Pj)wj must be calculated and split into the individual jump
terms. Using Equation (4.29), we obtain:
∑
i6=j
Kˆfm(Pi, Pj)wj =
α kmaj
{
(N(t) – 1)
∑
d2i m
– 12
i wi
+
[∑
d2i
∑
m
– 12
i wi –
∑
d2i m
– 12
i wi
]
+
[∑
d2i wi
∑
m
– 12
i –
∑
d2i m
– 12
i wi
]
+
[∑
d2i m
– 12
i
∑
wi –
∑
d2i m
– 12
i wi
]}
(4.30)
The indices are neglected here for clarity, however all sums act from i =
1 to N(t). For the freemolecular kernel, there are four λ(1)l λ
(2)
l terms, giv
ing Nh = 4. Note that a
∑
d2i m
– 12
i wi term is subtracted from the other terms
to remove selfcoagulations. The majorant jump rate of freemolecular co
agulation is therefore given by:
Rˆfm =
N(t)∑
i 6=j
Kˆfm(Pi, Pj)wj
Vsmp
(4.31)
48
Slipflow kernel
The slip flow kernel from Equation (2.24) is written in this chapter as:
Ksf(Pi, Pj) = β1
(
1
di
+
1
dj
+ β2
[
1
d2i
+
1
d2j
])
(di + dj) (4.32)
where β1 and β2 are constants, expressed as:
β1 =
2 kB T
3µ
, β2 = 1.257(2σ) (4.33)
The slipflow kernel does not need a majorant due to its simple form [121].
Thus, the sum can be directly expanded:∑
i 6=j
Ksf(Pi, Pj)wj =
β1
{
2(N(t) – 1)
∑
wi
+
[∑
di
∑
d–1i wi –
∑
wi
]
+
[∑
diwi
∑
d–1i –
∑
wi
]
+β2
(
(N(t) – 1)
∑
d–1wi
+
[∑
di
∑
d–2i wi –
∑
d–1wi
]
+
[∑
diwi
∑
d–2i –
∑
d–1wi
]
+
[∑
d–1i
∑
wi –
∑
d–1wi
])}
(4.34)
The jump rate of slipflow coagulation is thus given by:
Rsf =
N(t)∑
i 6=j
Ksf(Pi, Pj)wj
Vsmp
(4.35)
Particle selection algorithm
Equations (4.30) and (4.34) comprise the majorant rate expression for the
‘weighted majorant transition kernel’:
Kˆ(Pi, Pj)wj =
∑
i6=j Kˆfm(Pi, Pj)wj if Rˆfm < Rsf∑
i6=j Ksf(Pi, Pj)wj if Rˆfm ≥ Rsf
(4.36)
49
In a similar vein to [53], these expressions can be thought of as sums of
individual rate terms. For example, Equation (4.30) can be rewritten as:
Kˆfm(Pi, Pj)wj = Kˆ
fm,1 + Kˆfm,2 + Kˆfm,3 + Kˆfm,4 (4.37)
where, as an illustration, the first freemolecular jump rate term is given by:
Kˆfm,1 = α kmaj(N(t) – 1)
N(t)∑
i=1
d2i m
– 12
i wi (4.38)
Comparing Equation (4.30) to Equation (4.25), we see that each of these
individual rate terms Kˆ(fm,i) is composed of two parts λ(fm,1)1 λ
(fm,1)
2 . For
Kˆfm,1, these two terms are given by:
λ(fm,1)1 = (N(t) – 1)
λ(fm,1)2 =
N(t)∑
i=1
d2i m
– 12
i wi (4.39)
which yields the following selection properties h(fm,1)1 , h
(fm,1)
2 :
h(fm,1)1 = 1
h(fm,1)2 = d
2
i m
– 12
i wi (4.40)
A selection property of 1 indicates uniform selection probability; that is,
each computational particle has an equal chance for being selected, regard
less of its physical properties. The coagulating partners i and j are selected
according to the probabilities:
h(fm,1)1 (Pi, wi)
λ(fm,1)1
and
h(fm,2)2 (Pj, wj)
λ(fm,1)2
(4.41)
The particle selection properties and partial sums for the remainder of
the freemolecular and slipflow kernels are summarised in Table 4.1.
Coagulation jump process
If a coagulation event is selected to occur, two particles must be chosen from
the ensemble and physically joined. This process is represented by Equation
4.20. The weight transfer functions γ(Pi, wi, Pj, wj) define how the statistical
50
Table 4.1: Particle selection properties for the weighted transition kernel
Term Equation Particle 1 Particle 2
FM1 (N(t) – 1)
∑
d2i m
– 12
i wi Uniform d
2
i m
– 12
i wi
FM2
∑
d2i
∑
m
– 12
i wi –
∑
d2i m
– 12
i wi d
2
i m
– 12
i wi
FM3
∑
d2i wi
∑
m
– 12
i –
∑
d2i m
– 12
i wi m
– 12
i d
2
i wi
FM4
∑
d2i m
– 12
i
∑
wi –
∑
d2i m
– 12
i wi d
2
i m
– 12
i wi
SF1 2(N(t) – 1)
∑
wi Uniform wi
SF2
∑
di
∑
d–1i wi –
∑
wi di d–1i wi
SF3
∑
diwi
∑
d–1i –
∑
wi d–1i diwi
SF4 (N(t) – 1)
∑
d–1i wi Uniform d
–1
i wi
SF5
∑
di
∑
d–2i wi –
∑
d–1i wi di d
–2
i wi
SF6
∑
diwi
∑
d–2i –
∑
d–1i wi d
–2
i diwi
SF7
∑
d–1i
∑
wi –
∑
d–1i wi d
–1
i wi
Table 4.2: Definition of the coagulation weight transfer functions.
Name Formula
SWA1 wi
wj
wi + wj
SWA3 wi
mi
mi + mj
51
weights are manipulated through such a process and are summarised in
Table 4.2 [122].
The two best performing weight transfer functions (SWA1 and SWA3)
from the family derived in [122] were selected for the numerical tests in
this chapter.
In a stochastic coagulation process, two particles are first selected to co
agulate. For direct simulation, the mass of the second particle is added to
the first, and the ‘old’ second particle is deleted. As discussed in the intro
duction, the ensemble will eventually deplete if not somehow replenished.
This is done by doubling the ensemble when the number of computational
particles N(t) falls below half of the maximum capacity Nmax. When using
SWAs, the weight of the first particle is adjusted (using the functions in Table
4.2), then the mass of the second particle is added to the first, and the sec
ond particle is left unchanged. As no particles are deleted, there is no need
for a doubling algorithm. The ensemble properties cached in the binary tree
structure are updated at the end of the coagulation event [54, 122]. The
complete process in which a coagulation event is handled is provided in
Figure 4.1.
4.3 Numerical studies
While the stochastic weighted algorithms presented in the previous section
are general and may be applied to a stochastic population balance method
with any particle model, the multivariate silica particle model is used here
as a testcase. A simple case was defined in [143] in which numerical con
vergence studies were conducted. This case simulated a zerodimensional
batch reactor with 250 ppm of precursor (tetraethoxysilane, TEOS) in nitro
gen bath gas. The temperature was held constant at 900 ◦C and a residence
time of 0.8 s was chosen. The numerical parameters, model parameters and
process settings used for this test case are summarised in Table 4.3.
As in [143], the convergence of several key metrics was studied as a
function of the number of computational particles Nmax and runs L. These
are given in Table 4.4, along with their physical interpretations. For a sim
ulation with a given (Nmax, L), the runaverage (Equation (2.28)) and con
fidence interval (Equation (2.29)) of each metric was determined. These
results are given in the following section.
52
Start Select Pq,Pr
Calculate Kˆ tr
Perform LPDA updates
Calculate K tr
is Fictitious? wq ← γ(Pq,wq,Pr ,wr )wq
Pq ← Pq + Pr
Pq ← Pq + Pr
Delete Pr
N(t) < 12Nmax
Duplicate
ensemble
Update ensemble
statistics
Finish
TRUE
(DSA)
TRUE
(SWA)
FALSE
TRUE
FALSE
Figure 4.1: Comparison of the coagulation process implemented in direct
simulation algorithm and weighted particle methods.
53
Table 4.3: Numerical and model parameters used for the silica test case.
Description Symbol Value Ref.
Numerical parameters
Number of splits Nsplits 100 [143]
Timestep ∆t 0.0015 s [143]
Max. number of stochastic particles Nmax 64–131072 [143]
Number of runs L 1–2048 [143]
Model parameters
Maximum zeroth moment M0,max 1.6× 1018 #/m3 [143]
Sintering preexponential AS 1.1× 1016 s/m [143]
Sintering characteristic energy ES 1.2× 105 K [143]
Sintering minimum diameter dp,min 4.4 nm [143]
Surface smoothing factor σSR 1.0 [142]
Gasphase parameters   [142]
Freemolecular enhancement factor kfm 2.2 [66]
Process settings
Initial temperature T 900 ◦C [143]
Residence time τ 0.8 s [143]
Initial TEOS fraction yTEOS 250 ppm [143]
Initial total pressure Pi 1.0 atm [143]
54
Table 4.4: Summary of key process metrics (from Shekar et al. [144]) in
vestigated in this chapter.
X(t) Description Formula
M0(t) Particle number concentration
1
Vsmp
N(t)∑
i=1
wi
FV(t) Particle volume fraction
1
Vsmp
N(t)∑
i=1
wi Vi
µa(n)(t) Mean number of primaries per particle
1∑
wi
N(t)∑
i=1
wi ni
µa(dcol)(t) Mean collision diameter
1∑
wi
N(t)∑
i=1
wi dcol,i
µa(dpri)(t) Mean average primary diameter
1∑
wi
N(t)∑
i=1
widprii
µa(s)(t) Mean average sintering level
1∑
wi
N(t)∑
i=1
wi si
55
109
108
107
106
105
104
103
102
101
100
0 2 4 6 8 10 12 14 16
ke
rn
el
d
en
sit
y,
1
/n
m
d pri
, nm
109
108
107
106
105
104
103
102
101
100
0 100 200 300 400
ke
rn
el
d
en
sit
y,
1
/n
m
dcol, nm
DSA
SWA1
SWA3
Figure 4.2: PSDs predicted by the highprecision solutions at t = 0.8 s.
4.4 Results
4.4.1 Highprecision solution
Before investigating any useful convergence properties of the weighted meth
ods, one must ascertain whether the solutions when using SWAs converge
in the limit to those solutions when using direct simulation. Here, the high
precision solutions (Nmax = 131072, L = 10) are compared.
The particle size distributions (PSDs) predicted by the algorithms were
first examined. As the silica model tracks primary particles, a PSD of the
primary particles may be generated in addition to PSDs for the aggregates.
The highprecision primary and collision PSDs are given for direct simula
tion and the weighted methods in Figure 4.2. These were generated using a
Gaussian kernel density estimation procedure, and the normal distribution
approximation to estimate the bandwidth [146].
It is evident that the SWA1 and SWA3 algorithms perfectly reproduce
the PSD predicted by direct simulation. Further, SWA1/3 better resolve the
large (>300 nm) and rare particles.
The metrics in Table 4.4 also provide a useful tool to examine agreement
between algorithms. The temporal evolution of such metrics for the direct
simulation and weighted algorithms is presented in Figure 4.3.
Again, the SWA1 and SWA3 algorithms produce trajectories with excel
56
1014
1015
1016
1017
1018
1019
M
0
(t)
, 1
/m
3
DSA
SWA1
SWA3
109
108
F V
(t)
, 
0
20
40
60
80
100
µ a
(d c
o
l) (
t),
nm
0
1
2
3
4
5
6
7
8
9
µ a
(d p
ri
) (
t),
nm
0
10
20
30
40
50
60
70
80
90
0.0 0.2 0.4 0.6 0.8
µ a
(n)
(t)
, 
time, s
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
0.0 0.2 0.4 0.6 0.8
µ a
(s )
(t)
, 
time, s
Figure 4.3: Temporal evolution of Table 4.1’s key metrics for the high
precision solutions using DSA and SWAs. The two lines show
the upper and lowerbound of the 99.9% confidence interval.
57
lent agreement with those of DSA. The confidence intervals for these cases
are typically too narrow to be observed here. It is clear that the SWAs’ high
precision solutions are equivalent to those predicted through direct simula
tion, indicating that all algorithms approach the same solution in the limit
as Nmax →∞.
4.4.2 Convergence of algorithms
When using stochastic particle methods, one must choose L such that the sta
tistical uncertainty in the results is acceptably small. One must also choose
Nmax such that the systematic error appears close to zero, however in prac
tice it may be chosen such that the differences between successive values
of Nmax are not statistically significant. The multidimensional silica particle
model permits tracking of any number of functionals describing the particle
size and morphology. Here, the convergence behaviour of these metrics is
investigated.
The convergence studies reported here were conducted by varying Nmax
and L while holding their product Nmax × L constant at a value of 217. The
relationship between the final average value of these functionals and Nmax
is given in Figure 4.4. The uncertainty in each point is depicted using the
confidence interval I0.999(tf).
Particle number concentration
The final value of the SWA solutions for particle number concentra
tion M0 (or zeroth moment) has lower systematic error than direct
simulation.
Volume fraction
The volume fraction represents the ratio of total particle volume to
gas volume. The SWAs achieve a solution within the confidence in
terval of the highprecision solution for Nmax ≥ 128, indicating that
they accurately predict the volume fraction with very few stochastic
particles.
Mean number of primaries per particle
The mean number of primaries per particles µa(n) is a measure of both
the size of the aggregate and number of coagulation events. Figure
4.4 shows that much smaller aggregates are obtained when Nmax is
too low.
58
2.1
2.2
2.3
2.4
M
0
(t f)
, 1
01
4
1/
m
3
DSA
SWA1
SWA3
6.8
7.0
7.2
7.4
7.6
F V
(t f)
, 1
09

76
78
80
82
84
86
88
µ a
(d c
o
l) (
t f),
nm
8
8.5
9
9.5
10
µ a
(d p
ri
) (
t f),
nm
58
60
62
64
66
68
70
72
74
76
78
101 102 103 104 105 106
µ a
(n)
(t f
), 
Nmax, 
0.05
0.06
0.07
0.08
101 102 103 104 105 106
µ a
(s )
(t f
), 
Nmax, 
Figure 4.4: Runaveraged values of the functionals studied (Table 4.4)
and their corresponding confidence intervals at t = 0.8s with
Nmax × L = 217.
59
104
103
102
101
ε s
ta
t(M
0)
(t f)
, 
DSA
SWA1
SWA3
104
103
102
101
ε s
ta
t(F
V)
(t f)
, 
104
103
102
101
ε s
ta
t(µ
a
(d c
o
l))
(t f)
, 
104
103
102
101
ε s
ta
t(µ
a
(d p
ri
))
(t f)
, 
104
103
102
101
101 102 103 104 105
ε s
ta
t(µ
a
(n)
) (
t f),

Nmax, 
104
103
102
101
101 102 103 104 105
ε s
ta
t(µ
a
(s )
) (
t f),

Nmax, 
Figure 4.5: Relative statistical error of final values for Table 4.3’s metrics
with constant Nmax × L = 217.
60
Mean collision diameter
The mean collision diameter µa(dcol) is an important property which
characterises the size of an aggregate particle. SWAs appear to have
better accuracy at lower Nmax, but there is a common limit as Nmax
becomes large.
Mean primary diameter
The mean average primary diameter µa(dpri) estimates the average
size of the primary particles which compose an aggregate. Curiously,
despite SWAs predicting the final value of the collision diameter bet
ter than DSA, they have a larger systematic error than DSA for fewer
computational particles. This may be associated with SWAs preference
towards larger particles (with larger primaries).
Mean sintering level
The mean average sintering level µa(s) is a measure of the extent to
which the ensemble’s particles are sintered. SWAs and DSA show sim
ilar convergence properties for this metric.
In summary, the weighted methods show convergence for fewer stochas
tic particles for M0, FV, µa(dcol), and the sintering level. Figure 4.4 demon
strates that the SWAs consistently have a smaller confidence interval than
an equivalent DSA simulation with a given Nmax, L. This is consistent with
observations made in [122].
In any case, it is evident that Table 4.4’s functionals satisfactorily con
verge to a stable value. For this system, a choice of Nmax ≥ 16384 would be
recommended when using SWAs.
4.4.3 Computational efficiency
The computational time required to simulate the silica test system for a
given Nmax, L, or error is investigated in this section. These calculations
were conducted on a SGI Altix Cluster with nodes composed of two 3.00
GHz quad core Intel Harpertown processors and 8 GB of RAM.
The computational times of the DSA and SWA simulations for increasing
Nmax are compared in Figure 4.6. Typically, a simulation with a given Nmax
and L = 1 using SWAs required six to eight more times the CPU time than
61
100
101
102
103
104
105
106
101 102 103 104 105 106
CP
U
tim
e
pe
r r
un
, s
Nmax, 
DSA
SWA1
SWA3
Figure 4.6: Average CPU time required per run as a function of number of
stochastic particles.
the equivalent DSA simulation. Further, it was observed that the time re
quired per run increases approximately linearly with the maximum number
of stochastic particles.
It is also no coincidence that the additional time required by the SWAs as
compared to the DSA is proportional to the number of additional simulated
coagulation events, which are shown in Figure 4.7. Firstly, as SWAs empha
sise calculation of larger, rarer particles (e.g. [122] or Figure 4.2), more
time on average is required to determine the extent of sintering (Equations
(4.10) and (4.13) ) of the ensemble’s particles at the end of each timestep.
Secondly, as SWAs always have a full ensemble (N(t) = Nmax), there are
always more particles for which the process jump rates and sintering char
acteristics must be evaluated.
The additional numerical coagulation events that occur during SWA com
putations are a consequence of the redistributed statistical accuracy, that
also reduces the variance of the studied metrics. This is depicted in Fig
ure 4.5, where the confidence interval widths (or relative statistical error
stat) for the metrics’ final values are compared across algorithms. Gener
ally, a constant confidence interval is obtained, with a few notable excep
62
0
1000
2000
3000
4000
5000
6000
7000
8000
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
m
a
xi
m
um
n
um
be
r o
f c
oa
gu
la
tio
n
ev
en
ts
, 
time, s
DSA
SWA1
SWA3
Figure 4.7: Comparison of the maximum number of coagulation events un
dergone by particles for direct simulation and the SWAs for
Nmax = 131072.
63
tions. Those metrics which show steadily increasing or decreasing statistical
errors correspond to the metrics whose final values (at Nmax small) is sig
nificantly different (∼ 10%) from the highprecision solution. That is, the
presence of systematic error in the confidence interval when a simulation is
unconverged affects the estimation of statistical error.
The difference in size of DSA confidence intervals with respect to the
SWA intervals is dependent on the metric analysed. For the zeroth moment
they are the same, however for primary diameter they differ by a factor of
about 4. It is evident that the weighted methods offer superior variance
reduction for all of the unique properties tracked in the multidimensional
particle model.
4.5 Conclusions
This chapter has used majorant techniques to adapt the transition regime
coagulation kernel to weighted particle methods. This is an important de
velopment if SWAs are to be used to solve realworld coagulation problems.
This kernel is used in conjunction with a stateoftheart multidimensional
silica particle model to investigate the numerical behaviour of the popula
tion balance solution methodology.
Within the particle model, particles are described by the primary parti
cles and their connectivity, which are in turn described by the number of Si,
O and OH units of which they are composed. This permits tracking of many
detailed features of the particles and particle ensemble through time.
Two SWAs were tested for convergence to direct simulation in the limit.
Excellent agreement of all algorithms’ highprecision solutions was obtained,
indicating that the weighted algorithms are expected to converge to the di
rect simulation solution for sufficiently large number of stochastic particles.
The convergence of the functionals describing particle properties was
investigated with respect to the maximum number of computational parti
cles in the ensemble. SWAs generally achieved the correct value of these
functionals for fewer stochastic particles than DSA, except in the calculation
of number of primary particles and the primary diameter. For a calcula
tion with a given number of stochastic particles and runs, the confidence
intervals of the functionals were up to four times smaller, indicating that
weighted algorithms studied here achieve better variance reduction than di
64
rect simulation. However, this comes at the cost of increased computational
expense, with the algorithms requiring approximately seven times more CPU
time per run than direct simulation.
In this chapter, it was demonstrated that stochastic weighted algorithms
can be successfully used to solve population balances with realworld co
agulation kernels and complex particle models. This is particularly impor
tant for spatially inhomogeneous systems, where weighted methods offer
practical and numerical advantages over direct simulation. A new majorant
formulation of the transition kernel, as adapted to stochastic weighted al
gorithms was presented. A detailed numerical analysis of the convergence
properties and efficiency of these methods was conducted, showing that
weighted methods satisfactorily converge to a stable numerical solution.
65
Chapter 5
Solution of reactor networks
An algorithm to solve networks of reactors with population bal
ance models is presented. The convergence of the stochastic particle
method is evaluated as a function of network size, feedback and
numerical parameters.
Material from this chapter is published in [94].
Population balance modelling has traditionally been applied to mod
elling particle formation in a batch or a plugflow reactor [8, 92, 99]. In the
former case, the equations governing particle growth can be simplified to ex
clude terms for particle transport to and from the reactor. For the latter case,
the reactor model typically represents an axial streamline through the plug
flow reactor (PFR). However, there are many systems found in engineering
where advective or diffusive particle transport is important [82, 119]. This
typically requires use of computational fluid dynamics (CFD) or a reactor
network approach, coupled with a population balance model to accurately
capture particle dynamics.
Solutions for the population balance equation in these coupled systems
may be grouped into several classes, the most common of which are monodis
perse [18, 68], moment [3, 4, 48, 89] and sectional [130] methods. All
of these methods yield a set of differential equations which may be solved
within the framework of the differential equations describing fluid transport
and/or chemical reactions [82].
Many commercial packages now include particle population balance mod
els using the aforementioned methods, for example FLUENT [6], CHEMKIN
PRO [131] and gPROMS [10]. The major drawback of using these methods
66
to solve the population balance equations is that they either do not resolve
the particle size distribution (e.g. moments) or they utilise a simplified par
ticle model such as a spherical or surfacevolume type description. This can
cause errors when particle aggregates are formed through coagulation and
sintering [91].
Implementation of stochastic population balance solvers for flow sys
tems requires discretising the reactor geometry into a mesh composed of
reactor cells [82, 119]. These connected cells can be considered a type of
reactor network. Kruis et al. [82] developed a ‘cellbased weighted random
walk’ method to account for the movement of particles between cells in one
and twodimensional reactor geometries. While reasonably good agreement
with a moments code was reported, it was observed that the method was
sensitive to stochastic noise, a result of using too few computational parti
cles.
Patterson and Wagner [119] presented a stochastic method for coagulation
advection problems. The performance in accounting for coagulation of the
conventional DSA against SWAs was evaluated. It was reported that DSA is
particularly susceptible to stochastic noise, which eventually manifests itself
as a systematic variance from the true solution.
A study by Zhao and Zheng [169] gave a threedimensional adaptation
of a stochastic population balance solver (using an SWA [170]) coupled with
a CFD code. An algorithm was presented to decouple the fluid flow, particle
flow and particle dynamics by choosing an appropriate timestep size. In a
similar vein to [82], limits associated with computational cost restricted the
maximum number of computational particles, preventing precise agreement
with direct numerical simulation (DNS).
‘Multizonal’ [11, 12] or ‘compartmental’ models [5, 67, 84, 167] are
models in which the fine mesh of CFD is simplified to a coarse network of
reactors with a specified flow between each node. Such techniques typically
require orders of magnitude fewer reactor cells, thus reducing computa
tional cost, however this comes at the price of decreased spatial resolution
of particle dynamics.
Only recently have these models begun to be investigated using stochas
tic methods. Irizarry [67] recently presented an approach to solve a popu
lation balance compartmental model using a stochastic method. A new al
gorithm to move particles between cells (‘particle bundle flow’ method) was
67
reported, in which particles were restricted from jumping more than one
compartment by choosing a suitably small timestep. The numerical prop
erties of some test systems were evaluated as a function of the timestep,
however the effect of other numerical parameters was not addressed.
There are additional open questions in solving reactor networks with
stochastic population balance solvers. In order to maximise efficiency of
these algorithms, it is necessary to understand the convergence of these
algorithms with respect to the number of computational particles and inde
pendent runs [143]. Further, it is important to quantify and reduce stochas
tic noise in the system, as this was identified as a potential issue in [82].
Finally, the coupling of a such a system to a reacting gasphase has not yet
been investigated.
The purpose of this chapter is to extend the work of Irizarry [67], Kruis
et al. [82] and Zhao and Zheng [169] by developing an algorithm for solving
a fullycoupled gasphase ODE/particle population balance network. The
convergence properties of the stochastic population balance will first be in
vestigated, and methods to reduce statistical error will be considered. The
flexibility of stochastic methods in solving the population balance equation
will be demonstrated by use of a detailed multivariate particle model to
simulate silicon nanoparticle growth in a plugflow reactor.
The structure of this chapter is as follows. In Section 5.1, the reactor
models are presented. The stochastic solution methodology and model used
for generating reference solutions are discussed in Sections 5.1.1 and 5.1.2,
respectively. The numerical test cases and their results are discussed in Sec
tion 5.2. Finally, Section 5.3 illustrates one of the many ways in which
reactor networks can be applied to model a real reactor system.
5.1 Model
This work will study the numerical and physical behaviour of a constant
pressure network of perfectlystirred reactors (PSRs; or continuousstirred
tank reactors, CSTRs). It is assumed that each of these reactors has a char
acteristic residence time τPSR = Q/V associated with it, where Q is the
volumetric flowrate and V is the reactor volume. The inflow and outflow
terms modify Equation (2.13) to:
68
d
dt
n(Pq) = I(c, Pq) +K(n, Pq) +R(c, n, Pq) + S(n, Pq)
+
1
τPSR
Nin∑
j=1
f [j]n[j](Pq) – n(Pq)
– γ(c, n, T)n(Pq) (5.1)
where the reactor is assumed to have Nin inflow streams, each with fraction
f [j] of the total volumetric flow. The operators I(c, Pq), K(n, Pq),R(c, n, Pq)
and S(n, Pq) take the same form as in Equations (2.15), (2.17), (2.19) and
(2.21) respectively. For NGP gasphase species, the rate of change in concen
tration of species k is given by:
dck
dt
= ω˙k(c, T) + g˙k(c, n, T) +
1
τPSR
Nin∑
j=1
f [j]c[j]k – ck
– γ(c, n, T)ck (5.2)
where ω˙k(c, T) and g˙k(c, n, T) are the molar production rates of species k
at temperature T due to gasphase and particle reactions, respectively. The
gasphase expansion factor for a PSR is given by:
γ(c, n, T) =
1
ρ
NGP∑
k=1
[
ω˙k + g˙k
]
+
1
τPSR
Nin∑
j=1
f [j]ρ[j] – ρ
+ 1
T
dT
dt
(5.3)
where ρ is the molar density of the mixture. Ignoring energy changes due to
the particle population balance, the adiabatic energy balance of a constant
pressure PSR is:
ρCP
dT
dt
=
NGP∑
k=1
–Hˆkω˙k + 1τPSR
Nin∑
j=1
f [j]c[j]k
[
Hˆ[j]k – Hˆk
] (5.4)
where CP is the bulk heat capacity and Hˆk is the specific molar enthalpy of
species k.
5.1.1 Stochastic method
Typespace
The spherical typespace presented in Chapter 3 will be used for the numer
ical and case studies of this chapter. Here, the typespace will be written
as:
Pq = Pq (p(η1)) (5.5)
69
where computational particle q contains one primary p with an arbitrary
component η1. For this very simple particle model, the typespace can be
simplified to:
Pq = Pq (p(η1)) = Pq
(
v(q)
)
(5.6)
Particle processes
In the stochastic method, the process which create and change particles
are typically implemented as jump processes [91]. A very simple particle
mechanism will be used in this chapter to test the reactor network imple
mentation. It is briefly discussed here.
Inception
Particles are created in the population balance by inception:
χ1 + χ1 → Pq(v1) (5.7)
where χ1 is a nucleating gasphase species, and v1 is the size of the
smallest possible particle.
Coagulation
Particles coagulate and stick together according to the following reac
tion:
Pq
(
v(q)
)
+ Pr
(
v(r)
)
→ Ps
(
v(q) + v(r)
)
(5.8)
Both the constant kernel (Kconst = 1) and transition kernel (Ktr) are
used in this chapter.
Surface growth
A surface reaction event may occur in the following manner:
Pq
(
v(q)
)
+ χ1 → Pr
(
v(r) + δvA
)
+ χ2 (5.9)
where v(q) = v(r), χ1 is a reacting gasphase species, χ2 is a product
species, and δχA is the volume element added due to the reaction.
Algorithms for coagulation
The key challenge in stochastic particle methods is how to solve the Smolu
chowski coagulation equation [93]. In this chapter, DSA and a stochastic
weighted algorithm (SWA1 of Table 4.2) will again be compared.
70
Flow processes
Previous studies account for particle transport by moving computational
particles between cells [67, 169]. Alternative formations for accomplishing
this task are presented in the following section. An inflow process uniformly
selects a particle of index q from inflow stream j, identifying the particle to
be copied into the ensemble. This is represented by the reaction:
(Pq, wq)[j] → (Pr, wq), (Pr+1, wq), . . . , (Pr+nc , wq) (5.10)
where Pr is a copy of Pq, and nc copies of the particle are made. Note that
the particle weights remain unchanged in the inflow process. The number
of copies to add is determined by the scaling factor:
nc =
{
bFcc+ 1 for u < Fc – bFcc
bFcc otherwise
(5.11)
where u is a pseudorandomly chosen real number such that u ∈ [0, 1).
Adding multiple copies of particles in this manner has the advantage of
accelerating the computational particle count N(t) to the maximum value
Nmax, thus preventing depletion of the ensemble (when using the DSA)
should the outflow rate be similar. The quantity Fc represents the ratio of
the current to inflow ensemble scaling factors:
Fc =
Vsmp
V[j]smp
(5.12)
The rate of an inflow or outflow process is proportional to the number of
stochastic particles in the ensemble [13, 25]. For a constantpressure reactor
with constant volumetric inflow, the rate of inflow from stream j is given as:
R[j]inflow = f
[j] N
[j](t)
τPSR
(5.13)
where N[j](t) is the number of stochastic particles in the inflow stream. Two
possible implementations of an outflow are now presented. Firstly, upon
selection of an outflow jump process, particles can simply be deleted:
(Pq, wq)→ deleted (5.14)
The alternative requires consideration of flow processes as continuous rather
than as a stochastic jump process. Supposing particles are lost due to an out
flow process in time ∆t, the ensemble scaling factor Vsmp can be scaled by
71
factor Fs to represent a decrease in the real particle number concentration,
as M0 = N(t)/Vsmp [93]. The rescaling factor Fs is given by:
Fs =
1
1 – ∆t/τPSR
(5.15)
The new sample volume is obtained using Vnewsmp = FsVsmp. This form of
particle outflow is termed a ‘rescale’ outflow. It should also be observed that
the outflow process of an upstream reactor does not necessarily provide the
inflow term of a downstream reactor. The rate of particle outflow from a
reactor is given by:
Routflow =
N(t)
τPSR
. (5.16)
The inflow and outflow (in the case of particle deletion) may be included
in the stochastic particle method as conventional jump processes (e.g. in
Algorithm 2 of [143]). It should also be noted that copying and deletion of
computational particles makes tracking of a particle’s trajectory through the
network challenging.
Coupling to the gasphase
The stochastic particle method is fullycoupled to a gasphase ODE solver.
The technique of operator splitting is employed to couple the gas and
particlephases. The algorithm by which this is done (Strang splitting) has
been described in detail by Celnik et al. [25] and Shekar et al. [143].
5.1.2 Discrete Model
In order to assess the numerical properties of the stochastic method, a re
liable reference model is needed. In this work, the Discrete Model will be
used, which writes the population balance equation as a series of ordinary
differential equations (ODEs) [48]. These may be solved using conventional
differential equation packages. The Discrete Model is similar to a sectional
method, however every particle size is represented by its own differential
equation. Arbitrarily high precision can be obtained by choosing a suffi
ciently large maximum particle size, corresponding to the number of ODEs,
NODE.
For the test particle mechanism, constant inception (Kconst = 1) and
constant surface growth rates (RSG = 1) will be assumed. The rate of coag
ulation is determined using the Smoluchowski coagulation equation, and is
72
given in [48] as:
K
(
n, Pi(v
(i))
)
=
1
2
i–1∑
j=1
Kj,i–jnjni–j – ni
NODE∑
j=1
Ki,jnj (5.17)
where Ki,j is the coagulation kernel evaluated for particle pair (Pi, Pj). Sup
posing that particles of size i ∈ [1, NODE] are spherical with volume iv1, and
that coagulation occurs with the constant kernel, the population balance
equation is written as:
dni
dt
=
Kconst – Kconstni
NODE∑
j=1
nj – RSGni
+
1
τPSR
Nin∑
j=1
f [j]n[j]i – ni
i = 1
1
2
Kconst
i–1∑
j=1
njni–j – K
constni
NODE∑
j=1
nj + RSG
(
ni–1 – ni
)
+
1
τPSR
Nin∑
j=1
f [j]n[j]i – ni
i > 1
5.1.3 Solution of reactor networks
The formulation of the above flow processes permits linking of reactors in
a generic manner. To solve the transient response of the network, a simple
sequential procedure is adopted. This is termed ‘sequentialmodular simu
lation’, and similar approaches have been used in other applications [169].
The algorithm employed to solve the fullycoupled gasphase ODE system
with the stochastic population balance solver for Nreac reactors is given in
Figure 5.1.
The algorithm in Figure 5.1 illustrates the solution methodology for
Strang splitting, although in theory any form of coupling, such as predictor
corrector [26], could be used. The accuracy of this algorithm is clearly
dependent on choice of a sufficiently small splitting timestep ∆ts, which has
been investigated elsewhere [67, 169].
73
Start
j ← 0
i ← 0
i < Nreactj < tf
j ← j + 1
i ← i + 1
Finish
Loop over reactors
Loop over timesteps
Solve gasphase over
[
tj , tj + ∆ts2
]
Solve particlephase over
[
tj , tj +∆ts
]
Solve gasphase over
[
tj + ∆ts2 , tj +∆ts
]
TRUEFALSE
TRUE
FALSE
Figure 5.1: Algorithm used to solve the reactor network with Strang split
ting.
5.2 Numerical studies
As in Chapter 4, the statistical quantities from Section 2.4 will be used to
calculate the runaveraged mean and confidence intervals of a metric X(t).
Equations (2.31), (2.32) and (2.33) were additionally used to estimate the
total average normalised error (), normalised statistical error (stat(t)) and
average normalised statistical error (stat); respectively. The functionals
studied in this chapter are given in Table 5.1.
5.2.1 Comparison to Discrete Model
To test the validity of the stochastic particle algorithm in solving the popu
lation balance equations, a simple test system was devised: three reactors
connected in series, with possible countercurrent recycle streams opposing
the direction of flow. Material is recycled from reactors R2 and R3 with frac
tion fR, and the total volumetric flowrate Q is constant between all reactors.
It was assumed that the inflow to the network contained 1 #/m3 particles
of size i = 1, with each reactor having the same residence time τPSR. A
schematic of this system is provided in Figure 5.2(a).
This system was solved using a simple particle mechanism including
74
Table 5.1: Summary of process metrics studied in the present work.
X(t) Description Formula
M0(t) Particle number concentration
1
Vsmp
N(t)∑
i=1
wi
µa(v)(t) Average particle volume
1∑
wi
N(t)∑
i=1
v(i)
M2(t) Second mass moment
1
Vsmp
N(t)∑
i=1
wi
(
ρv(i)
)2
M3(t) Third mass moment
1
Vsmp
N(t)∑
i=1
wi
(
ρv(i)
)3
(1− fR)Q
Inflow
Outflow
fRQ
fRQ
R1 R2 R3 (1− fR)QQ Q
(a) Linear network with recycles
Q0
Inflow
Outflow
QNreacR1 R2 . . . RNreac
Q1
fQ0 fQ1 fQNreac−1
(b) Linear network with multiple inflows
Figure 5.2: Schematic of the test systems used in the present study.
75
terms for coagulation, surface growth and inception, all with rates given
by a constant. While these do not represent physical rate equations, they
are helpful for evaluating the performance of the stochastic particle method
against other solution methodologies. This was first undertaken by applying
the Discrete Model [48] and the stochastic particle method (using the SWA
[122]) to solve the test system. The results of this analysis are given in Fig
ure 5.3, for a linear (fR = 0) and countercurrent (with fR = 0.5) network;
and Nmax = 16384, L = 8.
The transient evolution of M0 illustrates the issue of statistical error
present in stochastic solutions to the population balance equation: the ran
dom nature of the processes in the stochastic method induce fluctuations in
functionals describing the system (e.g. the moments). The noise can be elim
inated by averaging the solution over multiple independent runs (i.e. with a
different random seed). The influence of this statistical error and systematic
variation will be addressed in the following sections.
5.2.2 Linear networks
It is useful to investigate how the error in a stochastic solution varies with
the number of computational particles (Nmax) and independent runs (L).
Further, it is likely that the different coagulation algorithms and two out
flow process types could contribute different amounts of systematic or sta
tistical error to the stochastic solution. Using the test case described in the
previous section, the stochastic solutions using the DSA and the SWA were
generated. It was assumed that reactors R1 and R3 used outflow rescaling
as their outflow process, while the outflow process type of R2 was varied.
The normalised average error in R2 of the stochastic solution with respect
to the Discrete Model was evaluated, with a constant Nmax × L = 10× 218,
and is given in Figure 5.4.
It is difficult to observe any meaningful differences between coagulation
algorithms or outflow processes in examining the error in the zeroth mo
ment (Figure 5.4, left panel). The approximately flat error indicates that
M0 has converged for all systems by Nmax = 256, after which the total er
ror is representative of the statistical error. Higher moments (e.g. M2, M3)
show analogous behaviour. In contrast to this, the error in average particle
volume (Figure 5.4, right panel), and similarly other average particle prop
erties, shows a clear difference between the convergence of the DSA and
76
0.75
0.8
0.85
0.9
0.95
1
1.05
M
0(t
), a
rb.
un
its
R1
SWA: linear
ODE: linear
SWA: fR=0.5ODE: fR=0.5
R2 R3
103
102
101
100
M
2(t
), a
rb.
un
its
103
102
101
100
0 2 4 6 8 10
M
3(t
), a
rb.
un
its
t/τPSR, s
0 2 4 6 8 10
t/τPSR, s
0 2 4 6 8 10
t/τPSR, s
Figure 5.3: Comparison of solutions for transient evolution of normalised
moments M0, M2 and M3 from the Discrete (ODE) Model [48]
and stochastic method (using the SWA) for the threereactor
system (Figure 5.2(a)) with fR = 0 and fR = 0.5. Confidence
intervals of the stochastic solution are given where visible.
77
103
102
101 102 103 104 105 106
ε (M
0),

Nmax, 
DSA:Delete
DSA:Rescale
SWA:Delete
SWA:Rescale
1 slope
103
102
101 102 103 104 105 106
ε (µ
a
(V
)),

Nmax, 
Figure 5.4: Convergence of the zeroth moment (M0) and average volume
(µa(V)) in R2 using the DSA, SWA; and outflow rescaling and
particle deletion.
78
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
M0 µa(V) M2 M3
ε s
ta
t ,

metric type
DSA:Delete
DSA:Rescale
SWA:Delete
SWA:Rescale
Figure 5.5: Comparison of the normalised average statistical error (stat) in
selected process metrics in R2 (fR = 0) using the DSA, SWA;
and outflow rescaling and particle deletion (Nmax = 16384,
L = 160).
SWA. Similar behaviour of the SWA was reported in [93].
The previous discussion focused on the reduction of systematic error as
a function of numerical parameters. Using Equation (2.33), the relative
amount of statistical error present in each case can also be compared. These
results are presented in Figure 5.5.
Only in the zeroth moment does DSA appear to have a slight advan
tage in terms of reducing statistical error. For all other cases (particularly
the higher moments), the SWA shows up to a 45% reduction in average
uncertainty, which is consistent with previous studies [93, 122]. The out
flow rescaling process type appears to contribute less statistical error than
its particle deletion counterpart. This is unsurprising, since the removal of
a particle from the ensemble deletes information that could otherwise be
preserved in a rescaling approach.
79
5.2.3 Networks with recycle loops
Chemical engineering process models may often include a ‘recycle stream’
[46], where material or energy is transferred from a downstream to an up
stream unit operation. The stochastic solution of the population balance
equation has only recently been studied with such networks [67]. As such,
this section aims to investigate some of the preliminary characteristics of the
stochastic particle algorithm where feedback loops are present.
The normalised average error in R2 with respect to the Discrete Model
was evaluated for the network in Figure 5.2(a) using the test particle mech
anism with fR = 0.5 (Nmax × L = 10 × 218), and is presented in Figure
5.6.
In contrast to the results of Figure 5.4, it is evident the system is more
sensitive to the number of computational particles. The linear decrease in
the metrics’ average error (slope of 1) with Nmax is consistent with other
studies of convergence of stochastic particle methods [25, 53, 99]. Only
for the zeroth moment does the DSA outperform the SWA for reduction of
error as a function of Nmax. Both methods appear to converge by the same
threshold; Nmax = 4096.
In addition to recycling material (gasphase and particles), a feedback
loop could also recycle statistical error. Using outflow rescaling, the tran
sient evolution of the statistical error in M3 was evaluated (Equation (2.32))
in R3 for different values of the recycle fraction. These results are displayed
in Figure 5.7.
As the recycle fraction increases, so does the ‘noise’ in the solution. Use
of the SWA is particularly advantageous in this system, in some cases show
ing more than an order of magnitude less statistical error than the DSA.
The same phenomenon is observed for all other metrics, with the exception
of M0. In the case of M0, the statistical error is so small that any poten
tial advantage of DSA is outweighed by the error accumulation in the other
metrics.
Although the statistical errors appear to stabilise as the reactor reaches
steadystate, the statistical error in the DSA solutions is concerningly large;
reaching more than 20% for fR = 0.8. In practical use of stochastic pop
ulation balance modelling, it is uncommon to use 160 independent runs
(8 is more common: [90, 143]), as the computational expense can be too
great. Since the statistical error is dependent on the number of runs, it fol
80
104
103
102
101
101 102 103 104 105 106
ε
(M
0),

Nmax, 
DSA:Delete
DSA:Rescale
SWA:Delete
SWA:Rescale
1 slope
104
103
102
101
101 102 103 104 105 106
ε
(µ a
(V
)),

Nmax, 
103
102
101
100
101
101 102 103 104 105 106
ε
(M
2),

Nmax, 
103
102
101
100
101
101 102 103 104 105 106
ε
(M
3),

Nmax, 
Figure 5.6: Convergence of the zeroth moment (M0), average volume
(µa(V)) and second and third mass moments (M2, M3) for the
second reactor (R2) in a countercurrent network with fR = 0.5.
81
103
102
101
100
101 100 101 102
ε s
ta
t(M
3),

t/τPSR, 
DSA
linear
fR=0.2fR=0.4fR=0.6fR=0.8
103
102
101
100
101 100 101 102
ε s
ta
t(M
3),

t/τPSR, 
SWA
Figure 5.7: Normalised statistical error in R3 as a function of dimensionless
reactor time and recycle fraction (Nmax = 16384, L = 160).
lows from these results that a SWA should be used in place of the DSA for
networks with strong feedback.
5.2.4 Reactors with multiple inlets
The previous analyses addressed linear reactor networks and those which
recycle information through the network for a system with three reactors.
It is also important to understand the effect of increasing the length of the
reactor network and the contribution of material from inflow streams to a
reactor. The Damköhler number is typically used to investigate the latter
phenomenon, and is defined here as:
Da =
1/τ incoag
1/τPSR
(5.18)
where τcoag is the characteristic coagulation time, given by 1/(KconstM0)
[91]. By varying the Damköhler number of a stream, the effect of an in
creased coagulation rate upon the test system can be assessed. To do this,
the system given in Figure 5.2(b) with length Nreac = 10 reactors and inflow
82
1
2
3
4
5
6
7
8
1 2 3 4 5 6 7 8 9 10
ε s
ta
t(µ
a
(V
)),

reactor index, 
Da = 1
×103
L=20
L=40
L=80
L=160
1
2
3
4
5
6
7
8
1 2 3 4 5 6 7 8 9 10
ε s
ta
t(µ
a
(V
)),

reactor index, 
Da = 10
×103
Figure 5.8: Normalised average statistical error in the mean particle vol
ume as function of reactor index (Nmax = 16384, SWA coagu
lation).
fraction f = 0.5 was used. The statistical error in each reactor as a function
of the inflow Damköhler number was determined and is given in Figure 5.8.
The error measurements increase through the length of the PSR chain,
before levelingoff after five reactors. Where more coagulation occurs (higher
Da, right panel), the error is increased. This analysis illustrates how error
can potentially propagate through a reactor chain, although in this particu
lar example the error is sufficiently low to be of no concern. It is suggested
that for systems where particle processes (e.g. coagulation, surface growth)
occur over much shorter characteristic times than the transport time, more
independent runs or computational particles could be used to offset the ac
cumulation of statistical error.
83
5.3 Applications
5.3.1 Approximation of a PFR
Many aerosol synthesis processes use a tubular hotwall reactor or similar
method to manufacture nanoparticles [75, 99, 141, 163]. These processes
are typically modelled using a plugflow reactor (PFR) model, where the
length coordinate of the reactor is translated into a time coordinate to per
mit solution of the population balance equations. However, this formulation
does not permit investigation of the transient behaviour of the flow reactor,
as the time coordinate actually refers to the distance through the reactor.
In development of the plugflow reactor model, it is typically assumed
that the reactor is composed of infinitesimally thin ‘plugs’. Supposing that
a PFR with constant volume V has a constant volumetric flowrate Q, the
residence time of each reactor is given by:
τPFR =
V
Q
(5.19)
Approximating the PFR as a series of Nreac perfectlystirred reactors, the
residence time of each reactor (or plug) is straightforwardly calculated:
τPSR =
τPFR
Nreac
(5.20)
Under this approximation, the ability to model networks of flow reactors
with the stochastic particle method permits the transient response of the
population balance equations in a PFR to be determined. That is, both the
time and the reactor length coordinates can be treated independently. This
is relevant to reactor startup and operation.
For practical purposes, it is useful to know how many PSRs are needed
to accurately approximate the true PFR solution. This will be a function of
the relative timescales of the particle processes with respect to the transport
timescale. The Damköhler number is again used to measure the relative in
fluence of coagulation; however the characteristic transport time τ is here
taken to be the residence time of the PFR being approximated by the PSR
network (i.e. Da = τPFR/τcoag). To investigate the difference of a network
solution to the PFR solution, the ‘average PFR error’ PFR is first defined.
For a reactor network composed of Nreac PSRs, each with residence time
τPSR(Da); the error in the network’s solution with respect to the PFR solu
84
0
0.2
0.4
0.6
0.8
1
1.2
0 2 4 6 8 10 12 14 16 18 20
ε P
FR
(M
3),

Nreac, 
increasing coagulation
Da=0.1
Da=0.32
Da=1.0
Da=3.2
Da=10
Figure 5.9: Normalised average error in M3 of the networkedPSR solution
with respect to the PFR solution for particles undergoing coag
ulation (transition kernel).
tion is given by:
(Da, Nreac)PFR =
1
Nreac
Nreac∑
i
∣∣∣ξ(Da)Ri (tSS) – X*PFR (iτPSR)∣∣∣
X*PFR
(
iτPSR
) (5.21)
where ξRi is the mean of metric X(t) in reactor i of the network over L
runs (Equation (2.28)), and X*PFR the PFR solution for that metric. The
number of reactors required to approximate a PFR was studied for a system
of particles undergoing a physical coagulation process (transition kernel,
Equation (2.22)) with no inception or surface growth. The results for a
variety of Damköhler numbers are given in Figure 5.9.
It is clear that when the relative influence of coagulation is small (low
Da), the PSR network approximation to a PFR is valid for a small number of
reactors, and the magnitude of the error is very small. This approximation
appears to become less useful as the coagulation process occurs faster, rel
ative to the transport time. It is again illustrated in Figure 5.10, where the
zeroth moment is plotted for various (Da, Nreac) pairs.
85
0.65
0.7
0.75
0.8
0.85
0.9
0.95
1
0 0.2 0.4 0.6 0.8 1
M
0/M
0,
i,

t/τPFR, 
Da = 1.0
PFR
Nreac=5Nreac=10Nreac=20
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0 0.2 0.4 0.6 0.8 1
M
0/M
0,
i,

t/τPFR, 
Da = 10
Figure 5.10: Temporal evolution of M0 normalised by its initial value (M0,i)
of a reactor chain compared to the PFR solution.
5.4 Conclusions
This work has presented an algorithm to solve a reactor network with a
fullycoupled gasphase and stochastic population balance solver. The algo
rithm uses sequential modular simulation to solve the network in a stepwise
fashion. Novel processes for accounting for particle in and outflow were
proposed.
The convergence properties of the network solution algorithm were as
sessed in order to identify suitable choices for numerical parameters. A
choice of 4096 computational particles was recommended to minimise sys
tematic error. A new algorithm for simulating particle outflow (termed vol
ume ‘rescaling’) was also shown to be advantageous in reducing statistical
error. Two algorithms for coagulation were considered: the direct simula
tion algorithm and a stochastic weighted algorithm. It was found that the
SWA algorithm showed superior reduction of systematic and statistical er
ror for almost all properties of the particle ensemble. The direct simulation
algorithm was found to be unsuitable for networks with strong feedback.
86
To illustrate potential applications of this methodology, a plug flow re
actor was approximated as a linear series of perfectlystirred reactors. The
number of PSRs required to approximate such a system as a function of co
agulation rate was assessed. It was evident that many reactors are likely
to be necessary when particle processes occur on a faster timescale than
particle transport.
87
Chapter 6
A detailed model for silicon
nanoparticle synthesis
A novel multivariate particle model to simulate the synthesis of
silicon nanoparticles across a wide range of process conditions is
presented.
Material from this chapter has been published in [90].
Silicon nanoparticles were traditionally considered as an unwanted con
taminant formed as a byproduct of semiconductor manufacture [41, 104]
during the thermal decomposition of silane (SiH4). Today, they have a va
riety potential applications as a new material, particularly when made to
a specific size and shape [7, 58, 86, 102]. As particle properties are pri
marily controlled by the complex decomposition mechanism of silane, it is
necessary to understand the processes involved in this mechanism in order
to make particles of a given diameter and morphology. In this chapter, a
new model will be presented to simulate these processes leading to particle
formation.
The thermal decomposition of silane begins with the collision and ac
tivation of the silane molecule [63, 125, 129, 159]. It then breaks down
and recombines into other silanes, silenes (doublybound silicon hydrides,
e.g. SiH2SiH2, suffixed ‘A’ [149]) and silylenes (radical species, e.g. SiH3SiH,
suffixed ‘B’), gradually increasing in size until it is treated as a particle or
deposited on a surface.
88
This ‘clustering’ mechanism of gasphase silicon hydrides has been in
vestigated in a number of ways. For example, the decomposition of monosi
lane (and higher silanes) has been analysed using shock tubes [77, 159],
rotatingdisk reactors, [63] plasma discharge [154] and gas chromatogra
phy [16, 129]. Accompanying these studies are the numerical models of
the decomposition process. Many models of varying degree of complex
ity have been proposed, from the ‘simple’ threestep reaction of Purnell and
Walsh [129]; to the popular mechanism of Swihart and Girshick [149]; who
proposed a mechanism containing 2615 reactions and 221 silicon hydride
species, up to Si20. Quantum chemical calculations and automatic mecha
nism generation have also been used to further develop these mechanisms
[1, 160]. An excellent review of the current understanding of the silane
pyrolysis mechanism and in particular, its pressure dependence, is given by
Petersen and Crofton [125].
The heterogeneous reaction of silicon hydrides at silicon surfaces has
also been studied in detail. Houf et al. [64] identified that silanes react with
a silicon surface according to a dissociationabsorptionhydrogen mecha
nism, assuming that silenes and silylenes collide with surfaces with a stick
ing coefficient of unity. Further, it was reported that the rate of evolution of
hydrogen from silicon hydride surfaces is proportional to the coverage of hy
drogen at the surface [64, 147]. These mechanisms were incorporated into
the kinetic mechanism of Ho et al. [63] in a detailed model of the chemical
vapour deposition process of silicon.
Like the gasphase mechanism, silicon particle synthesis has been widely
investigated. Early studies focused on the process conditions required to ob
serve the onset of nucleation of particles; with a view to suppressing parti
cle formation so as to yield higherquality semiconductors [41, 101]. Since
then, several groups have used a variety of experimental methodologies to
study silicon particle synthesis. Cannon et al. [22] used CO2 laser heating to
analyse the decomposition of silane and subsequent formation of particles
in inert gases. A series of papers [23, 44, 95, 96] expanded the initial work
to encompass the effect of process conditions upon the properties of silicon
films deposited through laserheating of silane.
The decomposition of silane is also sensitive to specific microwave fre
quencies. This was used in several studies to form photoluminescent silicon
nanoparticles, or ‘quantum dots’ in a tubular microwave reactor [51, 58,
89
72]. It was reported that control of the microwave power (and therefore
the temperature profile) allowed for fine adjustment of the primary particle
size.
The Flagan group conducted impressivelydetailed experimental stud
ies into nucleation of silicon in single [103, 162, 163] and multistage
[135, 161] tubular hotwall reactors. In all cases, different heating zones
were used to tightly control the temperature profile of the reacting gas. De
pending on the reactor configuration and the imposed temperature profile;
particles of different size and morphology could be obtained. In typical reac
tor operation, it was reported that coagulation and sintering would control
particle growth unless care was taken to reduce the number concentration
of particles as produced by homogeneous nucleation [103].
These process conditions were recently identified by Körmer et al. [75],
who produced narrowlydistributed silicon nanoparticles in a hotwall tubu
lar reactor with dilute silane mixtures (∼ 4 %) at low pressures (2500 Pa)
and high temperatures (900–1100◦C ). Population balance modelling was
used to confirm that heterogeneous reaction of silicon hydrides at the parti
cle surface was predominantly responsible for particle growth [56, 76, 92].
There have been a number of other attempts to use population balance
modelling to understand silicon nanoparticle formation. Nguyen and Flagan
[103] used a fineseed model to investigate the effect of pressure and tem
perature upon particle structure. Kruis et al. [79] employed homogeneous
nucleation theory to elucidate the process conditions leading to onset of nu
cleation. This was expanded upon by Nijhawan et al. [104], who coupled
the gasphase mechanism of Swihart and Girshick [149] with a population
balance model for a spray reactor. The kinetic mechanism of Swihart and
Girshick [149] was also used in in a twodimensional fullycoupled model
incorporating the effect of fluid and particle dynamics [31, 32]; illustrat
ing the relative importance of the various interconnected particle processes
involved in synthesis.
Much work has been undertaken on the chemical and physical prop
erties of particles obtained from silane thermal decomposition. Onischuk
et al. [110] used a flow reactor to produce particles with specific hydro
gen content and chemical bond structure. Subsequent works used kinetic
and population balance modelling approaches to analyse kinetic pathways
through which particles were formed [111, 112, 114]. Aggregate particles
90
were studied in [113], in which particles were synthesised at low tempera
tures under finiterate sintering conditions.
Silane decomposition at aboveatmospheric pressure was addressed the
oretically and experimentally by Odden et al. [107], where it was reported
that particle size was independent of initial silane pressure. Additional
mechanistic pathways (in comparison to the mechanism of Ho et al. [63])
were also proposed to account for the observed reaction rates at these pro
cess conditions [108, 109].
The aerosol dynamics models described previously are typically ‘vali
dated’ by comparison of theoretical predictions of particle size distributions
(PSDs) against those obtained experimentally. These are often obtained via
direct transmission electron microscope (TEM) imaging [75], or, in older
work, though differential mobility analysers and condensation nuclei coun
ters [103, 162]. However, despite the wealth of information available about
experimental configurations, PSDs and mechanistic pathways, no effort has
been made to unify the conclusions of these previous studies.
The key aim of this chapter therefore is to develop a new model with a
detailed typespace to accurately predict the size and morphology of silicon
nanoparticles synthesised from the thermal decomposition of silane across
a range of different process conditions. The stochastic particle method cou
pled with a gasphase ordinary differential equation solver is employed to
achieve this.
The structure of this chapter is as follows. The gasphase and particle
models are explained in Section 6.1. The parameter estimation procedure
applied to the gasphase mechanism is outlined in Section 6.2. The nu
merical results are discussed in Section 6.3.1, while the model predictions
are compared with experimental data in Section 6.3.2. The chapter is con
cluded by an evaluation of the model and discussion of further avenues for
research.
6.1 Model development
6.1.1 Gasphase kinetic model
Much study has been carried out into the decomposition of silanes. For
conditions of practical interest, the rate of decomposition of silane can be
91
written as the bimolecular expression [125]:
SiH4 + M→ SiH2 + H2 + M (6.1)
where M is a third body. It is wellunderstood that silane decomposition
proceeds through a series of intermediate gasphase species, such as silylene
(SiH2) and higher silenes/silanes (e.g. Si3H8). Thus, silicon decomposition
mechanisms are often presented with many subsequent reactions after the
initial decomposition step. Such examples include the work of Ho et al. [63]
and Swihart and Girshick [149]. The presence of these appended ‘clustering’
reactions accelerates precursor decomposition [125].
This work adopts the mechanisms proposed by Ho et al. [63] to describe
the chemical kinetics involved in silane decomposition and is given in full
in Table 6.1. Note that the reactions for the species “Si” are omitted. Where
appropriate, the lowpressure limit parameters (according to a Lindemann
form falloff reaction, Equation 2.11) or reverse parameters are given. This
mechanism was chosen in place of more recent [1, 125] or more detailed
[49, 52, 149] kinetics due to sound testing at atmospheric pressure [63],
relative simplicity and provision of pressuredependent reactions.
The system of ordinary differential equations describing the reactions
of Table 6.1 was solved using MOPS’ ODE solver and the inhouse devel
oped chemical kinetics library Sprog for a zerodimensional batch reactor
configuration. In order to gain better agreement with experimental data,
preexponential factors for five reactions were adjusted. The reasons for
this are discussed in Section 6.2.
6.1.2 Particle model
Typespace
The binarytree model from Section 3.1 is used for the particle model in this
chapter. Each primary particle is described by the number of silicon atoms
ηSi and hydrogen atoms ηH:
px = px
(
ηSi, ηH
)
(6.2)
The derived properties of the typespace are identical to those given for
the binarytree particle model in Table 3.1. Here, the value of the fractal
92
Table 6.1: Gasphase mechanism employed. Parameters in bold have been
adjusted from their original referenced values (see Section 6.2).
Preexponentials are given with units in terms of cm, mol, s.
No. Reaction A n EA Ref.
cm, mol, s  kcal/mol
1 SiH4 (+M) −⇀↽ SiH2 + H2 (+M) 3.2× 109 1.7 54.7 [63]
lowpressure parameters 4.0× 1012 0 45.1 [125]
2 Si2H6 (+M) −⇀↽ SiH4+ SiH2 (+M) 1.8× 1010 1.7 50.2 [63]
lowpressure parameters 1.8× 1035 10.4 56.0 [63]
3 Si2H6 (+M) −⇀↽ Si2H4B + H2 (+M) 9.1× 109 1.8 50.2 [63]
lowpressure parameters 7.8× 1040 7.8 59.0 [63]
4 Si3H8 (+M) −⇀↽ SiH2 + Si2H6 (+M) 7.0× 1012 1.0 52.7 [63]
lowpressure parameters 1.7× 1069 15.1 60.5 [63]
5 Si3H8 (+M) −⇀↽ Si2H4B + SiH4 (+M) 3.7× 1012 1.0 50.9 [63]
lowpressure parameters 4.4× 1066 17.3 59.3 [63]
6 Si2H4B (+M) −⇀↽ Si2H4A (+M) 2.5× 1013 0.2 5.4 [63]
lowpressure parameters 1.1× 1033 5.8 9.2 [63]
7 Si2H4B + H2 −⇀↽ SiH4 + SiH2 9.4× 1013 0 4.1 [63]
reverse parameters 9.4× 1010 1.1 5.7 [63]
8 Si2H4B + SiH4 −⇀↽ Si2H6 + SiH2 1.7× 1014 0.4 8.9 [63]
reverse parameters 5.2× 1011 0.1 8.5 [63]
93
dimension (Df) is assumed to be 1.59, taken from the silicon particle syn
thesis work of Rogak et al. [135]. Their work also provides an expression to
estimate the mobility diameter of particles in the freemolecular regime:
dfmmob = dpriq
√
0.802(nq – 1) + 1 (6.3)
6.1.3 Particle processes
Particles may be created and changed according to a number of processes.
Inception
Particles are created in the population balance according to a general
inception reaction, where silylene can collide with a silene or another
silylene to form particles:
SiiHjB + SikHlX → Pq
(
p1(ηSi = i + k, ηH = j + l), C
(q)
)
(6.4)
where Pq(p1(ηSi, ηH), C(q)) is the (N(t) + 1)th particle, with ηSi silicon
atoms and ηH hydrogen atoms and X is A (silene) or B (silylene). This
inception mechanism is consistent with that proposed by Giunta et al.
[52]. The rate of inception is given by:
Rincep =
1
2
N2AcSiiHjBcSikHlX K
tr (SiiHjB, SikHlX) dq ≥ d*
0 dq < d*
(6.5)
where dq is the incepting diameter of the particle Pq resulting from
collision of species SiiHjB and SikHl. cSiiHjB and cSikHlX refer to the gas
phase concentrations of the colliding species and Ktr is the transition
regime coagulation kernel (given in Equation (2.22)). The variable
d* is the critical nucleus diameter, representing the threshold above
which gasphase species may convert to particles. It is calculated using
the Kelvin equation [76, 79]:
d* =
4 γSi v1
kB T ln(SSi)
. (6.6)
v1 is the volume of silicon monomer and γSi is the surface energy of
silicon, given by [76]:
γSi = 1.152 – 1.574× 10–4T(K) N/m (6.7)
94
SSi is the supersaturation of silicon, given by the ratio of silicon partial
pressure to saturation vapour pressure of silicon:
SSi =
pSi
pSi,sat
(6.8)
In the present work, the silicon partial pressure is approximated using
the sum of partial pressure of silylenes. The following correlation is
used for the saturation vapour pressure [55, 76]:
pSi,sat = 10
7.5341– 23399T(K) atm (6.9)
When the incepting diameter dq exceeds the critical nucleus diameter,
the rate of inception is calculated using the transition regime coagula
tion kernel.
Coagulation
Coagulation of particles results in the usual change to the typespace
for the binarytree particle model:
Pq
(
p1, . . . , pnq , C
(q)
)
+ Pr
(
p1, . . . , pnr , C
(r)
)
→ Ps
(
p1, . . . , pnq , pnq+1, . . . , pnq+nr , C
(s)
) (6.10)
The rate of coagulation is calculated using the transition regime coag
ulation kernel.
Surface reaction
This work adopts a version of the ‘dissociationadsorptionhydrogen’
mechanism, originally proposed by Sinniah et al. [147]. Silanes (SiH4,
Si2H6, Si3H8) react with the particle surface, releasing hydrogen. A
primary px of particle Pq is transformed as:
Pq
(
. . . , px(ηSi, ηH), . . . , C
(q)
)
+ SiiH2i+2
→ Pr
(
. . . , px(ηSi + i, ηH + 2), . . . , C
(r)
)
+ iH2
(6.11)
where C(r) is the new connectivity matrix resulting from the addition
of surface area. The change in connectivity matrix (C(q) → C(r)) is
given by Equation (4.8). The rate of surface reaction is modelled as
an Arrhenius process; and is also proportional to the surface area of
particle, Sq:
RSR = ASR,SiiH2i+2 Sq cSiiH2i+2 exp
(
–
EA,SR
R T
)
(6.12)
95
Table 6.2: Heterogeneous reaction processes included in the present work.
Parameters in bold have been adjusted from their original
referenced values (see Section 6.2). Surface reaction pre
exponentials have units cm/mol s and hydrogen release pre
exponentials have units 1/s and condensation scaling factors are
dimensionless.
Reaction Type A EA Ref.
cm, mol, s kcal/mol
SiH4 + px → px(ηSi + 1, ηH + 2) + H2 S.R. 3.0× 1033 37.5 [63]
SiH2 + px → px(ηSi + 1, ηH + 2) cond 1.0  [63]
Si2H6 + px → px(ηSi + 2, ηH + 2) + 2H2 S.R. 3.0× 1034 37.5 [63]
Si2H4A + px → px(ηSi + 2, ηH + 4) cond. 1.0  [63]
Si2H4B + px → px(ηSi + 2, ηH + 4) cond. 1.0  [63]
Si3H8 + px → px(ηSi + 3, ηH + 2) + 3H2 S.R. 3.0× 1034 37.5 [63]
px → px
(
ηSi, ηH – 2
)
+ H2 H2 1.2× 1019 47.0 [147]
where A and EA are the Arrhenius parameters for silane species SiiH2i+2.
The Arrhenius parameters are estimated by the fitting procedure, ex
plained in detail in Section 6.2. These reactions and parameters are
given in Table 6.2. The Arrhenius parameters for Si2H6 and Si3H8 were
set to ten times that of SiH4, consistent with the observations of Buss
et al. [21] and as done in previous models [63, 104].
Condensation
Particles can grow by the deposition of silenes and silylenes on the
particle surface. The statespace of the particle is modified according
to:
Pq
(
. . . , px(ηSi, ηH), . . . , C
(q)
)
+ SiiHjX
→ Pr
(
. . . , px(ηSi + i, ηH + j), . . . , C
(r)
) (6.13)
This may be modelled as a collisional process, with rate given by:
Rcond = Acond cSiiHjX
√
pi kB T
2 mSiiHjX
(
dSiiHjX + dcol,q
)2
(6.14)
where Acond is the collisional efficiency of the process (assumed 1.0,
consistent with [63] for silylenes), mSiiHjX and dSiiHjX are the mass and
diameter of the colliding species respectively, and dcol,q is the collision
96
diameter of particle q. These processes are again summarised in Table
6.2.
Hydrogen release
The hydrogen content of a silicon nanoparticle is a strong function of
temperature and annealing time [95, 96, 147]. The present model
therefore includes a stochastic jump process to account for release of
hydrogen from particles:
Pq
(
. . . , px(ηSi, ηH), . . . , C
(q)
)
→ Pr
(
. . . , px(ηSi, ηH – 2), . . . , C
(r)
)
+ H2
(6.15)
A study of hydrogen desorption from a silicon hydride film found that
the rate of desorption was firstorder in coverage of hydrogen [147].
Here, the hydrogen coverage θ is estimated by the ratio of number of
hydrogen atoms to silicon atoms in each particle:
θ(Pq) =
∑nq
x=1 px(ηH)∑nq
x=1 px(ηSi)
(6.16)
It is assumed that the common surface element remains unchanged
due to the loss of hydrogen from the particle surface (C(q) = C(r)).
The rate of hydrogen desorption is proportional to the coverage, given
here by the ratio of number of hydrogen to number of silicon atoms:
RH2 = AH2 θq exp
(
–
EA,H2
R T
)
(6.17)
The preexponential parameter AH2 is adjusted according to the method
ology in Section 6.2, and EA,H2 is given in Table 6.2
Sintering
Sintering is implemented using the exponential excess surface area
decay formula given in Equation (3.18). In this model, the sintering
parameters suggested by Kruis et al. [78] are used. These parameters
account for sintering through a grainboundary diffusion process, with
the characteristic time given by:
τS = AS d
nS T exp
(
ES
R T
)
(6.18)
where AS, nS and ES are empirical parameters.
97
Table 6.3: Numerical and model parameters used in the silicon model.
Description Symbol Value Ref.
Numerical parameters
Maximum splitting timestep ∆ts 2.5× 10–4 s [143]
Max. number of stochastic particles Nmax 16,384 [143]
Number of runs L 8 [143]
Model parameters
Maximum zeroth moment M0,max 1.0× 1014 #/m3 
Sintering preexponential AS 1.15× 1013 s/m4 [78]
Sintering diameter power nS 4 [78]
Sintering characteristic energy ES 55.0 kcal/mol [78]
Surface smoothing factor σSR 1.0 [144]
Density of silicon ρSi 2329 kg/m
3 
Density of hydrogen ρH 101 kg/m3 
Freemolecular enhancement factor kfm 1.0 
Process settings
Temperature range T 800–1400 ◦C 
Residence time range τ 0.001–1.0 s 
Initial fraction of SiH4 range ySiH4 0.0004–1.0 
Total pressure range P 2500–101,325 Pa 
6.1.4 Coupling algorithm
The gasphase kinetic model and the particle population balance are cou
pled using an operatorsplitting algorithm with Strang splitting [25, 26].
The Linear Process Deferment Algorithm is used to accelerate linear pro
cesses such as surface reaction and condensation (see Section 2.2.2). The
full algorithm though which this occurs for an analogous particle model is
described by Shekar et al. [143]. The numerical and model parameters used
for the fullycoupled model are summarised in Table 6.3.
98
6.2 Parameter estimation
Previous studies have identified that, at low pressure, use of literature gas
phase mechanisms (e.g. that of Swihart and Girshick [149]) requires scaling
of either the gasphase preexponentials [56, 76, 92] or scaling of the nu
cleation rate [104]. In some cases, both the kinetic parameters and the nu
cleation rates (via usage of collision efficiencies) were empirically adjusted
to obtain model fits [113, 114]. Thus, it is necessary to choose appropriate
parameter values to ensure theoretical results are comparable with experi
mental results.
Of the reactions included in the gasphase mechanisms, Petersen and
Crofton [125] noted that the silane and silylene concentration was most
sensitive to Reactions 1 and 7 (given in Table 6.1). Onischuk et al. [114]
additionally reported that fitting of Reactions 1 and 3 was required for their
population balance model. Preliminary calculations indicated that adjust
ment of the preexponentials of Reaction 5 and Reaction 8 (reverse) was
necessary to achieve agreement with experimental work. As silane decom
position occurs at the lowpressure limit [125] up to 5 atm the lowpressure
preexponentials were adjusted, where available.
As in previous studies, the surface reaction [56, 76, 92] and hydrogen
release [114] rates were adjusted via the preexponential factors. In total,
there are seven parameters which require estimation. This system’s param
eter vector is therefore given by:
x =
(
A1,LP, A2,LP, A3,LP, A5,LP, A8,rev, ASR,SiH4 , AH2
)
(6.19)
To choose an appropriate set of parameters x, a staged optimisation
methodology analogous to those given in [92] was employed. Firstly, low
discrepancy (Sobol) sequences were used to located four parameter sets
close to minima, as defined by objective function Φ(x):
Φ(x) =
Nexp∑
i=1
(
φ
exp
i – φ
sim
i (x)
)2
(6.20)
where φi represents one of the datasets used in fitting, listed in Table
6.4. The superscript ‘exp’ refers to an experimentallysourced dataset, and
‘sim’ that obtained from simulation. A selection of studies using a range of
99
Table 6.4: Experimental data used in fitting the model. σa(dpri) refers to
the arithmetic standard deviation of the primary diameter.
i T ySiH4 P τ φi φ
exp
i , nm Ref.
1 1100 4.0 2.5 0.08 dpri,mode 26.7 [75]
2 1100 4.0 2.5 0.08 σa(dpri) 1.7 [75]
3 1100 4.0 2.5 0.08 dpri,10 24.1 [75]
4 1100 4.0 2.5 0.08 dpri,90 28.8 [75]
5 900 4.0 2.5 0.42 dpri,mode 21.2 [75]
6 580 5.0 39.0 0.53 dpri,mode 41.0 [113]
7 816 3.3 51.0 0.0026 dpri,mode 11.0 [49]
8 1047 3.3 51.0 0.0021 dpri,mode 11.0 [49]
process conditions were used, and equal weighting was applied to each of
the datasets used in the objective function. The bounds for the Sobol scans
were found empirically. The four best parameter sets from the initial param
eter space scan are further refined by applying the simultaneous stochastic
approximation algorithm (SPSA) [62]; the best of which x* is chosen ac
cording to:
x* = arg min
x
{Φ(x)} (6.21)
The optimal set of parameters is obtained via use of the Bayesian ap
proach of Mosbach et al. [100] and Braumann et al. [17]. Using an Inverse
Wishart prior (as the error in the experimental values is unknown), a Metropolis
Hastings algorithm is employed to estimate the posterior distribution of the
parameters. A fourthorder response surface was constructed locally around
the point x* in order to reduce the time needed for the many model evalu
ations required by the sampling algorithm. This methodology provides an
estimate of the modal value of each parameter in x and its corresponding
uncertainty.
6.3 Results and discussion
6.3.1 Numerical results
The model parameters and their credible intervals as obtained from the pa
rameter estimation procedure are given in Table 6.5. A qualitative under
100
Table 6.5: Optimal parameters obtained from the parameter estimation
procedure. The mode, upper and lower bounds are those of the
high probability density region for a 99.99% credible interval.
Parameter Units Lower bound Mode Upper bound
A1,LP cm3/mol.s 3.2× 1012 4.0× 1012 5.8× 1012
A2,LP cm3/mol.s 3.8× 1033 1.8× 1035 2.8× 1037
A3,LP cm3/mol.s 2.1× 1039 7.8× 1040 1.8× 1043
A5,LP cm3/mol.s 1.7× 1054 8.4× 1055 1.4× 1058
A8,rev cm3/mol.s 8.1× 1010 5.2× 1011 7.0× 1013
ASR,SiH4 cm/mol.s 1.1× 1033 3.0× 1033 1.3× 1034
AH2 1/s 1.2× 1019 3.08× 1020 6.3× 1020
standing of sensitivity of the objective function Φ(x) to each parameter can
be gained from this information. For example, the credible intervals of A1,LP
and ASR,SiH4 are within an order of magnitude of the modal value, indicating
that the accuracy of the model is strongly dependent on correct estimation
of these values. In comparison, the bounds for A3,LP span four orders of
magnitude, suggesting that Φ(x) is less sensitive to this value.
As the objective function is very sensitive to Reaction 1 and previous
studies [125] have identified it as the most important step in governing the
decomposition rate, it is worthwhile to discuss the value of A1,LP obtained
here through modelfitting. Table 6.6 lists the lowpressure parameters for
Reaction 1 for several other studies: one used in a mechanism to describe
chemical vapour deposition (CVD) of silane [63], another obtained through
detailed experimental studies [125], and the others used for coupling with
particle models [49, 76, 92, 123].
There is clearly a large variation in values of the preexponential, espe
cially when the additional temperature fitting power n is used. It is useful to
use a unimolecular rate constant (k1b) when comparing the reaction rates
predicted by each set of parameters [125]. Here, it is given by:
dcSiH4
dt
= –k1cSiH4cM = –k1bcSiH4 (6.22)
This rate constant is compared across the kinetic models of Table 6.6 in
Figure 6.1. It is evident that the rate constants in pure gasphase models
(e.g. [63, 125]) are three to four orders of magnitude greater than those
101
Table 6.6: Comparison of other kinetic expressions for the lowpressure
limit of Reaction 1.
Study A n EA
cm3/mol.s  kcal/mol
Ho et al. [63] 5.2× 1029 3.5 57.6
Frenklach et al. [49] 1.2× 1047 8.8 67.2
Petersen and Crofton [125] 7.2× 1015 0 45.1
Paur et al. [123] 1.3× 1012 0.68 58.0
Körmer et al. [76] 9.0× 1012 0 45.1
Menz et al. [92] 5.1× 1010 0 45.1
this work 4.0× 1012 0 45.1
coupled to particle models [49, 76, 92, 123]. Even within the set of those
coupled to particle models, the unimolecular rate constant spans three or
ders of magnitude. The rate expression found in this thesis (from parameter
estimation) appears to correspond well with that of Körmer et al. [76].
It is hypothesised that k1b reported here is a result of fitting an overall
inception rate to the experimental case studies. As already discussed, mod
els developed for silicon nanoparticle synthesis typically require scaling of
the nucleation rate (with a collision efficiency) or silane decomposition rate
(as in this case). This clearly indicates that a key challenge is obtaining the
correct nucleation rate of ‘seed’ particles. It is therefore hypothesised that
when combined with the Equation (6.5), the gasphase mechanism found
through fitting (Table 6.1) provides a correct overall nucleation rate for the
experimental cases studies. On its own, the gasphase mechanism is not
physically representative of the real decomposition rate.
After conversion from SURFACE CHEMKIN [30] to the typespace of this
particle model, the value of ASR,SiH4 shows good agreement with that pro
posed by Houf et al. [64] and used in the gassurface mechanism of Ho et al.
[63]. The conversion is given in Appendix B.1.
6.3.2 Case studies
There is an overwhelming amount of literature which describes the produc
tion of silicon nanoparticles under a variety of process conditions. Several
of these cases representing different types of particles or reactor configu
102
10−6
10−4
10−2
100
102
104
106
108
0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1
k 1
b,
1/
s
1/temperature, 1000/K
Ho et al., 1994
Frenklach et al., 1996
Petersen + Crofton, 2003
Paur et al, 2005
Körmer et al., 2010
Menz et al., 2012
this work
Figure 6.1: Dependence of the unimolecular rate constant (k1b) upon tem
perature at 2.5 kPa (conditions of Körmer et al. [75]) for a va
riety of literature kinetic models.
rations were selected. Those studied in the present work are displayed
in Table 6.7. Where available, published temperature profiles were used
(e.g. [163]), given in Figure 6.2.
All profiles (e.g. temperature versus length in [45]) were converted from
length coordinates to time coordinates assuming plug flow and volumetric
flowrates, unless a residence time was supplied. In the cases with no tem
perature profile (e.g. [113]), the time corresponds to a hotzone residence
time; and a constant temperature was used. PSDs are the numberdensity
distributions and were generated using kernel density estimation and the
normal distribution approximation to estimate the bandwidth [146].
Narrow PSDs
Narrowlydistributed particles are obtained when inception occurs for a very
short period of time, coagulation is negligible, and particles are grown by
heterogeneous deposition [76, 92]. A detailed study of the production of
uniform silicon nanoparticles was conducted by Körmer et al. [75], where
particle size could be controlled between 20 and 40 nm through choice of
103
Table 6.7: Experimental data and corresponding process conditions inves
tigated for model comparison.
Study Reactor T ySiH4 P (bath) τ
◦C % kPa s
Flint et al. [45] Laser 700–1000 0.21 20 (Ar) 0.005
Wu et al. [163] PFR 450–1200 1.0 101 (N2) 1.0
Nguyen and Flagan [103] PFR 500–1200 0.04–1.0 101 (N2) 0.9
Frenklach et al. [49] Shock tube 816–1307 3.3 51 (Ar) ∼0.001
Onischuk et al. [113] PFR 580 5.0 39 (Ar) 0.2–0.8
Körmer et al. [75] PFR 500–1100 2.0–13.0 2.5 (Ar) 0.08–0.42
200
400
600
800
1000
1200
1400
0.0001 0.001 0.01 0.1 1
te
m
pe
ra
tu
re
, o
C
time, s
Flint et al., 1986
Wu + Nguyen, 1987
Nguyen + Flagan, 1991
Körmer et al., 2010
Figure 6.2: Experimental temperature profiles used in the model calcula
tions for the case studies of Table 6.7.
104
process conditions.
Under these conditions, much of the precursor is lost to the wall due
to vapour deposition on the reactor surface. In a similar vein to [56], this
work assumes that only 16% of the initial precursor reacts to form parti
cles. A temperature profile similar to the optimised profile of Gröschel et al.
[56] is used for the 84 ms, 1100◦C cases, given in Figure 6.2. The model
predictions and experimental results for this system are compared in Figure
6.3.
As previously identified [56, 76, 92], coagulation should be negligible in
order to obtain a narrow size distribution. There is reasonable agreement
between the model and literature results, especially for greater silane frac
tions. The secondary and very small tertiary peaks for the ‘base case’ is a
result of particles coagulating once or twice. For the higher silane fraction
cases, much more silicon in the particle phase yields many more coagula
tion events, broadening the distribution and ‘smoothingout’ the additional
modal peaks.
At long residence times (τ ∼ 420 ms) the PSDs predicted by the model
are far too broad as compared to those observed experimentally [75]. It
is evident that in these cases, the model either overpredicts the rate of
coagulation or is not capturing some aspect of detail of the real process. It
is likely that transport processes should be included in the model for this
case study, as there is much precursor lost to the walls, particle size and
morphology has a radial dependency [75] and the system is particularly
sensitive to the temperature profile imposed [56].
Nucleation of silicon particles in a shock tube was studied by Frenklach
et al. [49]. Spherical particles with sizes between 10 and 40 nm and narrow
size distributions were produced in the shock tube at a range of process con
ditions. Experimental PSDs were reported for a 3.3% SiH4Ar mixture with
a total pressure of 0.5 atm, and are compared to the theoretical predictions
in Figure 6.4.
Very good agreement of theoretical and experimental PSDs is observed
for the 816◦C and 1307◦C cases. However, the modal diameter for the
1047◦C case is approximately double that reported in [49]. In a similar
vein to the 420 ms cases of Körmer et al. [75] described previously, the
broadening observed in this case is associated with too much coagulation.
105
0
0.05
0.1
0.15
0.2
0.25
0.3
15 20 25 30 35 40 45
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
base case: τ=84ms, T=1100°C, ySiH4=4%
0
0.05
0.1
0.15
0.2
0.25
0.3
10 20 30 40 50 60
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=84ms, T=1100°C, ySiH4=8%
0
0.05
0.1
0.15
0.2
0.25
0.3
10 20 30 40 50 60
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=84ms, T=1100°C, ySiH4=12.8%
0
0.05
0.1
0.15
0.2
0.25
0.3
10 15 20 25 30 35 40 45 50
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=420ms, T=1100°C, ySiH4=4%
exp.
sim. dpri
0
0.05
0.1
0.15
0.2
0.25
0.3
10 15 20 25 30 35 40 45 50
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=420ms, T=1000°C, ySiH4=4%
0
0.05
0.1
0.15
0.2
0.25
0.3
10 15 20 25 30 35 40 45 50
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=420ms, T=900°C, ySiH4=4%
Figure 6.3: Comparison of the experimental PSDs of Körmer et al. [75]
(crosses) at 25 kPa total pressure (in Ar) and those predicted
by the present model (lines).
106
0
0.05
0.1
0.15
0.2
0 10 20 30 40 50
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=1.6 ms, T=816 °C
0
0.05
0.1
0.15
0.2
0 10 20 30 40 50
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=2.6 ms, T=816 °C
exp. dpri
sim. dpri
0
0.05
0.1
0.15
0.2
0 10 20 30 40 50
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=1.2 ms, T=1047 °C
0
0.05
0.1
0.15
0.2
0 10 20 30 40 50
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=2.1 ms, T=1047 °C
0
0.05
0.1
0.15
0.2
0 10 20 30 40 50
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=1.0 ms, T=1307 °C
0
0.05
0.1
0.15
0.2
0 10 20 30 40 50
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
τ=1.8 ms, T=1307 °C
Figure 6.4: Comparison of theoretical (lines) and experimental PSDs
(crosses) for the shock tube reactor of Frenklach et al. [49] with
initial conditions ySiH4=3.3% (in Ar) and P=51 kPa.
107
0
0.001
0.002
0.003
0.004
0.005
0.006
0.007
0.008
0.009
0 100 200 300 400 500 600 700 800
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
exp. dmob
sim. dmob
Figure 6.5: Comparison of theoretical (lines) and experimental PSDs
(crosses) for the system of Wu et al. [163] with initial condi
tions ySiH4 = 1.0% (in N2) at atmospheric pressure.
Sintered particles
Wu et al. [163] used a tubular reactor to generate highlyspherical crys
talline silicon particles through coagulation and sintering of smaller parti
cles at a high temperature. A peak temperature of 1200◦C (temperature
profile given in Figure 6.2) and initial silane fraction of 1% at atmospheric
pressure (in nitrogen) was required to achieve these particles. The resulting
particle size distribution is given in Figure 6.5.
The model predictions show excellent agreement with the experimental
PSD. The particles predicted by the model are completely spherical, with a
geometric standard deviation of 1.50, indicating that the particle size distri
bution is at the selfpreserving limit.
Particle aggregates
Particle synthesis under conditions of finiterate sintering was experimen
tally and theoretically investigated by Nguyen and Flagan [103]. Using a
singlestage tubular reactor, they investigated formation of aggregates at
different silane concentrations (0.01–1.0%) with a lower peak reactor tem
108
0
20
40
60
80
100
120
140
160
180
200
0.01 0.1 1
di
am
et
er
, n
m
initial ySiH4, %
exp. dmob
sim. dmob
exp. dpri
sim. dpri
Figure 6.6: Dependence of particle size on initial precursor concentration
for the system of Nguyen and Flagan [103] with a peak reactor
temperature of 800◦C .
perature as compared to Wu et al. [163]. This reactor was operated at atmo
spheric pressure with nitrogen as the bath gas. Their results are compared
to those predicted by the present model in Figures 6.6 and 6.7.
It is evident that the theoreticallypredicted average mobility diameter
and the mobility diameter particle size distributions (Figure 6.7) are in good
agreement with the experimental data. However, the primary diameter is
too large for silane concentrations above 0.07% (Figure 6.6). This is at
tributed to the sintering kinetics used; which remain untested in population
balances [92].
The formation of particle aggregates was also identified by Onischuk
et al. [113] in a detailed examination of the nature of particle structures
obtained when sampling from the hot and coldzones of the reactor. Ag
gregates were sampled from inside the reactor, corresponding to various
hotzone residence times and their primary particle structure analysed. The
model and experimental results are compared in Figure 6.8 for a constant
temperature of 580◦C and initial conditions of 5% silane in argon at 39 kPa.
Excellent agreement is observed between the model predictions of the
temporal evolution of the average primary diameter and the experimental
109
0
0.005
0.01
0.015
0.02
0 50 100 150 200 250 300
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
ySiH4=0.10%
exp. dmob
sim. dmob
0
0.005
0.01
0.015
0.02
0 50 100 150 200
ke
rn
el
d
en
sit
y,
1
/n
m
diameter, nm
ySiH4=0.04%
exp. dmob
sim. dmob
Figure 6.7: Comparison of theoretical (lines) and experimental PSDs
(crosses) for the conditions of Nguyen and Flagan [103] with a
peak reactor temperature of 800◦C .
0
10
20
30
40
50
60
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
di
am
et
er
, n
m
time, s
exp. dpri
sim. dpri
Figure 6.8: Temporal evolution of the mean primary particle diameter for
the conditions of Onischuk et al. [113] with initial conditions
ySiH4=5% (in Ar) at 39 kPa and 580
◦C .
110
experiment simulation
Figure 6.9: An experimental TEM micrograph from the work of Onischuk
et al. [113] (left panel) compared with a computergenerated
TEMstyle image from the present work’s model (right panel)
for 0.53 s reaction time followed by 17.8 s in the reactor ‘cold
zone’. Particle images were rendered using POVRay [124].
results. In order to compare particle aggregates, a TEMstyle image was
generated from model results by assuming that primary particles randomly
adhere to each other (Figure 6.9). Here, it is clear that there is at least
qualitative agreement with aggregate size and structure.
Laser synthesis of particles
Silicon particles can also be synthesised in a laserdriven flame [22, 23,
44, 45]. These studies varied process conditions such as silane fraction and
cell pressure to produce highly spherical and looselyagglomerated particles.
The particle size and number density was tracked through the flame using
laserscattering measurements. To convert the temperaturelength profiles
given in [45], the velocity of the particles was estimated using the formula
provided in [44]:
v =
2 QSTP T/273
pi r2 P(atm)
(6.23)
111
0
10
20
30
40
50
60
0 0.001 0.002 0.003 0.004 0.005 0.006
di
am
et
er
, n
m
time, s
exp.
sim. dcol
sim. dpri
Figure 6.10: Evolution of mean particle diameter for the laserdriven flame
(case 630S) of Flint et al. [45] and Cannon et al. [23].
where QSTP is the volumetric silane flow at standard temperature and pres
sure and r is the radius of the flame. Thus, using the velocity through the
flame, the temperaturetime profiles were obtained. The experimental data
and model predictions for this system are supplied in Figure 6.10 for the
630S powder (see [45] for more information about different powder codes),
assuming 109 cm3/s of silane [45] and 400 cm3/s of argon [22] in the inlet
stream.
The new model gives reasonable agreement with the experimental tra
jectories when comparing the primary particle diameter to the experimen
tally observed diameter. The separation of dcol and dpri at approximately
0.003 s represents the peak of the temperature profile, where the reaction
has gone to completion and aggregates begin to form.
6.4 Conclusions
A novel fullycoupled gasphase and particle model for the synthesis of sili
con nanoparticles was presented. The mechanism of Ho et al. [63] was used
to describe the gasphase kinetics. A multidimensional particle model track
112
ing the number of silicon and hydrogen units of each primary particle and
the connectivity of the primaries in an aggregate was used in conjunction
with a population balance model to study nanoparticle synthesis.
By adjusting the preexponential factors of five gasphase and two het
erogeneous reactions, good fit to a variety of experimental studies was ob
tained. The model was applied to different reactor configurations, process
conditions and particle morphologies. Varying degrees of fit of theoreti
cal results with experimental results were obtained, with the model gener
ally better for atmospheric pressure cases. The robustness of this modelling
methodology was demonstrated by generating a TEMstyle image of silicon
aggregates in qualitative agreement with experiment.
It was hypothesised that, when combined with the inception process, the
new gasphase mechanism with adjusted parameters accurately described
the silicon particle inception rate. It was, however, noted that gasphase
mechanism is unlikely to represent the physical rate of decomposition of
silicon hydrides. A key challenge in developing a physical model for silicon
nanoparticle formation is resolving this issue.
It was reported that in some cases where aggregate particles were pro
duced, or where large amounts of coagulation occurred, the primary diam
eter was larger than that observed experimentally. This was attributed to
the finiterate sintering kinetics employed. As this area is generally poorly
characterised (particularly in terms of experimental work) for silicon, it is
suggested that the model could be significantly improved by using an accu
rate sintering expression. Despite this, the model successfully knitstogether
many of the historical studies of silicon nanoparticle synthesis and provides
a foundation on which the space of process conditions can be reliably ex
plored.
113
Chapter 7
Putting it all together
The new methods developed in Chapters 3–5 are applied to the
experimental case Wu et al. [163] with the new silicon particle
mechanism.
Material from this chapter has been published in [91, 94].
The advances in modelling presented in Chapters 3–5 have so far only
been applied to a general population balance model, or for the silica par
ticle model of Shekar et al. [144]. In this chapter, the methods developed
in Chapters 3–5 will be applied to the new model for silicon nanoparticle
synthesis presented in Chapter 6. It is intended that this should illustrate
potential uses for such methods.
The experimental system of Wu et al. [163] will be investigated using
the new model and methods. This case was selected as it represented a
relatively simple hotwall reactor system and showed good agreement with
experimental results (Figure 6.5). In Section 7.1, the case will be analysed
with different particle models, while in Section 7.2, a linear network of
reactors will be used to investigate the unsteadystate response of the hot
wall reactor.
7.1 Particle model analysis
In Chapter 3, the use of different particle models was investigated. It was
identified that use of a particular particle model can strongly alter the solu
tion of the population balance model. The calculations done in that analysis
were for a generic implementation of the binarytree particle model, and
114
did not include factors such as a reactive gasphase or surface growth on
particles.
Here, the analysis of Chapter 3 will be expanded by combining the spher
ical and surfacevolume particle models with the new mechanism for silicon
nanoparticle formation given in Chapter 6. The key aim is to assess how
different particle models behave in a ‘reallife’ modelling application.
To do this, characteristic times for coagulation and sintering of the par
ticle ensemble are first defined. The characteristic coagulation time is esti
mated using the average properties of the ensemble:
τC(t) =
1
M0(t)Ktr
(
µa
[
dcol
]
(t),µa [m] (t)
) (7.1)
where the average mass and average diameter are substituted into Equation
(2.22) for Ktr. The expression of Equation (6.18) for sintering is nonlinear.
Thus, the three particle models will begin to diverge when sintering moves
from effectively instantaneous (very small particles) to finiterate. This tran
sition is best approximated by use of the maximum primary particle diame
ter, hence the characteristic sintering time is given by:
τS(t) = τS
(
max
[
dpri
]
(t), T(t)
)
(7.2)
The evolution of important particle and ensemble characteristics for this
model, as applied to the case of Wu et al. [163], is given in Figure 7.1. Note
that a temperature profile (increasing from 500◦C to 1250◦C ) is imposed
across this simulation (Figure 6.2).
All the models predict that spherical particles are obtained at the end
of the process time. However, neither the surfacevolume nor the spherical
particle model perform comparably to the binarytree model. The number
density (M0) is double the binarytree model’s prediction, and the average
collision diameter (µa(dcol)) 20 % lower. The surfacevolume and spherical
models predict the occurrence of a secondary nucleation period (at t = 0.6s
in M0 panel of Figure 7.1), consistent with other spherical particle modelling
efforts [103] of this system.
The predictions of the surfacevolume model are also very similar to
those of the spherical particle model, despite containing a finite sintering
rate. It is hypothesised that this is due to the high polydispersity of particles
(σg > 3.0, Figure 7.1 lower panels) for this case study, resulting in a similar
115
1014
1015
1016
1017
1018
0 0.2 0.4 0.6 0.8 1
M
0,
1/
m
3
time, s
bintree
surfvol
spherical
109
108
107
106
105
0 0.2 0.4 0.6 0.8 1
F V
,

time, s
0
50
100
150
200
250
0 0.2 0.4 0.6 0.8 1
µ a
(d c
o
l),
nm
time, s
0
20
40
60
80
100
120
140
160
0 0.2 0.4 0.6 0.8 1
µ a
(d p
ri),
nm
time, s
1.0
2.0
3.0
4.0
5.0
6.0
7.0
0 0.2 0.4 0.6 0.8 1
σ
g(d
co
l),

time, s
1.0
2.0
3.0
4.0
5.0
6.0
7.0
0 0.2 0.4 0.6 0.8 1
σ
g(d
pr
i),

time, s
bintree µa[σg(dpri)]
Figure 7.1: Comparison of model predictions for the case of Wu et al. [163]
using Chapter 6’s silicon model. Depicted are the zeroth mo
ment (M0), volume fraction (FV), arithmetic mean (µa) and
geometric standard deviation (σg) of the collision and primary
diameters.
116
t/τ
C,

t/τS, 
5 % error
bintree
surfvol
spherical
t=0.40s
t=0.75s
t=1.0s
t=1.0s
t=0.75s
0
2
4
6
8
10
12
14
102 101 100 101 102 103
direction
of time
Spherical
zone
Nearspherical
zonePartialsintering zone
Figure 7.2: Trajectory of the case of Wu et al. [163] using Chapter 6’s silicon
model through the t/τC, t/τS space. The shading depicts the 5%
error margins as described in Figure 3.3. Specific process times
for the binarytree and surfacevolume models from Figure 7.1
are labelled, t < 0.2s omitted for clarity.
‘convergence’ of the two models as depicted in Figure 3.4. To further under
stand the differences between these models, this case is plotted atop Figure
3.3’s coagulationsintering map in Figure 7.2.
The system begins with small spherical particles at high t/τS. The three
models diverge from each other after passing through the nearspherical
zone. The surfacevolume and spherical models subsequently show a large
increase in coagulation rate, while it steadily decreases for the binarytree
model. This is due to the secondary nucleation period causing the number
density to remain steady (Figure 7.1). The secondary nucleation period
is absent for the binary tree model, as the increased surface area of the
aggregate particles causes faster depletion of the precursor.
The minimum sintering rate, reached at approximately 0.75 s, repre
sents the point at which the sintering time’s rate of decrease is balanced by
the contributions from the rate of increase in particle diameter and rate of
117
increase in temperature. Past this point, the increasing temperature pro
file causes coalescence of aggregate structures, returning all models to the
spherical zone.
Finally, it is also important to note that the history of the particle ensem
ble through the t/τC, t/τS space is as important as its location in the space
(Figure 7.2). When the particle system moves from one zone to another, it
does not imply that all particles have the structure represented by the zone
type. Rather, it implies that should the system remain in that zone, as the
ensemble’s particles will change with time to eventually resemble that of the
zone.
7.2 Unsteadystate response
In Chapter 5, the reactor network infrastructure was used to approximate
a PFR by consecutive PSRs. A SWA was used with the new formulation of
the transition regime coagulation kernel developed in Chapter 4 to achieve
this with reasonable precision. The same reactor network model will now
be used to simulate the unsteadystate response of the hotwall PFR of Wu
et al. [163].
A network of 30 consecutive PSRs was composed, assuming that each
PSR is isothermal and has a temperature which corresponds to its location
along the PFR (i.e. Ti = T(iτPSR)). This was done to account for the in
creasing temperature profile through the PFR (Figure 6.2). The initial state
of the network is taken to have no particles, but is filled with nitrogen. The
transient response of the PSR network was evaluated, allowing sufficient
time for the system to reach steady state. The results of this analysis are
given in Figure 7.3, where the transient and axial evolution of the particle
properties are plotted.
Figure 7.3 (left panel) shows the behaviour of the primary polydisper
sity, given by the average of the geometric standard deviation of primary
diameters in a particle. The ‘time through reactor’ coordinate refers to the
axial position in the reactor, whereas the ‘time’ axis is the transient response.
There is good agreement between the network and PFR solutions in early
reactor times (i.e. initial section of PFR length). This agreement appears to
decay as coagulation accelerates, similar to Figures 5.9 and 5.10.
The degree of sintering (Figure 7.3, right panel) shows better agreement
118
time through reactor, s
(axial position)
0.0
0.2
0.4
0.6
0.8
1.0
tran
sie
nt r
esp
ons
e ti
me
, s
0.0
0.2
0.4
0.6
0.8
1.0
1.2
pr
im
ar
y
po
ly
di
sp
er
si
ty
, 
1.0
1.2
1.4
1.6
1.8
2.0
PFR Solution
time through reactor, s
(axial position)
0.0
0.2
0.4
0.6
0.8
1.0
tran
sie
nt r
esp
ons
e ti
me
, s
0.0
0.2
0.4
0.6
0.8
1.0
1.2
av
g.
si
nt
er
in
g
le
ve
l,

0.0
0.2
0.4
0.6
0.8
1.0
Figure 7.3: Mesh plots depicting the transient evolution of the primary
polydispersity (given here by the Equation (3.21)) and the de
gree of sintering. Contour projections of the network solution
are compared to the PFR solution (thick line).
119
AB
C
D
A
B
C
D
A B C D
200 nm200 nm10 nm10 nm
Figure 7.4: Contour plots of transient evolution of average collision diame
ter and particle number concentration. Insets illustrate the par
ticle ensemble state with TEMstyle images at the given reactor
time (i.e. axial position) and transient response time coordi
nates, generated using POVRay [124].
than the primary polydispersity with PFR solution. In this particular exam
ple, the growth by coagulation and sintering under similar timescales [91]
causes aggregates of partiallysintered primaries to form, thus depression in
the sintering level in both time trajectories is observed.
Information stored in the binarytree particle model can be used to gen
erate images of sample particle structures. This is illustrated in Figure 7.4,
where sample TEMstyle images are generated for specific PFR axial position
times and transient response times.
The collision diameter (Figure 7.4 left panel) and number concentra
tion (Figure 7.4 right panel) are presented to give representative values of
120
particle size and quantity as a function of reactor position and transient re
sponse time. Panel A represents the initial generation of ‘seed’ particles,
which are transported through the network until point B, where the tem
perature profile causes further rapid nucleation and surface growth. The
spike in number concentration at B increases the coagulation rate, causing
aggregates to form (Panel D). If particles are transported to the outlet of the
reactor, the final temperature increase sinters the particles back to spheres
(Panel C).
7.3 Conclusions
The case study in Section 7.1 illustrates the issue with parameter fitting in
population balance modelling: what if the surfacevolume model was used
in the place of the binarytree model for the fitting in Section 6.2? How
would the parameters have been different, and would the resultant model
have predicted a fundamentally different mechanism of particle formation?
These questions have a wider impact than the particular case of the sil
icon model analysed here. For example, the sintering parameters typically
used in many titania particle models were obtained from modelfitting with
the surfacevolume typespace [141, 165]. It was not until these empirical
expressions were augmented with insight from molecular dynamics (MD)
calculations that any reasonable confidence can be held in the sintering
parameters [19]. In essence, any fitted parameters should be used with
caution until they can be be evaluated with another modelling tool.
Further, the study in Section 7.1 highlights an example where experi
mental insight should be used with caution: if spherical particles are ob
tained experimentally, it does not necessarily justify the use of a spherical
particle model.
The second case study presented an approach for simulating the unsteady
state response of a hotwall reactor. The transient and axial evolution of the
quantities illustrates the power of the methodology developed in this thesis;
using a fullycoupled model incorporating gasphase kinetics and a multi
dimensional particle model enables an unprecedented amount of detail of
particle composition, size and structure to be captured.
121
Chapter 8
Conclusions
8.1 Findings of this work
The principal motivations of this thesis was to develop new methods and
models for simulating silicon nanoparticle synthesis. This was accomplished
using stochastic numerical methods to solve population balance models.
The methods developed include a new expression for Brownian coagula
tion, compatible with stochastic weighted algorithms; and an approach to
solve networks of linked reactors. Following formulation of a new multi
variate model for silicon nanoparticle synthesis, these tools and algorithms
were applied to the new silicon model in order to demonstrate how they
might be used.
In the first part of this thesis, three existing particle models were com
pared. The particle properties predicted by these models as a function of
coagulation rate and sintering rate were analysed, and different zones cor
responding to a particular particle structure were identified. It was found
that the models substantially differ from each other where coagulation and
sintering occur on similar timescales. This analysis illustrated the impor
tance of using multivariate particle models to accurately capture the size
and structure of aggregate particles.
Stochastic methods for solving the population balance equations have
only recently been extended to spatially inhomogeneous systems. Stochas
tic weighted algorithms (SWAs) were proposed in order to control the accu
mulation of statistical and systematic error in the solution of such systems
[119]. The transition regime coagulation kernel was adapted for use with
122
these methods, allowing for efficient calculation of the coagulation rate us
ing the technique of majorant rates and fictitious jumps. Recognising the
importance of multidimensional particle models, the new kernel was tested
using a detailed silica particle model. It was found that the new coagulation
algorithm, while more computationally expensive than conventional direct
simulation, exhibited the superior variancereduction expected of SWAs.
An algorithm was subsequently presented for solving the population bal
ance equations with the stochastic numerical method for networks of linked
reactors. The particle population balance was also fullycoupled with a gas
phase ODE model. The reactor network could represent a mesh of cells
in a computational fluid dynamics (CFD) simulation or a ‘compartmental
model’ approximating such a system. The convergence behaviour of the al
gorithm was tested as a function of network size, degree of feedback and
numerical parameters. It was demonstrated that use of SWAs was essential
in controlling the accumulation of error, particularly where a large quantity
of material was recycled.
A new model for silicon nanoparticle synthesis was then proposed. A
systematic parameter estimation procedure was used to determine the pre
exponential factors of five gasphase chemical reactions and two hetero
geneous surface reactions, fitting to a variety experimental studies over a
range of process conditions. The model was tested against six different
experimental configurations, with excellent agreement observed for the ma
jority of cases. This new model represents a major step forward in uniting
decades of research into silicon nanoparticle formation.
The new silicon model was analysed using the particle formation zones
identified in the first part of the thesis. The analysis highlighted that param
eters found through an estimation procedure are tied to the particle model
with which they were developed, unless they can be confirmed using other
sources of data. A reactor networking approach was employed to gain fur
ther insight into this silicon system, by studying the transient evolution of
a onedimensional tubular reactor. Here, the particle model was used to vi
sualise the changes in particle concentration, size and structure through the
reactor with time.
As stochastic methods are unrivalled in their ability to capture detail of
particle structure and morphology, it is anticipated that their use will play
a stronger role in process modelling. This thesis has contributed to the
123
understanding of practical and numerical aspects of such approaches, and
provided a new foundation from which the nature of silicon nanoparticle
synthesis may be explored.
8.2 Suggestions for future work
A more detailed particle model?
The binarytree particle model adapted to the silicon system in this
thesis is currently unrivalled in the amount of detail it can provide
about particle internal structure. However, it currently employs a
number of assumptions to simplify the calculation of particle prop
erties. For example, a constant fractal dimension is assumed to deter
mine the collision diameter. Alternatively, the fractal dimension could
be calculated dynamically from the structural information provided by
the particle model. This would yield a more complete description of
coagulation.
Influence of surface growth
In Chapter 3, the effect of coagulation and sintering on particle models
was assessed. However, it is welldocumented that surfacegrowth can
cause rounding of aggregate particles. It could therefore be useful to
append a fourth axis representing a characteristic surfacegrowth rate
to Figure 3.3, which could potentially elucidate the suitability of par
ticle models in accounting for surfacerounding. This would require
further investigation into the ‘surface smoothing factor’ (σSR) for the
binarytree particle model, as its current value is somewhat arbitrary
and it is not clear if this parameter can be linked to a physical property
of a particle system.
Compartmental models
The reactor networking tools presented in this thesis have many poten
tial applications. One of the most immediate areas for future investiga
tion is for the development of compartmental models. As discussed in
Chapter 5 and mentioned above, these models coarsely approximate
a CFD grid with linked reactors representing zones of flow. Since the
binarytree particle model is likely to be too computationally expensive
to directly couple to CFD calculations, a compartmental model could
124
be a useful middleground for applying a detailed particle model to a
system where fluid and particle transport is challenging to model.
Silicon colloids
The model for silicon nanoparticle synthesis was developed on the as
sumption that silane was the precursor gas. It has, however, been
demonstrated that disilane and higher silanes could be used to gen
erate silicon materials, such as the ‘silicon colloids’ of Fenollosa et al.
[42]. Confidence in the silicon model’s physical relevance could be
further improved if the model was evaluated for these higher silane
systems.
Silicon and silica models
It is wellknown that silicon will rapidly oxidise in contact with a
source of oxygen. Indeed, typical production of silicon nanoparticles
can yield particles with a thin oxide layer [75]. The silicon model pre
sented in this thesis could therefore be combined with the multivariate
silica particle model of Shekar et al. [144]. With an appropriate kinetic
expression to describe the conversion of silicon hydride to silicon diox
ide, the two models could be harmonised to provide a more realistic
view of industrial manufacture of particles.
125
Appendix A
Appendix to Chapter 3
A.1 Collision of particles in the nearspherical zone
Suppose two spherical particles with volume v1 and kv1 collide. How do the
derived properties of the resulting particles differ? The total surface are of
the resulting particle is given by:
S = pi
(
6v1
pi
) 2
3 [
1 + k
2
3
]
(A.1.1)
The equivalent spherical diameter is given by:
dsph =
(
6v1
pi
) 1
3 [
1 + k
] 1
3 (A.1.2)
Thus, the collision diameter as predicted by the surfacevolume model may
be calculated:
dsurfvolcol =
1
2
(
6v1
pi
) 1
3
[
(1 + k)
1
3 +
√
1 + k
2
3
]
(A.1.3)
The binarytree model requires estimation of the average primary diameter:
dpri =
1
2
(
6v1
pi
) 1
3 [
1 + k
1
3
]
(A.1.4)
Using this information, and calculating the ‘reduced number of primary par
ticles’ (given by the term inside the 1Df exponent) the binarytree model
calculates the collision diameter of the aggregate particle as:
dbintreecol =
1
2
(
6v1
pi
) 1
3 [
1 + k
1
3
]
[
1 + k
2
3
]3
[
1 + k
]2
1
Df
(A.1.5)
126
0.85
0.9
0.95
1
1.05
1.1
1.15
1.2
0.1 1 10
d c
o
l b
in
tr
ee
/ d
co
ls
u
rf
vo
l ,

k, 
Df=1.8Df=1.56
Figure A.1.1: Ratio of collision diameters predicted for the surfacevolume
and binary tree models as a function of relative colliding par
ticle size.
The two models are compared by plotting the ratio of dbintreecol to d
surfvol
col
using a fractal dimension Df = 1.8 and Df = 1.56 . It is evident that for like
sized particles (k ≈ 1) the collision diameter as predicted by the surface
volume model will be smaller than that predicted by the binarytree particle
model to a maximum of 10% underestimation.
A.2 Comparison of maximum sinterable surface area
Now suppose that a particle Pq consists of n equallysized primary particles
of volume v1. The equivalent spherical surface area is given by:
Ssph = pi
(
6nv1
pi
) 2
3
(A.2.1)
127
The initial surface area of such a particle as calculated by the surfacevolume
model is the sum of the surface area of each of the primaries:
Ssurfvol = npi
(
6v1
pi
) 2
3
(A.2.2)
In the binarytree model, there are n–1 connections between particles, each
of these with a surface area given by the sum of the surface area of the two
touching primaries. The total surface area considered by the binary tree
model is therefore given by:
Sbintree = 2(n – 1)pi
(
6v1
pi
) 2
3
(A.2.3)
If sintering occurs with the same characteristic time τS, the relative rates of
sintering can be compared between models by comparing the difference be
tween the initial and spherical surface area, termed the maximum sinterable
surface area:
∆S = S – Ssph (A.2.4)
To compare these between the binary tree and surfacevolume models, the
ratio of the maximum sinterable surface area (r) for between models is used:
r =
∆Sbintree
∆Ssurfvol
=
2(n – 1) – n2/3
n – n2/3
(A.2.5)
For n ≥ 2, r ≥ 1.0 and as such, more sinterable surface area exists in the
binary tree model, as compared to the surfacevolume model.
128
Appendix B
Appendix to Chapter 6
B.1 Conversion of SURFACE CHEMKIN parameters
SURFACE CHEMKIN was developed by Coltrin et al. [30] to model the com
plex surface chemistry in chemical vapour deposition processes. It consists
of a set of packages to interpret general heterogeneous chemical mecha
nisms in a standardised manner. As the gassurface mechanism proposed by
Houf et al. [64] and Ho et al. [63] was formulated in terms of a SURFACE
CHEMKIN mechanism, the form of the rate expression requires ‘conversion’
so that is compatible with the typespace of the present work’s model.
In a similar vein to Equation (2.1) An irreversible heterogeneous reac
tion may be represented in the general form [30]:
NGP∑
k=1
v′k χk →
NGP∑
k=1
v′′k χk (B.1.1)
Assuming only a single irreversible heterogeneous reaction is present, the
rate of progress variable q (in mol/cm2s) is given by:
q = kf
NGP∏
k=1
c
v′k
χk (B.1.2)
where cχk is the gasphase (mol/cm
3) or surface (mol/cm2) concentration
of species k. The forward rate constant kf adopts a standard Arrhenius form:
kf = A
* exp
(
–
EA
R T
)
(B.1.3)
129
Table B.1: Surface reaction mechanism for silanes proposed by Ho et al.
[63]. The (S) and (B) subscripts refer to surface and bulk
species, respectively.
A* EA
Reaction cm5/mol2s kcal/mol
SiH4 + 2Si(S) → 2SiH(S) + Si(B) + H2 8.39× 1026 37.45
Si2H6 + 2Si(S) → 2SiH(S) + 2Si(B) + 2H2 8.39× 1027 37.45
Si3H8 + 2Si(S) → 2SiH(S) + 3Si(B) + 3H2 8.39× 1027 37.45
To adapt the mechanisms for silane surface reactions [63, 64] to the
present work, Equations B.1.1–B.1.3 must be applied to the chemical reac
tions given in Table B.1.
As an example, this will be done for SiH4. Using Equation B.1.1, we
observe that the forward rate coefficients are v′SiH4 = 1 and v
′
Si(S)
= 2.
Thus, the rate of progress q for this reaction is given by:
q = A*c1SiH4 c
2
Si(S)
exp
(
–
EA
R T
)
(B.1.4)
Writing Equation B.1.4 in this form enables deduction of the units of A*.
Here, they must be cm5/mol2s in order to ensure dimensional consistency.
Assuming that the concentration of surface Si(S) sites–given by λSi (#/cm
2)–
is constant, we write:
λSi = NA cSi(S) (B.1.5)
To ‘translate’ Equation B.1.4 to the typespace of the present work, we
must ensure that the surface reaction rate RSR is dimensionally consistent
with the other rates:
RSR = S NA q (B.1.6)
Substituting in for λSi and rearranging yields:
RSR = A cSiH4 S
(
–
EA
R T
)
(B.1.7)
where A is the alternative preexponential factor used in the present work,
given by:
A =
A* λ2Si
NA
(B.1.8)
130
Table B.2: Estimation of ASR,SiH4 based on different site densities.
Parameter Lower bound This work Upper bound
λSi (#/cm
2) 8.0× 1014 [140]  1.7× 1015
ASR,SiH4 (cm/mol.s) 8.9× 1032 3.0× 1033 4.0× 1033
By estimating the site density λSi, it is possible to convert the value of the
preexponential factor A* into the A used in present work. An upperbound
on λSi can be estimated by using the bulk density of silicon. A lower bound
on λSi is taken for the {111} surface of silicon at full saturation [140].
We observe that Equation B.1.7 is the same as that used by Menz et al.
[92]; and that the parameters for silane in Table B.2 shows good agreement
with their preexponential factors calculated via optimisation (ASR,SiH4 =
3.0× 1033) in the present work.
131
Nomenclature
Uppercase Roman
A Arrhenius preexponential factor
C Connectivity matrix
Da Damköhler number
Df Fractal dimension
EA Activation energy
Fc Ratio of downstream to upstream scaling factors
FV Volume fraction of particles
Fs Outflow rescaling factor
G Gasphase statespace
Hˆ Specific molar enthalpy
K Coagulation kernel
Kc Chemical equilibrium constant
Kfm Freemolecular coagulation kernel
Ksf Slipflow coagulation kernel
Ktr Transition regime coagulation kernel
Kn Knudsen number
L Number of independent runs
M Thirdbody collision species symbol
Mi(t) Moment i of the particle ensemble
NA Avogadro’s number
NGP Number of gasphase chemical species
Nh Number of factorised terms in majorant kernel
Nrxn Number of gasphase chemical reactions
132
Nmax Maximum number of computational particles
N(t) Number of computational particles at time t
P Computational particle type
P Pressure
Pr Reduced pressure
Q Volumetric flowrate
R Ideal gas constant
Ri Rate of particle process i
S Surface area
SSi Supersaturation of silicon
Ssph Equivalent spherical surface area
T Temperature (in Kelvin)
V Volume
Vsmp Sample volume for stochastic particle method
X(t) Functional describing some particle ensemble property
Z Particle statespace
Lowercase Roman
c Molar concentration of chemical species
dcol Particle collision diameter
dmob Particle mobility diameter
dpri Primary particle diameter
f Fractional inflow
g˙ Molar production rate due to particle processes
h(j)i Upper bound for majorant kernel
m Mass of particle
n Number of primary particles
n(t, P) Number concentration of particles
k Chemical rate coefficient
k0 Lowpressure limit rate coefficient
k∞ Highpressure limit rate coefficient
kB Boltzmann constant
133
kmaj Majorant constant
t Time
p Primary particle
pSi Partial pressure of silicon
pSi,sat Saturation vapour pressure of silicon
q Chemical rate of progress variable
s Sintering level
v Primary particle volume
w Statistical weight
x Parameter vector
y Mole fraction
Uppercase Greek
Φ Objective function
Lowercase Greek
α Freemolecular coagulation parameter
β Arrhenius temperature exponent
β1,β2 Slipflow coagulation parameters
γ(Pi, wi, Pj, wj) Weight transfer function [122]
γsi Surface energy of silicon
δv Change in volume due to surface reaction
Relative numerical error
stat Relative numerical error due to statistical fluctuations
ηi Number of chemical units of species i
θ Molar coverage of species on surface
λ Partial sum in particle binary tree
µ(X)(t) Average quantity of functional X(t) over number of particles
ν Stochiometric coefficient
ξ Average quantity (over runs)
ρ Density (mass or molar)
134
σ Mean free path
σSR Surface smoothing factor
ς(P) Function representing sintering
τ Residence or characteristic time
φ Data point for fitting in objective function
χ Chemical symbol
ψ Gasphase expansion rate
ω˙ Molar production rate due to gasphase reactions
Superscripts
(q) Particle property index
[j] Stream index
Subscripts
a Arithmetic
g Geometric
S Denotes sintering property
SR Denotes property of a surface reaction
Abbreviations
CFD Computational fluid dynamics
CSTR Continuousstirred tank reactor
DSA Direct simulation algorithm
MD Molecular dynamics
ODE Ordinary differential equation
PFR Plugflow reactor
PSD Particle size distribution
PSR Perfectlystirred reactor
SWA Stochastic weighted algorithm
135
Bibliography
[1] A. J. Adamczyk, M. F. Reyniers, G. B. Marin, and L. J. Broadbelt. Ex
ploring 1, 2hydrogen shift in silicon nanoparticles: Reaction kinet
ics from quantum chemical calculations and derivation of transition
state group additivity database. The Journal of Physical Chemistry A,
113(41):10933–10946, 2009.
[2] M. K. Akhtar, Y. Xiong, and S. E. Pratsinis. Vapor synthesis of titania
powder by titanium tetrachloride oxidation. AIChE Journal, 37(10):
1561–1570, 1991.
[3] J. Akroyd, A. J. Smith, L. R. McGlashan, and M. Kraft. Numerical in
vestigation of DQMoMIEM as a turbulent reaction closure. Chemical
Engineering Science, 65(6):1915–1924, 2010.
[4] J. Akroyd, A. J. Smith, R. Shirley, L. R. McGlashan, and M. Kraft. A
coupled CFDpopulation balance approach for nanoparticle synthesis
in turbulent reacting flows. Chemical Engineering Science, 66(17):
3792–3805, 2011.
[5] A. H. Alexopoulos, P. Pladis, and C. Kiparissides. Nonhomogeneous
mixing population balance model for the prediction of particle size
distribution in large scale emulsion polymerization reactors. Indus
trial & Engineering Chemistry Research, 52:12285–12296, 2013.
[6] ANSYS, Inc. ANSYS FLUENT 12. 0 User’s Guide, 2009.
[7] H. Antoniadis. Silicon ink high efficiency solar cells. In Photovoltaic
Specialists Conference (PVSC), 2009 34th IEEE, pages 650 –654, 2009.
[8] M. Balthasar and M. Kraft. A stochastic approach to calculate the par
ticle size distribution function of soot particles in laminar premixed
flames. Combustion and Flame, 133(3):289–298, 2003.
136
[9] J. F. Banfield and H. Zhang. Nanoparticles in the environment. Re
views in Mineralogy and Geochemistry, 44(1):1–58, 2001.
[10] S. K. Bermingham, P. J. T. Verheijen, and H. J. M. Kramer. Opti
mal design of solution crystallization processes with rigorous models.
Chemical Engineering Research and Design, 81(8):893–903, 2003.
[11] F. Bezzo and S. Macchietto. A general methodology for hybrid multi
zonal/CFD models: Part II. Automatic zoning. Computers & Chemical
Engineering, 28(4):513–525, 2004.
[12] F. Bezzo, S. Macchietto, and C. C. Pantelides. A general methodology
for hybrid multizonal/CFD models: Part I. Theoretical framework.
Computers & Chemical Engineering, 28(4):501–511, 2004.
[13] A. Bhave and M. Kraft. Partially stirred reactor model: Analytical
solutions and numerical convergence study of a PDF/Monte Carlo
method. SIAM Journal on Scientific Computing, 25(5):1798–1823,
2004.
[14] G. Biskos, V. Vons, C. U. Yurteri, and A. SchmidtOtt. Generation and
sizing of particles for aerosolbased nanotechnology. KONA Powder
Particle J, 26:13–35, 2008.
[15] G. Blanquart and H. Pitsch. Analyzing the effects of temperature on
soot formation with a joint volumesurfacehydrogen model. Com
bustion and Flame, 156(8):1614–1626, 2009.
[16] M. Bowrey and J. H. Purnell. The pyrolysis of disilane and rate con
stants of silene insertion reactions. Proceedings of the Royal Society of
London A, 321(1546):341–359, 1971.
[17] A. Braumann, P. L. W. Man, and M. Kraft. The inverse problem in
granulation modeling–Two different statistical approaches. AIChE
Journal, 57(11):3105–3121, 2011.
[18] V. S. Buddhiraju and V. Runkana. Simulation of nanoparticle syn
thesis in an aerosol flame reactor using a coupled flame dynamics–
monodisperse population balance model. Journal of Aerosol Science,
43(1):1–13, 2012.
137
[19] B. Buesser, A. J. Gröhn, and SE Pratsinis. Sintering rate and mech
anism of TiO2 nanoparticles by molecular dynamics. The Journal of
Physical Chemistry C, 115(22):11030–11035, 2011.
[20] A. Burcat. Thermochemical data for combustion calculations. In Com
bustion Chemistry, pages 455–473. Springer, 1984.
[21] R. J. Buss, P. Ho, W. G. Breiland, and M. E. Coltrin. Reactive sticking
coefficients for silane and disilane on polycrystalline silicon. Journal
of Applied Physics, 63(8):2808–2819, 1988.
[22] W. R. Cannon, S. C. Danforth, J. H. Flint, J. S. Haggerty, and R. A.
Marra. Sinterable ceramic powders from laserdriven reactions: I,
Process description and modeling. Journal of the American Ceramic
Society, 65(7):324–330, 1982.
[23] W. R. Cannon, S. C. Danforth, J. S. Haggerty, and R. A. Marra. Sinter
able ceramic powders from laserdriven reactions: II, Powder char
acteristics and process variables. Journal of the American Ceramic
Society, 65(7):330–335, 1982.
[24] Cantera Developers. Cantera: An objectoriented software toolkit for
chemical kinetics, thermodynamics, and transport processes, 2013.
URL http://code.google.com/p/cantera.
[25] M. Celnik, R.I.A. Patterson, M. Kraft, and W. Wagner. Coupling a
stochastic soot population balance to gasphase chemistry using op
erator splitting. Combustion and Flame, 148(3):158–176, 2007.
[26] M. Celnik, R.I.A. Patterson, M. Kraft, and W. Wagner. A predictor
corrector algorithm for the coupling of stiff ODEs to a particle popu
lation balance. Journal of Computational Physics, 228(8):2758–2769,
2009.
[27] D. Chen, Z. Zainuddin, E. Yapp, J. Akroyd, S. Mosbach, and M. Kraft.
A fully coupled simulation of PAH and soot growth with a population
balance model. Proceedings of the Combustion Institute, 34:1827–
1835, 2013.
[28] Z. Chen, H. Meng, G. Xing, H. Yuan, F. Zhao, R. Liu, X. Chang, X. Gao,
T. Wang, and G. Jia. Agerelated differences in pulmonary and car
138
diovascular responses to SiO2 nanoparticle inhalation: nanotoxicity
has susceptible population. Environmental Science & Technology, 42
(23):8985–8992, 2008.
[29] K.Y. Cheng, R. Anthony, U. R. Kortshagen, and R. J. Holmes. Hy
brid silicon nanocrystal organic lightemitting devices for infrared
electroluminescence. Nano Letters, 10(4):1154–1157, 2010.
[30] M. E. Coltrin, R. J. Kee, F. M. Rupley, and E. Meeks. SURFACE
CHEMKIN III: A Fortran Package for Analyzing Heterogeneous Chem
ical Kinetics at a SolidSurfaceGasPhase Interface. Sandia National
Laboratories, 1996.
[31] H. Dang and M. T. Swihart. Computational modeling of silicon
nanoparticle synthesis: I. A general twodimensional model. Aerosol
Science and Technology, 43(3):250–263, 2009.
[32] H. Dang and M. T. Swihart. Computational modeling of silicon
nanoparticle synthesis: II. A twodimensional bivariate model for sil
icon nanoparticle synthesis in a laserdriven reactor including finite
rate coalescence. Aerosol Science and Technology, 43(6):554–569,
2009.
[33] E. Debry, B. Sportisse, and B. Jourdain. A stochastic approach for the
numerical simulation of the general dynamics equation for aerosols.
Journal of Computational Physics, 184(2):649–669, 2003.
[34] R. E. L. DeVille, N. Riemer, and M. West. Weighted Flow Algorithms
(WFA) for stochastic particle coagulation. Journal of Computational
Physics, 230:8427–8451, 2011.
[35] K. E. Drexler. Engines of Creation. Fourth Estate, 1986.
[36] M. L. Eggersdorfer, D. Kadau, H. J. Herrmann, and S. E. Pratsinis.
Aggregate morphology evolution by sintering: Number & diameter
of primary particles. Journal of Aerosol Science, 46:7–19, 2012.
[37] A. Eibeck and W. Wagner. An efficient stochastic algorithm for study
ing coagulation dynamics and gelation phenomena. SIAM Journal on
Scientific Computing, 22(3):802–821, 2001.
139
[38] A. Eibeck and W. Wagner. Stochastic particle approximations for
Smoluchoski’s coagulation equation. Annals of Applied Probability,
pages 1137–1165, 2001.
[39] D. Erhardt. Materials conservation: Notsonew technology. Nature
Materials, 2(8):509–510, 2003.
[40] F. Erogbogbo, K.T. Yong, I. Roy, G. Xu, P. N. Prasad, and M. T. Swi
hart. Biocompatible luminescent silicon quantum dots for imaging of
cancer cells. ACS Nano, 2(5):873–878, 2008.
[41] F. C. Eversteijn. Gasphase decomposition of silane in a horizontal
epitaxial reactor. Philips Res. Repts, 26(2):134–144, 1971.
[42] R. Fenollosa, F. RamiroManzano, M. Tymczenko, and F. Meseguer.
Porous silicon microspheres: synthesis, characterization and appli
cation to photonic microcavities. Journal of Materials Chemistry, 20
(25):5210–5214, 2010.
[43] T. Finsterbusch and A. Mankertz. Porcine circoviruses–small but pow
erful. Virus Research, 143(2):177–183, 2009.
[44] J. H. Flint and J. S. Haggerty. A model for the growth of silicon
particles from laserheated gases. Aerosol Science and Technology, 13
(1):72–84, 1990.
[45] J. H. Flint, R. A. Marra, and J. S. Haggerty. Powder temperature,
size, and number density in laserdriven reactions. Aerosol Science
and Technology, 5(2):249–260, 1986.
[46] H. S. Fogler. Elements of Chemical Reaction Engineering. Pearson, New
Jersey, 2006.
[47] M. Frenklach. Method of moments with interpolative closure. Chem
ical Engineering Science, 57(12):2229–2239, 2002.
[48] M. Frenklach and S. J. Harris. Aerosol dynamics modeling using the
method of moments. Journal of Colloid and Interface Science, 118(1):
252–261, 1987.
[49] M. Frenklach, L. Ting, H. Wang, and M. J. Rabinowitz. Silicon par
ticle formation in pyrolysis of silane and disilane. Israel Journal of
Chemistry, 36(3):293–303, 1996.
140
[50] A. Gallagher, A. A. Howling, and C. H. Hollenstein. Anion reactions
in silane plasma. Journal of Applied Physics, 91(9):5571–5580, 2002.
[51] B. Giesen, H. Wiggers, A. Kowalik, and P. Roth. Formation of Si
nanoparticles in a microwave reactor: Comparison between experi
ments and modelling. Journal of Nanoparticle Research, 7(1):29–41,
2005.
[52] C. J. Giunta, R. J. McCurdy, J. D. ChappleSokol, and R. G. Gordon.
Gasphase kinetics in the atmospheric pressure chemical vapor depo
sition of silicon from silane and disilane. Journal of Applied Physics,
67(2):1062–1075, 1990.
[53] M. Goodson and M. Kraft. An efficient stochastic algorithm for simu
lating nanoparticle dynamics. Journal of Computational Physics, 183
(1):210–232, 2002.
[54] M. Goodson and M. Kraft. Simulation of coalescence and break
age: an assessment of two stochastic methods suitable for simulating
liquidliquid extraction. Chemical Engineering Science, 59(18):3865–
3881, 2004.
[55] D. E. Gray, editor. American Institute of Physics Handbook. McGraw
Hill, New York, 1972.
[56] M. Gröschel, R. Körmer, M. Walther, G. Leugering, and W. Peuk
ert. Process control strategies for the gas phase synthesis of silicon
nanoparticles. Chemical Engineering Science, 73:181–194, 2012.
[57] F. Guias¸. A stochastic numerical method for diffusion equations and
applications to spatially inhomogeneous coagulation processes. In
Monte Carlo and QuasiMonte Carlo Methods 2004, pages 147–161.
Springer, 2006.
[58] A. Gupta and H. Wiggers. Surface chemistry and photoluminescence
property of functionalized silicon nanoparticles. Physica E, 41(6):
1010–1014, 2009.
[59] Z. Haibo, Z. Chuguang, and X. Minghou. MultiMonte Carlo ap
proach for general dynamic equation considering simultaneous parti
141
cle coagulation and breakage. Powder Technology, 154(23):164–178,
2005.
[60] S. J. Harris and I. M. Kennedy. The coagulation of soot particles with
van der waals forces. Combustion Science and Technology, 59(46):
443–454, 1988.
[61] M. C. Heine and S. E. Pratsinis. Polydispersity of primary particles in
agglomerates made by coagulation and sintering. Journal of Aerosol
Science, 38(1):17–38, 2007.
[62] T. Hirokami, Y. Maeda, and H. Tsukada. Parameter estimation using
simultaneous perturbation stochastic approximation. Electrical Engi
neering in Japan, 154:30–39, 2006.
[63] P. Ho, M. E. Coltrin, and W. G. Breiland. Laserinduced fluorescence
measurements and kinetic analysis of Si atom formation in a rotat
ing disk chemical vapor deposition reactor. The Journal of Physical
Chemistry, 98(40):10138–10147, 1994.
[64] W. G. Houf, J. F. Grcar, and W. G. Breiland. A model for low pressure
chemical vapor deposition in a hotwall tubular reactor. Materials
Science and Engineering: B, 17(13):163–171, 1993.
[65] M. J. Hounslow. A discretized population balance for continuous
systems at steady state. AIChE Journal, 36(1):106–116, 1990.
[66] A. A. Howling, J. L Dorier, C. H. Hollenstein, U. Kroll, and F. Finger.
Frequency effects in silane plasmas for plasma enhanced chemical
vapor deposition. Journal of Vacuum Science & Technology A, 10(4):
1080–1085, 1992.
[67] R. Irizarry. Fast compartmental monte carlo simulation of population
balance models: Application to nanoparticle formation in nonhomo
geneous conditions. Industrial & Engineering Chemistry Research, 51
(47):15484–15496, 2012.
[68] T. Johannessen, S. E. Pratsinis, and H. Livbjerg. Computational fluid
particle dynamics for the flame synthesis of alumina particles. Chem
ical Engineering Science, 55(1):177–191, 2000.
142
[69] K. C. Kannenberg and I. D. Boyd. Strategies for efficient particle
resolution in the direct simulation monte carlo method. Journal of
Computational Physics, 157(2):727–745, 2000.
[70] A. Kazakov and M. Frenklach. Dynamic modeling of soot particle
coagulation and aggregation: Implementation with the method of
moments and application to highpressure laminar premixed flames.
Combustion and Flame, 114(3):484–501, 1998.
[71] Robert J Kee, Fran M Rupley, Ellen Meeks, and James A Miller.
CHEMKINIII: A FORTRAN chemical kinetics package for the analysis
of gasphase chemical and plasma kinetics. Sandia National Laborato
ries, 1996.
[72] J. Knipping, H. Wiggers, B. Rellinghaus, P. Roth, D. Konjhodzic, and
C. Meier. Synthesis of high purity silicon nanoparticles in a low pres
sure microwave reactor. Journal of Nanoscience and Nanotechnology,
4(8):1039–1044, 2004.
[73] W. Koch and SK Friedlander. The effect of particle coalescence on the
surface area of a coagulating aerosol. Journal of Colloid and Interface
Science, 140(2):419–427, 1990.
[74] A. Kolodko and K. Sabelfeld. Stochastic particle methods for Smolu
chowski coagulation equation: variance reduction and error estima
tions. Monte Carlo Methods Applied, 9(4):315–339, 2003.
[75] R. Körmer, M. P. M. Jank, H. Ryssel, H. J. Schmid, and W. Peuk
ert. Aerosol synthesis of silicon nanoparticles with narrow size
distribution–Part 1: Experimental investigations. Journal of Aerosol
Science, 41(11):998–1007, 2010.
[76] R. Körmer, H. J. Schmid, and W. Peukert. Aerosol synthesis of silicon
nanoparticles with narrow size distribution–Part 2: Theoretical anal
ysis of the formation mechanism. Journal of Aerosol Science, 41(11):
1008–1019, 2010.
[77] M. Koshi, S. Kato, and H. Matsui. Unimolecular decomposition of
silane, fluorosilane, and difluorosilane at high temperatures. The
Journal of Physical Chemistry, 95(3):1223–1227, 1991.
143
[78] F. E. Kruis, K. Kusters, S. E. Pratsinis, and B. Scarlett. A simple model
for the evolution of the characteristics of aggregate particles under
going coagulation and sintering. Aerosol Science and Technology, 19:
514–526, 1993.
[79] F. E. Kruis, J. Schoonman, and B. Scarlett. Homogeneous nucleation
of silicon. Journal of Aerosol Science, 25(7):1291–1304, 1994.
[80] F. E. Kruis, H. Fissan, and A. Peled. Synthesis of nanoparticles in the
gas phase for electronic, optical and magnetic applications–a review.
Journal of Aerosol Science, 29(5):511–535, 1998.
[81] F. E. Kruis, A. Maisels, and H. Fissan. Direct simulation Monte Carlo
method for particle coagulation and aggregation. AIChE Journal, 46
(9):1735–1742, 2000.
[82] F. E. Kruis, J. Wei, T. van der Zwaag, and S. Haep. Computational
fluid dynamics based stochastic aerosol modeling: Combination of
a cellbased weighted random walk method and a constantnumber
MonteCarlo method for aerosol dynamics. Chemical Engineering Sci
ence, 70:109–120, 2012.
[83] P. Lavvas, M. Sander, M. Kraft, and H. Imanaka. Surface chemistry
and particle shape: Processes for the evolution of aerosols in titan’s
atmosphere. The Astrophysical Journal, 728(2):80, 2011.
[84] J. Li, B. Freireich, C. Wassgren, and J. D. Litster. A general
compartmentbased population balance model for particle coating
and layered granulation. AIChE Journal, 58(5):1397–1408, 2012.
[85] X. Li, Y. He, and M. T. Swihart. Surface functionalization of silicon
nanoparticles produced by laserdriven pyrolysis of silane followed
by HFHNO3 etching. Langmuir, 20(11):4720–4727, 2004.
[86] Z. F. Li and E. Ruckenstein. Watersoluble poly (acrylic acid) grafted
luminescent silicon nanoparticles and their use as fluorescent biolog
ical staining labels. Nano Letters, 4(8):1463–1467, 2004.
[87] F. A. Lindemann, S. Arrhenius, I. Langmuir, N. R. Dhar, J. Perrin, and
W. C. McC. Lewis. Discussion on “the radiation theory of chemical
action”. Transactions of the Faraday Society, 17:598–606, 1922.
144
[88] L. Mangolini. Synthesis, properties, and applications of silicon
nanocrystals. Journal of Vacuum Science & Technology B, 31(2):
020801, 2013.
[89] D. L. Marchisio and R. O. Fox. Solution of population balance equa
tions using the direct quadrature method of moments. Journal of
Aerosol Science, 36(1):43–73, 2005.
[90] W. J. Menz and M. Kraft. A new model for silicon nanoparticle syn
thesis. Combustion and Flame, 160:947–958, 2013.
[91] W. J. Menz and M. Kraft. The suitability of particle models in captur
ing aggregate structure and polydispersity. Aerosol Science & Technol
ogy, 47:734–745, 2013.
[92] W. J. Menz, S. Shekar, G. Brownbridge, S. Mosbach, R. Körmer,
W. Peukert, and M. Kraft. Synthesis of silicon nanoparticles with
a narrow size distribution: A theoretical study. Journal of Aerosol
Science, 44:46–61, 2012.
[93] W. J. Menz, R. I. A. Patterson, W. Wagner, and M. Kraft. Application of
stochastic weighted algorithms to a multidimensional silica particle
model. Journal of Computational Physics, 248:221–234, 2013.
[94] W. J. Menz, J. Akroyd, and M. Kraft. Stochastic solution of population
balance equations for reactor networks. Journal of Computational
Physics, 256:615–629, 2014.
[95] M. Meunier, J. H. Flint, J. S. Haggerty, and D. Adler. Laserinduced
chemical vapor deposition of hydrogenated amorphous silicon. I.
Gasphase process model. Journal of Applied Physics, 62(7):2812–
2821, 1987.
[96] M. Meunier, J. H. Flint, J. S. Haggerty, and D. Adler. Laserinduced
chemical vapor deposition of hydrogenated amorphous silicon. II.
Film properties. Journal of Applied Physics, 62(7):2821–2829, 1987.
[97] X. Michalet, F. F. Pinaud, L. A. Bentolila, J. M. Tsay, S. Doose, J. J. Li,
G. Sundaresan, A. M. Wu, S. S. Gambhir, and S. Weiss. Quantum dots
for live cells, in vivo imaging, and diagnostics. Science, 307(5709):
538–544, 2005.
145
[98] N. M. Morgan, C. G. Wells, M. Kraft, and W. Wagner. Modelling
nanoparticle dynamics: coagulation, sintering, particle inception and
surface growth. Combustion Theory and Modelling, 9(3):449–461,
2005.
[99] N. M. Morgan, C. G. Wells, M. J. Goodson, M. Kraft, and W. Wag
ner. A new numerical approach for the simulation of the growth of
inorganic nanoparticles. Journal of Computational Physics, 211(2):
638–658, 2006.
[100] S. Mosbach, A. Braumann, P. L. W. Man, C. A. Kastner, G. P. E. Brown
bridge, and M. Kraft. Iterative improvement of bayesian parame
ter estimates for an engine model by means of experimental design.
Combustion and Flame, 159(3):1303–1313, 2012.
[101] T. Murthy, N. Miyamoto, M. Shimbo, and J. Nishizawa. Gasphase
nucleation during the thermal decomposition of silane in hydrogen.
Journal of Crystal Growth, 33(1):1–7, 1976.
[102] O. M. Nayfeh, D. A. Antoniadis, K. Mantey, and M. H. Nayfeh. Mem
ory effects in metaloxidesemiconductor capacitors incorporating
dispensed highly monodisperse 1 nm silicon nanoparticles. Applied
Physics Letters, 90(15):153105–153105, 2007.
[103] H. V. Nguyen and R. C. Flagan. Particle formation and growth in
singlestage aerosol reactors. Langmuir, 7(8):1807–1814, 1991.
[104] S. Nijhawan, P. H. McMurry, M. T. Swihart, S. M. Suh, S. L. Girshick,
S. A. Campbell, and J. E. Brockmann. An experimental and numerical
study of particle nucleation and growth during lowpressure thermal
decomposition of silane. Journal of Aerosol Science, 34(6):691–711,
2003.
[105] B. Nowack and T. D. Bucheli. Occurrence, behavior and effects of
nanoparticles in the environment. Environmental Pollution, 150(1):
5–22, 2007.
[106] T. R. Nurkiewicz, D. W. Porter, A. F. Hubbs, J. L. Cumpston, B. T.
Chen, D. G. Frazer, and V. Castranova. Nanoparticle inhalation aug
ments particledependent systemic microvascular dysfunction. Parti
cle and Fibre Toxicology, 5(1), 2008.
146
[107] J. O. Odden, P. K. Egeberg, and A. Kjekshus. From monosilane to
crystalline silicon, Part I: Decomposition of monosilane at 690–830
K and initial pressures 0. 1–6. 6 MPa in a freespace reactor. Solar
Energy Materials and Solar Cells, 86(2):165–176, 2005.
[108] J. O. Odden, P. K. Egeberg, and A. Kjekshus. From monosilane to
crystalline silicon, Part III. Characterization of amorphous, hydrogen
containing silicon products. Journal of NonCrystalline Solids, 351
(14):1317–1327, 2005.
[109] J. O. Odden, P. K. Egeberg, and A. Kjekshus. From monosilane to
crystalline silicon, Part II: Kinetic considerations on thermal decom
position of pressurized monosilane. International Journal of Chemical
Kinetics, 38(5):309–321, 2006.
[110] A. A. Onischuk, V. P. Strunin, R. I. Samoilova, A. V. Nosov, M. A.
Ushakova, and V. N. Panfilov. Chemical composition and bond struc
ture of aerosol particles of amorphous hydrogenated silicon forming
from thermal decomposition of silane. Journal of Aerosol Science, 28
(8):1425–1441, 1997.
[111] A. A. Onischuk, V. P. Strunin, M. A. Ushakova, and V. N. Panfilov.
On the pathways of aerosol formation by thermal decomposition of
silane. Journal of Aerosol Science, 28(2):207–222, 1997.
[112] A. A. Onischuk, V. P. Strunin, M. A. Ushakova, and V. N. Panfilov.
Studying of silane thermal decomposition mechanism. International
Journal of Chemical Kinetics, 30(2):99–110, 1998.
[113] A. A. Onischuk, A. I. Levykin, V. P. Strunin, K. K. Sabelfeld, and V. N.
Panfilov. Aggregate formation under homogeneous silane thermal
decomposition. Journal of Aerosol Science, 31(11):1263–1281, 2000.
[114] A. A. Onischuk, A. I. Levykin, V. P. Strunin, M. A. Ushakova, R.. I
Samoilova, K. K. Sabelfeld, and V. N. Panfilov. Aerosol formation un
der heterogeneous/homogeneous thermal decomposition of silane:
experiment and numerical modeling. Journal of Aerosol Science, 31
(8):879–906, 2000.
147
[115] J.H. Park, L. Gu, G. von Maltzahn, E. Ruoslahti, S. N. Bhatia, and
M. J. Sailor. Biodegradable luminescent porous silicon nanoparticles
for in vivo applications. Nature Materials, 8(4):331–336, 2009.
[116] S. H. Park and S. N. Rogak. A novel fixedsectional model for the
formation and growth of aerosol agglomerates. Journal of Aerosol
Science, 35(11):1385–1404, 2004.
[117] R. I. A. Patterson and M. Kraft. Models for the aggregate structure of
soot particles. Combustion and Flame, 151(12):160–172, 2007.
[118] R. I. A. Patterson and W. Wagner. A stochastic weighted particle
method for coagulation–advection problems. SIAM Journal on Sci
entific Computing, 34(3):B290–B311, 2012.
[119] R. I. A. Patterson and W. Wagner. A stochastic weighted particle
method for coagulation–advection problems. SIAM Journal on Sci
entific Computing, 34(3):290–311, 2012.
[120] R. I. A. Patterson, J. Singh, M. Balthasar, M. Kraft, and J. R. Norris.
The linear process deferment algorithm: A new technique for solving
population balance equations. SIAM Journal on Scientific Computing,
28(1):303, 2006.
[121] R. I. A. Patterson, J. Singh, M. Balthasar, M. Kraft, and W. Wagner.
Extending stochastic soot simulation to higher pressures. Combustion
and Flame, 145(3):638–642, 2006.
[122] R. I. A. Patterson, W. Wagner, and M. Kraft. Stochastic weighted
particle methods for population balance equations. Journal of Com
putational Physics, 230:7456–7472, 2011.
[123] H. R. Paur, W. Baumann, H. Mätzing, and H. Seifert. Formation of
nanoparticles in flames; measurement by particle mass spectrometry
and numerical simulation. Nanotechnology, 16:S354, 2005.
[124] Persistence of Vision Pty. Ltd. Persistence of Vision Raytracer (Version
3. 6), 2004. URL http://www.povray.org/.
[125] E. L. Petersen and M. W. Crofton. Measurements of hightemperature
silane pyrolysis using SiH4 IR emission and SiH2 laser absorption.
The Journal of Physical Chemistry A, 107(50):10988–10995, 2003.
148
[126] W. Phadungsukanan, S. Shekar, R. Shirley, M. Sander, R. H. West, and
M. Kraft. Firstprinciples thermochemistry for silicon species in the
decomposition of tetraethoxysilane. The Journal of Physical Chemistry
A, 113(31):9041–9049, 2009.
[127] S. E. Pratsinis. Simultaneous nucleation, condensation, and coagula
tion in aerosol reactors. Journal of Colloid and Interface Science, 124
(2):416–427, 1988.
[128] S. E. Pratsinis and P. T. Spicer. Competition between gas phase and
surface oxidation of TiCl4 during synthesis of TiO2 particles. Chemi
cal Engineering Science, 53(10):1861–1868, 1998.
[129] J. H. Purnell and R. Walsh. The pyrolysis of monosilane. Proceedings
of the Royal Society of London A, 293(1435):543–561, 1966.
[130] J. Pyykönen and J. Jokiniemi. Computational fluid dynamics based
sectional aerosol modelling schemes. Journal of Aerosol Science, 31
(5):531–550, 2000.
[131] Reaction Design. CHEMKIN Software Theory Manual, 2011.
[132] Reaction Design. CHEMKIN, retrieved 8 May 2013. URL http://
www.reactiondesign.com/chemkin/.
[133] N. Riemer, M. West, R. A. Zaveri, and R. C. Easter. Simulating the
evolution of soot mixing state with a particleresolved aerosol model.
Journal of Geophysical Research, 114:D09202, 2009.
[134] S. Rjasanow and W. Wagner. A stochastic weighted particle method
for the Boltzmann equation. Journal of Computational Physics, 124
(2):243–253, 1996.
[135] S. N. Rogak, R. C. Flagan, and H. V. Nguyen. The mobility and struc
ture of aerosol agglomerates. Aerosol Science and Technology, 18(1):
25–47, 1993.
[136] M. Sander, R. H. West, M. S. Celnik, and M. Kraft. A detailed model
for the sintering of polydispersed nanoparticle agglomerates. Aerosol
Science and Technology, 43(10):978–989, 2009.
149
[137] M. Sander, R. I. A. Patterson, A. Braumann, A. Raj, and M. Kraft.
Developing the PAHPP soot particle model using process informatics
and uncertainty propagation. Proceedings of the Combustion Institute,
33(1):675–683, 2011.
[138] D. W. Schaefer and A. J. Hurd. Growth and structure of combustion
aerosols: fumed silica. Aerosol Science and Technology, 12(4):876–
890, 1990.
[139] H. J. Schmid, S. Tejwani, C. Artelt, and W. Peukert. Monte carlo sim
ulation of aggregate morphology for simultaneous coagulation and
sintering. Journal of Nanoparticle Research, 6(6):613–626, 2004.
[140] G. Schulze and M. Henzler. Adsorption of atomic hydrogen on clean
cleaved silicon (111). Surface Science, 124(2):336–350, 1983.
[141] T. Seto, M. Shimada, and K. Okuyama. Evaluation of sintering of
nanometersized titania using aerosol method. Aerosol Science and
Technology, 23(2):183–200, 1995.
[142] S. Shekar, M. Sander, R. C. Riehl, A. J. Smith, A. Braumann, and
M. Kraft. Modelling the flame synthesis of silica nanoparticles from
tetraethoxysilane. Chemical Engineering Science, 70:54–66, 2011.
[143] S. Shekar, W. J. Menz, A. J. Smith, M. Kraft, and W. Wagner. On a
multivariate population balance model to describe the structure and
composition of silica nanoparticles. Computers & Chemical Engineer
ing, 43:130–147, 2012.
[144] S. Shekar, A. J. Smith, W. J. Menz, M. Sander, and M. Kraft. A mul
tidimensional population balance model to describe the aerosol syn
thesis of silica nanoparticles. Journal of Aerosol Science, 44:83–98,
2012.
[145] R. Shirley, W. Phadungsukanan, M. Kraft, J. Downing, N. E. Day,
and P. MurrayRust. Firstprinciples thermochemistry for gas phase
species in an industrial rutile chlorinator. The Journal of Physical
Chemistry A, 114(43):11825–11832, 2010.
[146] B. W. Silverman. Density estimation for statistics and data analysis.
Chapman & Hall/CRC, 1986.
150
[147] K. Sinniah, M. G. Sherman, L. B. Lewis, W. H. Weinberg, J. T. Yates Jr,
and K. C. Janda. Hydrogen desorption from the monohydride phase
on Si (100). The Journal of Chemical Physics, 92:5700, 1990.
[148] P. H. Stewart, D. M. Golden, and C. W. Larson. Pressure and tem
perature dependence of reactions proceeding via a bound complex.
Combustion and Flame, 75(1), 1989.
[149] M. T. Swihart and S. L. Girshick. Thermochemistry and kinetics of
silicon hydride cluster formation during thermal decomposition of
silane. The Journal of Physical Chemistry B, 103(1):64–76, 1999.
[150] Thermodynamics Research Center, NIST Boulder Laboratories. Ther
modynamics source database, retrieved 8 May 2013. URL http:
//webbook.nist.gov.
[151] J. Troe. Theory of thermal unimolecular reactions in the falloff
range. I. Strong collision rate constants. Berichte der Bunsenge
sellschaft für physikalische Chemie, 87(2):161–169, 1983.
[152] S. Tsantilis and S. E. Pratsinis. Evolution of primary and aggregate
particlesize distributions by coagulation and sintering. AIChE Jour
nal, 46(2):407–415, 2000.
[153] S. Vemury and S. E. Pratsinis. Selfpreserving size distributions of
agglomerates. Journal of Aerosol Science, 26(2):175–185, 1995.
[154] S. Veprˇek, K. Schopper, O. Ambacher, W. Rieger, and MGJ Veprˇek
Heijman. Mechanism of cluster formation in a clean silane discharge.
Journal of The Electrochemical Society, 140:1935, 1993.
[155] C. G. Wells, N. M. Morgan, M. Kraft, and W. Wagner. A new method
for calculating the diameters of partiallysintered nanoparticles and
its effect on simulated particle properties. Chemical Engineering Sci
ence, 61(1):158–166, 2006.
[156] J. Z. Wen, M. J. Thomson, S. H. Park, S. N. Rogak, and M. F. Light
stone. Study of soot growth in a plug flow reactor using a moving
sectional model. Proceedings of the Combustion Institute, 30(1):1477–
1484, 2005.
151
[157] E. Werwa, A. A. Seraphin, L. A. Chiu, C. Zhou, and K. D. Kolen
brander. Synthesis and processing of silicon nanocrystallites using
a pulsed laser ablation supersonic expansion method. Applied Physics
Letters, 64(14):1821–1823, 1994.
[158] H. Wiggers, R. Starke, and P. Roth. Silicon particle formation by py
rolysis of silane in a hot wall gasphase reactor. Chemical Engineering
& Technology, 24(3):261–264, 2001.
[159] D. Woiki, L. Catoire, and P. Roth. Hightemperature kinetics of Si
containing precursors for ceramic processing. AIChE Journal, 43
(S11):2670–2678, 1997.
[160] H. W. Wong, X. Li, M. T. Swihart, and L. J. Broadbelt. Detailed kinetic
modeling of silicon nanoparticle formation chemistry via automated
mechanism generation. The Journal of Physical Chemistry A, 108(46):
10122–10132, 2004.
[161] J. J. Wu and R. C. Flagan. Onset of runaway nucleation in aerosol
reactors. Journal of Applied Physics, 61(4):1365–1371, 1987.
[162] J. J. Wu, R. C. Flagan, and O. J. Gregory. Submicron silicon powder
production in an aerosol reactor. Applied Physics Letters, 49(2):82–
84, 1986.
[163] J. J. Wu, H. V. Nguyen, and R. C. Flagan. A method for the synthesis
of submicron particles. Langmuir, 3(2):266–271, 1987.
[164] J. S. Wu, W. J. Hsiao, Y. Y. Lian, and K. C. Tseng. Assessment of
conservative weighting scheme in simulating chemical vapour depo
sition with trace species. International Journal for Numerical Methods
in Fluids, 43(1):93–114, 2003.
[165] Y. Xiong and S. E. Pratsinis. Formation of agglomerate particles by
coagulation and sintering–Part I. A twodimensional solution of the
population balance equation. Journal of Aerosol Science, 24(3):283–
300, 1993.
[166] M. Zachariah and M. Carrier. Molecular dynamics computation of
gasphase nanoparticle sintering: a comparison with phenomenolog
ical models. Journal of Aerosol Science, 30:1139–1152, 1999.
152
[167] R. Zauner and A. G. Jones. On the influence of mixing on crystal pre
cipitation processes–application of the segregated feed model. Chem
ical Engineering Science, 57(5):821–831, 2002.
[168] H. Zhao and C. Zheng. A new eventdriven constantvolume method
for solution of the time evolution of particle size distribution. Journal
of Computational Physics, 228(5):1412–1428, 2009.
[169] H. Zhao and C. Zheng. A population balanceMonte Carlo method for
particle coagulation in spatially inhomogeneous systems. Computers
& Fluids, 71:196–207, 2013.
[170] H. Zhao, F. E. Kruis, and C. Zheng. A differentially weighted Monte
Carlo method for twocomponent coagulation. Journal of Computa
tional Physics, 229(19):6931–6945, 2010.
153