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 sequential-modular 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 gas-phase 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 start-up of a plug-flow 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 Fredrich-Alexander-
Universität Erlangen-Nü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 Gas-phase chemical kinetics . . . . . . . . . . . . . . . . . . . 5
2.1.1 Rate expressions . . . . . . . . . . . . . . . . . . . . . 6
2.1.2 Three-body reactions . . . . . . . . . . . . . . . . . . . 7
2.1.3 Falloff reactions . . . . . . . . . . . . . . . . . . . . . . 7
2.1.4 Thermochemistry . . . . . . . . . . . . . . . . . . . . . 8
2.1.5 Gas-phase 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 Gas-particle 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 Type-space . . . . . . . . . . . . . . . . . . . . . . . . 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 Type-space . . . . . . . . . . . . . . . . . . . . . . . . 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 High-precision 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 Gas-phase 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 Unsteady-state 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 near-spherical 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 highly-spherical 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 binary-tree 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 three-reactor 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 computer-generated TEM
images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.10 Temporal evolution of particle diameter in a laser-driven 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 Gas-phase 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 parameter-fitting 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 surface-reaction pre-exponential 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 sub-divided and mul-
tiplied in size, like doomsday-bringing ‘grey goo’. One of the biggest, and
possibly most well-known 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 by-product 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 nano-sized 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’ length-scale, 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 gas-phase. 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 gas-phase 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 better-account 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 gas-phase 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 tied-together
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 start-up of a plug-flow 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 gas-phase 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 Gas-phase 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 gas-phase
chemical species and reactions will be referred to throughout this text as
the ‘gas-phase mechanism’ or the gas-phase state-space, 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 gas-phase 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
best-known form for the forward rate constant is the modified-Arrhenius
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 Three-body 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 high-pressure limit, where the rate is depen-
dent on temperature. However, three-body 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 low-pressure limit [71].
Expressions for the rate constant which span the low- and high-pressure
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 modified-Arrhenius
parameters are first defined: those for the high-pressure limit (k∞) and the
low-pressure 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 on-line [150] or in other peer-
reviewed literature (e.g. [149]). Alternatively, it is possible to derive this
data through use of quantum-mechanical 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 Gas-phase 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 gas-phase chemical kinetics library,
which supplies the rate of progress and mixture properties to an ODE solver.
This work uses an in-house built library named Sprog to determine chem-
ical reaction rates and gas-phase 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 type-space of a
8
particle, denoted as P. The collection of all particles and the information
used to describe them is called the particle state-space, 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 gas-phase species concentrations and n is
the vector of particle concentrations. Each of these processes represents a
reaction which alters type-space 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 gas-phase 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 gas-phase species (condensation) or a heterogeneous
reaction of a gas-phase 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
computationally-efficient 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 gas-phase or particle species. Kfm is the free-
molecular kernel and Ksf is the slip-flow kernel. The dominant kernel is
dependent on the Knudsen number (Kn). The free-molecular 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 mean-free-path 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 exponentially-distributed
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 cross-terms, that is, they cannot be written as linear sums of
particle properties (e.g.
√
a + b). However, as the binary-tree pre-calculates
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 free-molecular 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 over-estimation 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 single-particle 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 non-linear 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 non-trivial,
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 gas-phase
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 plug-flow 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 gas-phase reactions
ω˙k(c, T) is given by Equation (2.2). The gas-phase 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 constant-volume 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 Gas-particle coupling
Equations 2.15, 2.19, 2.25 and 2.26 imply a two-way relationship between
the particle- and gas-phase. 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 particle-phases. 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
gas-phase 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 high-precision 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
Gas-phase, 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: Highly-spherical silicon nanoparticles made in a hot-wall 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 hot-wall 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 wall-growth 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 parallel-plate 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 Si-H 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
light-emission components. Cheng et al. [29] suggested they could be intro-
duced into organic light-emitting 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
next-generation circuits or displays [7].
22
Chapter 3
Choosing a particle model
Three particle models commonly-used 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 type-space.
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 type-space of any
level of detail can be used.
Some of the most commonly-used 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 one-dimensional 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 surface-volume (surf-vol)
particle model extends the basic spherical model by tracking the volume and
total surface area of a particle [78, 165]. Another two-dimensional 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 surface-volume model
was also extended in the surface-volume-hydrogen model of Blanquart and
Pitsch [15], which considered the number of active sites as an additional
independent variable.
The surface-volume 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 binary-tree 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 surface-volume 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 over-simplification 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 Type-space
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 type-space is very simple:
Pq = Pq
(
p(q)
)
(3.2)
Here the super-script is used to denote ownership of the primary particle
p(q) by computational particle Pq. The surface-volume model tracks the
surface area S(q) in addition to the chemical composition of primaries:
Pq = Pq
(
p(q), S(q)
)
(3.3)
The surface-volume 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 binary-tree model.
This model represents a particle as a binary tree of primary particles with
connectivity stored in the non-leaf nodes:
Pq = Pq
(
p1, p2, . . . , pnq , C
(q)
)
(3.4)
where Pq contains nq primary particles and C is a lower-diagonal 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 Surf-vol Binary tree
type-space 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 type-space 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 surface-volume 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 surface-volume 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 common-surface 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 surface-volume model is an estimate based-on the surface
area and volume. The primary particle diameter for the binary-tree 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 common-surface 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 surf-vol binary-tree
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 surface-volume and binary tree
models, the exponential excess surface area decay formula–as popularised
by Koch and Friedlander [73]–is employed. In surface-volume 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 surface-volume 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 surface-volume 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 (number-based) 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 run-averages (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 based-on 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 type-space 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 surface-volume
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 surface-volume 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)
bin-tree – µg(dcol)
surf-vol
µg(dcol)bin-tree
≤ 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 near-spherical 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 surface-volume model
predictions.
33
τ/τ
C,
i, 
-
τ/τS, -
intersection
5 % error
10-2
10-1
100
101
102
10-5 10-4 10-3 10-2 10-1 100 101 102 103
Partial-sintering zone
Near-spherical zone
Spherical
zone
Aggregate
zone
Figure 3.3: A contour plot depicting where the binary tree model is equiv-
alent to the surface-volume model. The shading corresponds to
within 5% deviation (with respect to collision diameter) of the
surface-volume 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 similarly-sized particles,
the surface-volume model will under-estimate 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 free-molecular regime), thus increasing the maximum error.
In the partial-sintering zone, particles coagulate to form partially-sintered
aggregates with more than 2–3 primary particles. Here, the surface-volume
model over-predicts 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 surface-volume 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 partial-sintering and near-spherical 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 un-sintered 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 surface-volume 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 surface-volume and binary-tree 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), -
bin-tree
surf-vol
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), -
bin-tree µ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 surface-volume and binary-tree 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 binary-tree model appears to predict larger aggregates
with smaller primaries than the surface-volume model. This is attributed to
the phenomenon discussed in Section 3.2.1 for the prediction of collision
diameter in the near-spherical 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 self-preserving 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
quad-core Intel Harpertown processors and 8 GB of RAM.
In the stochastic formulation of the population balance, it is evident that
the surface-volume 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
10-5 10-4 10-3 10-2 10-1 100 101 102
t C
PU
/L
, s
τ/τS, -
Partial-sintering zone
N
ea
r-s
ph
er
ica
l
zo
n
e
Sp
he
ric
al
 z
on
e
Aggregate zone
bin-tree
surf-vol
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 binary-tree model require more time, and in this case
it requires orders of magnitude more. In the partial-sintering 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 partial-sintering 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 surface-volume model with respect to a model which tracks
full aggregate structure (binary-tree 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 surface-volume 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 well-documented
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
surface-volume model has been used, when alternative two-dimensional 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 well-characterised 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 constant-number 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 type-space 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 fully-coupled gas-phase and particle model is used to simulate the for-
mation of silica from tetraethoxysilane. The kinetic model describing the
gas-phase decomposition is described in [126, 142] and consists of 58 re-
versible gas-phase reactions among 27 chemical species.
4.1.1 Type-space
The binary-tree model presented in Section 3.1 is adapted here to the par-
ticular case of silica. A detailed mathematical formulation of this model’s
type-space 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 gas-phase 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 newly-incepted 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 point-connection [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 non-neighbouring
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 surface-rounding factor, a smoothing factor such that 0 ≤
σ ≤ 2 [137]. Cxy represents the element of the common-surface 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 gas-phase 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 viscous-flow 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.
Intra-particle 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 intra-particle 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
whole-particle 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 gas-phase must be carefully considered. A
process which removes one molecule of χ1 from the gas-phase and release
one molecule of χ2 into the gas-phase 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
Free-molecular kernel
The free-molecular 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 free-molecular 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 self-coagulations. The majorant jump rate of free-molecular co-
agulation is therefore given by:
Rˆfm =
N(t)∑
i 6=j
Kˆfm(Pi, Pj)wj
Vsmp
(4.31)
48
Slip-flow 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 slip-flow 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 slip-flow 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 free-molecular 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 free-molecular and slip-flow 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 test-case. A simple case was defined in [143] in which numerical con-
vergence studies were conducted. This case simulated a zero-dimensional
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 run-average (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 pre-exponential 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]
Gas-phase parameters - - [142]
Free-molecular 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
10-9
10-8
10-7
10-6
10-5
10-4
10-3
10-2
10-1
100
 0  2  4  6  8  10  12  14  16
ke
rn
el
 d
en
sit
y,
 1
/n
m
d pri

, nm
10-9
10-8
10-7
10-6
10-5
10-4
10-3
10-2
10-1
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 high-precision solutions at t = 0.8 s.
4.4 Results
4.4.1 High-precision 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 high-precision 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
10-9
10-8
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 lower-bound 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 high-precision 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
0-9
 
-
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: Run-averaged values of the functionals studied (Table 4.4)
and their corresponding confidence intervals at t = 0.8s with
Nmax × L = 217.
59
10-4
10-3
10-2
10-1
ε- s
ta
t(M
0) 
(t f)
, -
DSA
SWA1
SWA3
10-4
10-3
10-2
10-1
ε- s
ta
t(F
V) 
(t f)
, -
10-4
10-3
10-2
10-1
ε- s
ta
t(µ
a
(d c
o
l)) 
(t f)
, -
10-4
10-3
10-2
10-1
ε- s
ta
t(µ
a
(d p
ri

)) 
(t f)
, -
10-4
10-3
10-2
10-1
101 102 103 104 105
ε- s
ta
t(µ
a
(n)
)  (
t f),
 -
Nmax, -
10-4
10-3
10-2
10-1
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 high-precision 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 real-world coagulation problems.
This kernel is used in conjunction with a state-of-the-art 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’ high-precision 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 real-world 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 plug-flow 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 surface-volume 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 ‘cell-based weighted random
walk’ method to account for the movement of particles between cells in one-
and two-dimensional 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 three-dimensional adaptation
of a stochastic population balance solver (using an SWA [170]) coupled with
a CFD code. An algorithm was presented to de-couple 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 gas-phase 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 fully-coupled gas-phase 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 plug-flow 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 perfectly-stirred reactors (PSRs; or continuous-stirred
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 gas-phase 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 gas-phase and particle reactions, respectively. The
gas-phase 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
Type-space
The spherical type-space presented in Chapter 3 will be used for the numer-
ical and case studies of this chapter. Here, the type-space 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 type-space 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 gas-phase 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 gas-phase 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 constant-pressure 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 gas-phase
The stochastic particle method is fully-coupled to a gas-phase ODE solver.
The technique of operator splitting is employed to couple the gas- and
particle-phases. 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 ‘sequential-modular simu-
lation’, and similar approaches have been used in other applications [169].
The algorithm employed to solve the fully-coupled gas-phase 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 gas-phase over
[
tj , tj + ∆ts2
]
Solve particle-phase over
[
tj , tj +∆ts
]
Solve gas-phase 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 run-averaged 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
10-3
10-2
10-1
100
M
2(t
), a
rb.
 un
its
10-3
10-2
10-1
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 three-reactor
system (Figure 5.2(a)) with fR = 0 and fR = 0.5. Confidence
intervals of the stochastic solution are given where visible.
77
10-3
10-2
101 102 103 104 105 106
ε- (M
0),
 -
Nmax, -
DSA:Delete
DSA:Rescale
SWA:Delete
SWA:Rescale
-1 slope
10-3
10-2
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 (gas-phase 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
steady-state, 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
10-4
10-3
10-2
10-1
101 102 103 104 105 106
ε-  
(M
0),
 -
Nmax, -
DSA:Delete
DSA:Rescale
SWA:Delete
SWA:Rescale
-1 slope
10-4
10-3
10-2
10-1
101 102 103 104 105 106
ε-  
(µ a
(V
)),
 -
Nmax, -
10-3
10-2
10-1
100
101
101 102 103 104 105 106
ε-  
(M
2),
 -
Nmax, -
10-3
10-2
10-1
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
10-3
10-2
10-1
100
10-1 100 101 102
ε s
ta
t(M
3),
 -
t/τPSR, -
DSA
linear
fR=0.2fR=0.4fR=0.6fR=0.8
10-3
10-2
10-1
100
10-1 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
×10-3
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
×10-3
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 leveling-off 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 hot-wall reactor or similar
method to manufacture nanoparticles [75, 99, 141, 163]. These processes
are typically modelled using a plug-flow 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 plug-flow 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 perfectly-stirred 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 start-up 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 networked-PSR 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
fully-coupled gas-phase 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 perfectly-stirred 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 by-product 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 (doubly-bound 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 gas-phase 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],
rotating-disk 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’ three-step 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 dissociation-absorption-hydrogen 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 gas-phase 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 higher-quality 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 laser-heating 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 impressively-detailed experimental stud-
ies into nucleation of silicon in single- [103, 162, 163] and multi-stage
[135, 161] tubular hot-wall 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 narrowly-distributed silicon nanoparticles in a hot-wall 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 fine-seed 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 gas-phase 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 two-dimensional fully-coupled 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 finite-rate sintering conditions.
Silane decomposition at above-atmospheric 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 type-space 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 gas-phase ordinary differential equation solver is employed to
achieve this.
The structure of this chapter is as follows. The gas-phase and particle
models are explained in Section 6.1. The parameter estimation procedure
applied to the gas-phase 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 Gas-phase 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 well-understood that silane decomposition
proceeds through a series of intermediate gas-phase 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 low-pressure 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 pressure-dependent reactions.
The system of ordinary differential equations describing the reactions
of Table 6.1 was solved using MOPS’ ODE solver and the in-house devel-
oped chemical kinetics library Sprog for a zero-dimensional batch reactor
configuration. In order to gain better agreement with experimental data,
pre-exponential factors for five reactions were adjusted. The reasons for
this are discussed in Section 6.2.
6.1.2 Particle model
Type-space
The binary-tree 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 type-space are identical to those given for
the binary-tree particle model in Table 3.1. Here, the value of the fractal
92
Table 6.1: Gas-phase mechanism employed. Parameters in bold have been
adjusted from their original referenced values (see Section 6.2).
Pre-exponentials 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]
low-pressure parameters 4.0× 1012 0 45.1 [125]
2 Si2H6 (+M) −⇀↽ SiH4+ SiH2 (+M) 1.8× 1010 1.7 50.2 [63]
low-pressure parameters 1.8× 1035 -10.4 56.0 [63]
3 Si2H6 (+M) −⇀↽ Si2H4B + H2 (+M) 9.1× 109 1.8 50.2 [63]
low-pressure parameters 7.8× 1040 -7.8 59.0 [63]
4 Si3H8 (+M) −⇀↽ SiH2 + Si2H6 (+M) 7.0× 1012 1.0 52.7 [63]
low-pressure parameters 1.7× 1069 -15.1 60.5 [63]
5 Si3H8 (+M) −⇀↽ Si2H4B + SiH4 (+M) 3.7× 1012 1.0 50.9 [63]
low-pressure parameters 4.4× 1066 -17.3 59.3 [63]
6 Si2H4B (+M) −⇀↽ Si2H4A (+M) 2.5× 1013 0.2 5.4 [63]
low-pressure 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 free-molecular 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 gas-phase 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 type-space
for the binary-tree 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 ‘dissociation-adsorption-hydrogen’
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 state-space 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 first-order 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 pre-exponential 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 grain-boundary 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 pre-exponential 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 -
Free-molecular 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 gas-phase kinetic model and the particle population balance are cou-
pled using an operator-splitting 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 fully-coupled 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 gas-phase pre-exponentials [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 gas-phase 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 pre-exponentials of Reaction 5 and Reaction 8 (reverse) was
necessary to achieve agreement with experimental work. As silane decom-
position occurs at the low-pressure limit [125] up to 5 atm the low-pressure
pre-exponentials were adjusted, where available.
As in previous studies, the surface reaction [56, 76, 92] and hydrogen
release [114] rates were adjusted via the pre-exponential 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 experimentally-sourced 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 fourth-order 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 model-fitting. Table 6.6 lists the low-pressure 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 pre-exponential, 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 gas-phase 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 low-pressure
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 gas-phase mechanism found
through fitting (Table 6.1) provides a correct overall nucleation rate for the
experimental cases studies. On its own, the gas-phase mechanism is not
physically representative of the real decomposition rate.
After conversion from SURFACE CHEMKIN [30] to the type-space 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 gas-surface 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 hot-zone residence
time; and a constant temperature was used. PSDs are the number-density
distributions and were generated using kernel density estimation and the
normal distribution approximation to estimate the bandwidth [146].
Narrow PSDs
Narrowly-distributed 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 ‘smoothing-out’ 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 over-predicts 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% SiH4-Ar 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 highly-spherical 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 self-preserving limit.
Particle aggregates
Particle synthesis under conditions of finite-rate sintering was experimen-
tally and theoretically investigated by Nguyen and Flagan [103]. Using a
single-stage 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 theoretically-predicted 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 cold-zones of the reactor. Ag-
gregates were sampled from inside the reactor, corresponding to various
hot-zone 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 computer-generated
TEM-style 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 POV-Ray [124].
results. In order to compare particle aggregates, a TEM-style 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 laser-driven flame [22, 23,
44, 45]. These studies varied process conditions such as silane fraction and
cell pressure to produce highly spherical and loosely-agglomerated particles.
The particle size and number density was tracked through the flame using
laser-scattering measurements. To convert the temperature-length 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 laser-driven 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 temperature-time 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 fully-coupled gas-phase and particle model for the synthesis of sili-
con nanoparticles was presented. The mechanism of Ho et al. [63] was used
to describe the gas-phase 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 pre-exponential factors of five gas-phase 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 TEM-style image of silicon
aggregates in qualitative agreement with experiment.
It was hypothesised that, when combined with the inception process, the
new gas-phase mechanism with adjusted parameters accurately described
the silicon particle inception rate. It was, however, noted that gas-phase
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 finite-rate 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 knits-together
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 hot-wall 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 unsteady-state 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 binary-tree particle model, and
114
did not include factors such as a reactive gas-phase or surface growth on
particles.
Here, the analysis of Chapter 3 will be expanded by combining the spher-
ical and surface-volume 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 ‘real-life’ 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 non-linear.
Thus, the three particle models will begin to diverge when sintering moves
from effectively instantaneous (very small particles) to finite-rate. 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 surface-volume nor the spherical
particle model perform comparably to the binary-tree model. The number
density (M0) is double the binary-tree model’s prediction, and the average
collision diameter (µa(dcol)) 20 % lower. The surface-volume 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 surface-volume 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
bin-tree
surf-vol
spherical
10-9
10-8
10-7
10-6
10-5
 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
bin-tree µ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
bin-tree
surf-vol
spherical
t=0.40s
t=0.75s
t=1.0s
t=1.0s
t=0.75s
 0
 2
 4
 6
 8
 10
 12
 14
10-2 10-1 100 101 102 103
direction
of time
Spherical
zone
Near-spherical
zonePartial-sintering 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 binary-tree and surface-volume 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 coagulation-sintering 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 near-spherical
zone. The surface-volume and spherical models subsequently show a large
increase in coagulation rate, while it steadily decreases for the binary-tree
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 Unsteady-state 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 unsteady-state response of the hot-wall 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 TEM-style images at the given reactor
time (i.e. axial position) and transient response time coordi-
nates, generated using POV-Ray [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 partially-sintered primaries to form, thus depression in
the sintering level in both time trajectories is observed.
Information stored in the binary-tree particle model can be used to gen-
erate images of sample particle structures. This is illustrated in Figure 7.4,
where sample TEM-style 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 surface-volume model was used
in the place of the binary-tree 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 model-fitting with
the surface-volume type-space [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 hot-wall reactor. The transient and axial evolution of the
quantities illustrates the power of the methodology developed in this thesis;
using a fully-coupled model incorporating gas-phase 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 variance-reduction 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 fully-coupled 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 gas-phase 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 one-dimensional 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 binary-tree 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 well-documented that surface-growth can
cause rounding of aggregate particles. It could therefore be useful to
append a fourth axis representing a characteristic surface-growth rate
to Figure 3.3, which could potentially elucidate the suitability of par-
ticle models in accounting for surface-rounding. This would require
further investigation into the ‘surface smoothing factor’ (σSR) for the
binary-tree 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
binary-tree particle model is likely to be too computationally expensive
to directly couple to CFD calculations, a compartmental model could
124
be a useful middle-ground 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 well-known 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 near-spherical 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 surface-volume model may
be calculated:
dsurf-volcol =
1
2
(
6v1
pi
) 1
3
[
(1 + k)
1
3 +
√
1 + k
2
3
]
(A.1.3)
The binary-tree 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 binary-tree model
calculates the collision diameter of the aggregate particle as:
dbin-treecol =
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 surface-volume
and binary tree models as a function of relative colliding par-
ticle size.
The two models are compared by plotting the ratio of dbin-treecol to d
surf-vol
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 binary-tree particle
model to a maximum of 10% under-estimation.
A.2 Comparison of maximum sinterable surface area
Now suppose that a particle Pq consists of n equally-sized 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 surface-volume
model is the sum of the surface area of each of the primaries:
Ssurf-vol = npi
(
6v1
pi
) 2
3
(A.2.2)
In the binary-tree 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:
Sbin-tree = 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 surface-volume models, the
ratio of the maximum sinterable surface area (r) for between models is used:
r =
∆Sbin-tree
∆Ssurf-vol
=
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 surface-volume 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 gas-surface 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 type-space 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 gas-phase (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 type-space 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 pre-exponential 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
pre-exponential factor A* into the A used in present work. An upper-bound
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 pre-exponential factors calculated via optimisation (ASR,SiH4 =
3.0× 1033) in the present work.
131
Nomenclature
Upper-case Roman
A Arrhenius pre-exponential 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 Gas-phase state-space
Hˆ Specific molar enthalpy
K Coagulation kernel
Kc Chemical equilibrium constant
Kfm Free-molecular coagulation kernel
Ksf Slip-flow coagulation kernel
Ktr Transition regime coagulation kernel
Kn Knudsen number
L Number of independent runs
M Third-body collision species symbol
Mi(t) Moment i of the particle ensemble
NA Avogadro’s number
NGP Number of gas-phase chemical species
Nh Number of factorised terms in majorant kernel
Nrxn Number of gas-phase 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 state-space
Lower-case 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 Low-pressure limit rate coefficient
k∞ High-pressure 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
Upper-case Greek
Φ Objective function
Lower-case Greek
α Free-molecular coagulation parameter
β Arrhenius temperature exponent
β1,β2 Slip-flow 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
ψ Gas-phase expansion rate
ω˙ Molar production rate due to gas-phase 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 Continuous-stirred tank reactor
DSA Direct simulation algorithm
MD Molecular dynamics
ODE Ordinary differential equation
PFR Plug-flow reactor
PSD Particle size distribution
PSR Perfectly-stirred reactor
SWA Stochastic weighted algorithm
135
Bibliography
[1] A. J. Adamczyk, M. F. Reyniers, G. B. Marin, and L. J. Broadbelt. Ex-
ploring 1, 2-hydrogen 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 DQMoM-IEM 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 CFD-population 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. Schmidt-Ott. Generation and
sizing of particles for aerosol-based 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 volume-surface-hydrogen 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 laser-driven 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 laser-driven reactions: II, Powder char-
acteristics and process variables. Journal of the American Ceramic
Society, 65(7):330–335, 1982.
[24] Cantera Developers. Cantera: An object-oriented 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 gas-phase 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. Age-related 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 light-emitting 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 Solid-Surface-Gas-Phase Interface. Sandia National
Laboratories, 1996.
[31] H. Dang and M. T. Swihart. Computational modeling of silicon
nanoparticle synthesis: I. A general two-dimensional 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 two-dimensional bivariate model for sil-
icon nanoparticle synthesis in a laser-driven 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: Not-so-new 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. Gas-phase decomposition of silane in a horizontal
epitaxial reactor. Philips Res. Repts, 26(2):134–144, 1971.
[42] R. Fenollosa, F. Ramiro-Manzano, 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 laser-heated 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 laser-driven 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. Chapple-Sokol, and R. G. Gordon.
Gas-phase 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 nano-particle 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
liquid-liquid 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 Quasi-Monte 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. Multi-Monte Carlo ap-
proach for general dynamic equation considering simultaneous parti-
141
cle coagulation and breakage. Powder Technology, 154(2-3):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(4-6):
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. Laser-induced 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 hot-wall tubular reactor. Materials
Science and Engineering: B, 17(1-3):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 high-pressure laminar premixed flames.
Combustion and Flame, 114(3):484–501, 1998.
[71] Robert J Kee, Fran M Rupley, Ellen Meeks, and James A Miller.
CHEMKIN-III: A FORTRAN chemical kinetics package for the analysis
of gas-phase 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 cell-based weighted random walk method and a constant-number
Monte-Carlo 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
compartment-based 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 laser-driven pyrolysis of silane followed
by HF-HNO3 etching. Langmuir, 20(11):4720–4727, 2004.
[86] Z. F. Li and E. Ruckenstein. Water-soluble 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. Laser-induced
chemical vapor deposition of hydrogenated amorphous silicon. I.
Gas-phase process model. Journal of Applied Physics, 62(7):2812–
2821, 1987.
[96] M. Meunier, J. H. Flint, J. S. Haggerty, and D. Adler. Laser-induced
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. Gas-phase
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 metal-oxide-semiconductor 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
single-stage 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 low-pressure 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 particle-dependent 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 free-space 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 Non-Crystalline 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 fixed-sectional 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(1-2):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 high-temperature
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. First-principles 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 particle-resolved 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 PAH-PP 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
nanometer-sized 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. Murray-Rust. First-principles 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 fall-off
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
particle-size distributions by coagulation and sintering. AIChE Jour-
nal, 46(2):407–415, 2000.
[153] S. Vemury and S. E. Pratsinis. Self-preserving 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 partially-sintered 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. High-temperature 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 two-dimensional 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
gas-phase 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 event-driven constant-volume 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 balance-Monte 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 two-component coagulation. Journal of Computa-
tional Physics, 229(19):6931–6945, 2010.
153