Global sensitivity analysis of a model for silicon nanoparticle synthesis

Global sensitivity analysis of a model for silicon nanoparticle synthesis

Journal of Aerosol Science 76 (2014) 188–199 Contents lists available at ScienceDirect Journal of Aerosol Science journal homepage:

1MB Sizes 0 Downloads 33 Views

Journal of Aerosol Science 76 (2014) 188–199

Contents lists available at ScienceDirect

Journal of Aerosol Science journal homepage:

Global sensitivity analysis of a model for silicon nanoparticle synthesis William J. Menz, George P.E. Brownbridge, Markus Kraft n Department of Chemical Engineering and Biotechnology, University of Cambridge, New Museums Site, Pembroke Street, Cambridge CB2 3RA, United Kingdom

a r t i c l e in f o


Article history: Received 29 January 2014 Received in revised form 23 June 2014 Accepted 27 June 2014 Available online 7 July 2014

This paper presents a global sensitivity analysis of the detailed population balance model for silicon nanoparticle synthesis of Menz & Kraft [2013a, A new model for silicon nanoparticle synthesis, Combustion & Flame, 160:947–958]. The model consists of a gas-phase kinetic model, fully coupled with a particle population balance. The sensitivity of the model to its seven adjusted parameters was analysed in this work using a High Dimensional Model Representation (HDMR). An algorithm is implemented to generate response surface polynomials with automatically selected order based on their coefficient of determination. A response surface is generated for 19 different experimental cases across a range of process conditions and reactor configurations. This enables the sensitivity of individual experiments to certain parameters to be assessed. The HDMR reveals that particle size was most sensitive to the heterogeneous growth process, while the particle size distribution width is also strongly dependent on the rate of nucleation. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Silicon Parameter estimation High dimensional model representation Population balance

1. Introduction Knowledge of the properties of silicon has underpinned the revolution of information technology. A great deal of study has therefore been undertaken into its synthesis, purification and resultant properties. Silicon nanoparticles were first made in the late 1970s (Murthy et al., 1976) and since then, have been the subject of considerable work towards improving their manufacture and identifying potential applications. Gas-phase, laser and plasma synthesis of particles are the most common methods with which silicon nanoparticles are manufactured (Mangolini et al., 2013). In general, these processes begin with silane (SiH4), which may be decomposed by thermal, laser or microwave radiation (Cannon et al., 1982; Knipping et al., 2004; Nguyen & Flagan, 1991). The decomposition of silane forms reactive silicon hydrides, which combine with each other to nucleate into silicon nanoparticles (Swihart & Girshick, 1999). Various models—from gas-phase kinetic to particle population balances—have been proposed to describe the formation of nanoparticles from silane (Menz et al., 2012; Onischuk et al., 2000; Swihart & Girshick, 1999). In all of these models, it is common to find model parameters which have uncertain values. Examples include empirical expressions for sintering (Chen et al., 2013; Menz et al., 2012; Shekar et al., 2012).


Corresponding author. Tel.: þ44 1223 762784; fax: þ44 1223 334796. E-mail address: [email protected] (M. Kraft). URL: (M. Kraft). 0021-8502/& 2014 Elsevier Ltd. All rights reserved.

W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199


In some cases, it is possible to estimate these values by trial-and-error (Körmer et al., 2010). However, for a non-linear model with many unknown parameters, this process becomes difficult. Systematic parameter estimation can be used to s move from by-hand guesswork to computer-based solutions. For example, a model implemented in Excel could use the GOALSEEK or SOLVER functionality to arrive at better input parameter values. There are a range of optimisation techniques which can be used to improve model values (Kraft & Mosbach, 2010). Lowdiscrepancy sampling can be used to evaluate the model response over a space of parameter values (Braumann et al., 2011). Gradient-search methods such as the simultaneous perturbation stochastic approximation (SPSA) algorithm can locate local and global minima (Chen et al., 2013; Menz et al., 2012). Response surfaces, or surrogate models, can also be generated from low-discrepancy sampling, yielding a computationally efficient approximations of the true model (Kastner et al., 2013). Then, Markov Chain Monte Carlo (MCMC) sampling and Bayesian analyses are often applied to assess the credible regions in which the optimal parameter values may lie (Kastner et al., 2013; Mosbach et al., 2012). Recently, a detailed model for silicon nanoparticle synthesis was presented (Menz & Kraft, 2013a). The model incorporates a gas-phase kinetic model, fully coupled with a particle population balance. In its original development, systematic parameter estimation was used to optimise the model's seven parameters with respect to experimental data across a range of different process conditions. However, the relative importance of each parameter in the objective function was not quantitatively addressed in this work. Nor was the influence of particular experiments on the objective function. These open questions should be addressed in order to further assess the physical relevance of the model as well as its sensitivity to input parameters. Due to the variety of solution methods for population balance models, there are a range of different approaches for studying the sensitivity of the model's parameters. Vikhansky & Kraft (2004, 2006) demonstrated use of a gradient search method to assess the sensitivity of stochastic (Monte Carlo) solution methodologies as applied to population balance equations. In complex models, a simple scan of the parameter space can also be used (Lavvas et al., 2011; Shekar et al., 2012). An excellent review of sensitivity analysis methods is given by Tomlin (2013). Of these methods, the high dimensional model representation (HDMR) is one which is yet to be applied to population balance modelling. In creating a HDMR of a model, one obtains a response surface for each of the supporting experiments of the model. While potentially useful for a low-computational effort optimisation, the surfaces serve a dual purpose: their coefficients reveal the model response's sensitivity to its input parameters. This information can then be used to refine the model or gain further insight into the factors which most strongly drive the results. HDMR has further benefits over conventional brute-force approaches for sensitivity analysis, as it averages over many points—this reduces the method's sensitivity to local variations in gradient (Azadi et al., 2014). The purpose of this work is to investigate the estimated parameters in the model of Menz & Kraft (2013a) for silicon nanoparticle synthesis. A global sensitivity analysis will be conducted using HDMR response surfaces in order to elucidate the influence of each parameter in the optimisation. This will illustrate an alternative approach through which sensitivity analyses can be conducted, not as yet used in this community. Finally, the sensitivity analysis will also be applied to gain additional physical insight into the model. The structure of this paper is as follows. A brief description of the model is given in Section 2, including the gas-phase (Section 2.1) and the population balance (Section 2.2). The techniques used for parameter estimation and sensitivity analysis are presented in Section 3. The results from the sensitivity analysis are given and discussed in Section 4.

2. Model The model for silicon nanoparticle synthesis is composed of a gas-phase kinetic model and a particle population balance model. A full formulation of the model is given in Menz & Kraft (2013a), however, a brief description is given here. In the original development of the model, parameter estimation was used to optimise gas-phase and particle-phase parameters with respect to experimental results from the literature. This process is illustrated in Fig. 1.

2.1. Gas-phase model The bulk decomposition of silane can be described as the bimolecular expression (Petersen & Crofton, 2003): SiH4 þ M-SiH2 þ H2 þM


where M is a third body. It is well-understood that silane decomposition proceeds through a series of intermediate gasphase species, such as silylene (SiH2) and higher silenes/silanes. For this model, the mechanism of Ho et al. (1994) is adopted. The mechanism has eight core reactions, mostly described by Lindemann falloff expressions. Pre-exponential factors for five of the reactions were adjusted from their initial values. The equations describing the rate of change of chemical species due to these reactions are solved using a conventional ordinary differential equation (ODE) solver.

W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199

gas−phase clustering

H3SiSiH + SiH 4 A A6 −SiH 2 8 + H2 H2SiSiH2 A3 Si 2H6

+ SiH 4

experimental data SiH

4 thermal decomposition 2

x x

Si 3H8

A SR,SiH 4 surface reaction


PSDs f(T, P, y SiH ,τ)

x x


particle population balance




new parameter estimation methodology

+SiH 4

+SiH 2




+SiH 2 (rev) A 4


x x x x

inception hydrogen release A H2


silicon nanoparticle synthesis model −H2 A 1 SiH 2 (rev) A 7 + SiH 4 −H2 +SiH 2 (rev) A


new model results frequency


PSDs coagulation

sintering diameter

Fig. 1. Schematic of the information transfer among the models, the parameter estimation methodology and experimental sources. The quantities highlighted in red refer to the adjusted parameters. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

2.2. Particle population balance model 2.2.1. Population balance The nucleation and growth of particles is 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 mathematically describe a particle. This representation is called the type-space of a particle, denoted as P. The creation and interaction of particles with each other in these reactors is described by a population balance equation. Excluding transport of particles to-and-from the reactor, the evolution of the number density of particles nðP q Þ in ensemble vector n of type Pq is given by           d   n P q ¼ I c; P q þ K n; P q þ R c; n; P q þS n; P q þ ψ ðc; n; T Þn P q dt


where the operator I ðc; P q Þ defines the inception and nucleation of particles, Kðn; P q Þ is the coagulation operator, Rðc; n; P q Þ is the operator for surface growth or oxidation and Sðn; P q Þ the operator for sintering. The gas-phase expansion factor is given by symbol ψ and the gas-phase species concentrations are given by vector c. Each of the operators contains a set of processes which may create, destroy and transform particles. For example, operator Rðc; n; P q Þ would contain a surface growth process describing the transformation P q -P r by adding a gas-phase component from c. Further definition of these quantities is given in Celnik et al. (2007) and Menz et al. (2014). A stochastic numerical method is employed to solve the population balance model. Stochastic methods solve the population balance equations by modelling how the particle processes alter an ensemble of computational particles. This method has been well-documented (Balthasar & Kraft, 2003; Morgan et al., 2006; Shekar et al., 2012) and takes advantage of various features such as linear process deferment (Patterson et al., 2006), a binary tree cache (Goodson & Kraft, 2002) and the concept of majorant rates with fictitious jumps (Eibeck & Wagner, 2001) to accelerate the solution of the population balance model. Coupling with the gas-phase is accomplished using the technique of operator splitting (Celnik et al., 2007; Shekar et al., 2012). Performance of the model and its numerical convergence behaviour have been addressed in Celnik et al. (2007), Shekar et al. (2012), and Menz et al. (2013, 2014). 2.2.2. Particle model The binary-tree particle model of Sander et al. (2009) is used in this work. Each particle in the ensemble, Pq, is represented as P q ¼ P q ðp1 ; …; pnq ; CÞ


where particle Pq contains nq primary particles px. C is a lower-diagonal matrix representing the common surface area between two primary particles. Each primary particle is described by the number of silicon atoms ηSi and hydrogen atoms ηH : px ¼ px ðηSi ; ηH Þ


W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199


Table 1 Particle processes included in the population balance model. Process


Inception Surface reaction

Sii Hj B þ Sik Hl -PðηSi ¼ i þ k; ηH ¼ j þ l; CÞ Pð…; px ðηSi ; ηH Þ; …; CÞ þ Sii H2i þ 2 -Pð…; px ðηSi þ i; ηH þ 2Þ; …; C0 Þ þ iH2

Condensation ðj a iþ 2Þ

Pð…; px ðηSi ; ηH Þ; …; CÞ þ Sii Hj

Hydrogen release

Pð…; px ðηSi ; ηH Þ; …; CÞ -Pð…; px ðηSi ; ηH  2Þ; …; CÞþ H2

-Pð…; px ðηSi þ i; ηH þ jÞ; …; C0 Þ

Coagulation P q ðp1 ; …; pnq ; Cq Þþ P r ðp1 ; …; pnr ; Cr Þ -P s ðp1 ; …; pnq ; pnq þ 1 ; …; pnq þ ns ; Cs Þ Pðp1 ; …; pn ; CÞ-Pðp1 ; …; pn0 ; C0 Þ

Sintering ðn0 r nÞ

All other properties of the particle are derived from this representation. A full mathematical formulation of the model is given in Menz & Kraft (2013a) and Shekar et al. (2012). Particle processes can change the type-space of a particle in a variety of different ways. The processes included in the model for silicon nanoparticle synthesis are listed in Table 1 and described below: Inception: A collision of a silylene (e.g. SiH2) species with another silylene or silene species will form a particle. The rate of inception dependent on the critical nucleus diameter, a quantity which may be estimated using macroscopic properties of silicon and the process conditions (Menz & Kraft, 2013a). When the diameter of the particle to be incepted is greater than the critical diameter, the rate of inception is calculated using the transition regime coagulation kernel (Menz et al., 2013). Surface reaction: In this context, surface reactions refer to the heterogeneous reaction of silanes (SiH4, Si2H6, Si3H8) on the particle surface. Silanes must proceed through an energy barrier in order for the reaction to proceed (Houf et al., 1993). Hence, the rate of surface reaction is represented by an Arrhenius rate constant (Ho et al., 1994). It is also proportional to the particle surface area (Menz & Kraft, 2013a). In the population balance model, a surface reaction event also causes rounding of joined primaries. This is represented by the transformation C-C0 . Condensation: Silenes (e.g. H2SiSiH2) and silylenes (e.g. SiH2) react with particle surfaces through a condensation process. That is, their rate of reaction is given by a collision kernel, here the free-molecular kernel. It is assumed that they stick with probability 1.0 and that there is no energy barrier to reaction. Hydrogen release: Particles must release hydrogen in order to maintain a stable crystal structure when growing. The rate of hydrogen desorption is proportional to the coverage of hydrogen on the particle surface as well as an Arrhenius rate constant (Sinniah et al., 1990). In this model, the coverage is approximated by the ratio of number of hydrogen atoms to number of silicon atoms in each particle. Coagulation: The rate of coagulation is dependent on the Knudsen number, and is given by the transition kernel. In a coagulation event, the particle trees are added to each other, preserving the particle size and connectivity information held in the state-space. Sintering: Adjacent primary particles may sinter through a grain-boundary diffusion process. The common surface element of connected primaries px and py, Cxy, will decrease until it is sufficiently close to the equivalent spherical area of the two primaries. Then, the primaries are merged into a single primary particle. This occurs as a continuous process throughout the tree of primaries in the computational particle. 3. Parameter estimation In the original development of the model, parameter estimation was used to obtain suitable fit of the model response to experimental results. Seven parameters were identified for estimation. The physical and the numerical meaning of these parameters are given in Fig. 1 and Table 2 respectively. This results in a parameter vector given by

θ ¼ ðA1;LP ; A2;LP ; A3;LP ; A5;LP ; A8;rev ; ASR;SiH4 ; AH2 Þ n


The optimal set of parameters ðθ ¼ ðA1;LP ; …; AH2 Þ, the star denoting optimal) is reliant on measuring the distance of the model response from the experimental results. To do this, a least-squares objective function, ΦðθÞ, was formulated: Nexp

2 ΦðθÞ ¼ ∑ ðηexp  ηsim i ðθÞÞ i i¼1





W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199

Table 2 Overview of adjusted parameters. j




1 2 3 4 5 6 7

A1;LP A2;LP A3;LP A5;LP A8;rev ASR;SiH4 AH2

Gas Gas Gas Gas Gas Particle Particle

Low-pressure pre-exponential factor for Reaction Low-pressure pre-exponential factor for Reaction Low-pressure pre-exponential factor for Reaction Low-pressure pre-exponential factor for Reaction Reverse pre-exponential factor for Reaction 8 Pre-exponential factor for SiH4surface reactions Pre-exponential factor for H2 release

1 2 3 5

where ηexp is an experimental response, and ηsim is obtained from the model for the corresponding experiment. The number i i of experiments used in the optimisation is given by Nexp . The parameters were ‘coded’ using logarithmic coordinates. Such transformations are common to simplify calculations (Kastner et al., 2013; Mosbach et al., 2014), particularly where the bounds on the parameters can span several orders of magnitude:   log θ  12 log θub þ log θlb 0   θ ¼ ð7Þ 1 2 log θ ub log θ lb

θ0 ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi logðθ= θub θlb Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðlog θub =θlb Þ

ð8Þ 0

where θub is the upper bound of the parameter, θlb is the lower bound, θ is the actual value (in original units) and θ is the 0 coded value ðθ A ½  1; 1Þ. A staged optimisation procedure was used to find the optimal set of parameters θ^ , with numerical results given in Menz & Kraft (2013a). In that work, a least-squares objective function was used to select parameters, based on specific experimental studies. In the present work, a least-squares objective function is also used; however, a larger set of experimental values will be used, and an objective function including a measure of the particle size distribution (PSD) width will be considered. This will extend the previous work's analysis by quantifying the sensitivity of parameters on model response. 3.1. High dimensional model representation In this work, we are interested in the sensitivity of the model and its parameters near the optimal found in previous work. The global sensitivities can be calculated from a High Dimensional Model Representation (HDMR) (Rabitz & Aliş, 1999). A HDMR is a form of response surface, or surrogate model. Its main feature is the decomposition of the full function for a response ηsim into a sum of functions that only depend on the input variables such that Nparam

Nparam Nparam

ηsim ðθÞ ¼ f ðθÞ ¼ f 0 þ ∑ f i ðθi Þ þ ∑ i¼1

∑ f ij ðθi ; θj Þ þ ⋯ þf 12‥Nparam ðθ1 ; θ2 ; …; θNparam Þ


i ¼ 1 j ¼ iþ1

where Nparam is the number of parameters, i and j index the input parameters and f0 is the mean value of f ðθÞ. In practical applications, it is possible to truncate the above expression to second-order terms (Rabitz & Aliş, 1999). The order of a group of terms is given the symbol d, with the maximum order considered dmax . Thus, a response is approximated by Nparam

ηsim ðθÞ  f ðθÞ ¼ f 0 þ

Nparam Nparam

∑ f i ðθ i Þ þ ∑ ∑ f ij ðθi ; θj Þ i¼1 i ¼ 1 j ¼ iþ1 |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

firstorder terms


secondorder terms

For example, if one were to write the form of a response surface with maximum interaction order dmax ¼ 2, and maximum polynomial orders M max;d ¼ 1 ¼ 2, M max;d ¼ 1 ¼ 1, the functional form could appear similar to f ðθÞ ¼ f

ðd ¼ 0Þ

N param

ðd ¼ 1;M ¼ 1Þ

þ ∑ fi i¼1

N param N param

¼ 1;M ¼ 2Þ 2 θi þf ðd θi þ ∑ i

The challenge is to then determine the coefficients f

i ¼ 1 j ¼ iþ1



ðd ¼ 2;M ¼ 1Þ

þ f ij

ðd ¼ 1;M ¼ 1Þ



ðd ¼ 1;M ¼ 2Þ

ð11Þ ,f

ðd ¼ 2;M ¼ 1Þ


3.1.1. Calculation of sensitivities The global sensitivities can be determined by considering the variance of the model response. The first-order variance, σ 2ηsim ;i , and the second-order variance, σ 2ηsim ;ij , both contribute to this quantity, given by Z 1 σ 2ηsim ;i ¼ f i ðθi Þ2 dθi ð12Þ 1

W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199

σ 2ηsim ;ij ¼


1 1


1 1


f ij ðθi ; θj Þ2 dθi dθj


where the integrals are taken over the upper- and lower-bounds chosen for the parameters, here assumed to be normalised to be ½  1; 1. This is consistent with other experimental design literature (Braumann et al., 2010; Man et al., 2010). The total variance is then given by Nparam

N param N param

σ 2ηsim ¼ ∑ σ 2ηsim ;i þ ∑ i¼1

∑ σ 2ηsim ;ij


i ¼ 1 j ¼ iþ1

Then, the first- and second-order normalised sensitivities are given by   σ 2ηsim ;i Sηsim θi ¼ 2

σ ηsim


  σ 2ηsim ;ij Sηsim θi ; θj ¼ 2


σ ηsim

That is, Sηsim ðθi Þ represents the first-order sensitivity of model response ηsim to parameter θi, and so on. These quantities can be used to compare the influence of parameters within a particular model response. To compare across model responses, the sensitivity should be converted to an absolute sensitivity. Here, the absolute sensitivity is measured by the standard deviation of the first-order model response: Snηsim ðθi Þ ¼ σ ηsim ;i


where Sηsim ðθi Þ has dimensions ½ηsim . Physically, this quantity represents the standard deviation in output response ηsim achievable from varying parameter θi within its specified bounds. n

3.1.2. Automatic order selection The functions f ðθÞ must be determined in order to evaluate the sensitivities. These functions are decomposed into orthogonal basis functions (Li et al., 2002), which are taken as Legendre polynomials. The coefficients of these polynomials are determined through a recursive fitting procedure. This has the advantage of being able to tailor the polynomial order to each of the interaction order functions. Similar approaches have been adopted in the work of Ziehn & Tomlin (2008). The algorithm to conduct this whole fitting process is given in Fig. 2. The user specifies a maximum interaction order ðdmax Þ and a maximum polynomial order for each interaction order ðM max;d Þ and the algorithm yields the first polynomial 2 2n 2 2 with R above a tolerance R . A tolerance on the minimum improvement in R , R min is applied to increase the efficiency of


Fig. 2. Algorithm for recursive fitting of the polynomials and selection of R .


W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199

Table 3 Experimental datasets used in the present work. dtype refers to the physical quantity being measured from the system: either primary ðdpri Þ or mobility ðdmob Þ diameter. The statistical quantity measured for the particle size is given in the column μtype . An asterisk (n) denotes an estimated quantity. i (–)

T (1C)

ySiH4 (%)

P (kPa)

τ (ms)

dtype (–)

μtype (–)

ηexp;μ i

ηexp;σ i


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

1100 1100 1100 1100 1100 1100 1100 900 1000 816 1047 1307 1200 1050 1150 1000 800 800 600

0.04 0.04 0.125 0.128 0.02 0.08 0.04 0.04 0.04 0.033 0.033 0.033 0.01 0.214 0.09 0.06 0.001 0.0004 0.05

2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 51 51 51 101 20 20 20 101 101 39

80 192 192 80 80 80 420 420 420 2600 2100 1800 1000 6 15 53 900 900 870

dpri dpri dpri dpri dpri dpri dpri dpri dpri dpri dpri dpri dmob dpri dpri dpri dmob dmob dpri

Mode Mean Mean Mode Mode Mode Mode Mode Mode Mode Mode Mode Mode Mean Mean Mean Mode Mode Mean

26.7 26.0 38.0 31.0 41.0 24.0 32.5 21.2 28.5 11.0 11.0 15.0 127 43.4 55.4 23.0 89.0 51.0 52.0

1.07 1.11 1.28 1.3 1.06 1.25 1.10 1.11 1.11 1.50n 1.55n 1.54n 1.38n – – – 1.61n 1.55n –

Körmer et al. (2010) Körmer et al. (2010) Körmer et al. (2010) Körmer et al. (2010) Körmer et al. (2010) Körmer et al. (2010) Körmer et al. (2010) Körmer et al. (2010) Körmer et al. (2010) Frenklach et al. (1996) Frenklach et al. (1996) Frenklach et al. (1996) Wu et al. (1987) Flint et al. (1986) ‘630S’ Flint et al. (1986) ‘631S’ Flint et al. (1986) ‘654S’ Nguyen & Flagan (1991) Nguyen & Flagan (1991) Onischuk et al. (2000)

the process. If a new iteration is found to be worse than the previous one, the new iteration is discarded and the algorithm moves onto the next interaction order. 3.2. Analysis procedure There is a wealth of experimental conditions at which the silicon nanoparticle synthesis model can be evaluated (Flint et al., 1986; Frenklach et al., 1996; Körmer et al., 2010; Nguyen & Flagan, 1991; Wu et al., 1987). These are given in Table 3. There, experiments are described by their process conditions: representative temperature (T), initial mole fraction of silane ðySiH4 Þ, total pressure (P) and residence time (τ). The resultant PSDs from such experiments, where available, are described exp;μ σ by a modal or average particle size ðηi Þ, and a geometric standard deviation ðηexp; Þ. i A low-discrepancy (Sobol) sequence was used to evaluate the model at points within the bounds defined for the parameters. Choice of bounds was guided by ranges quoted in previous work and by ensuring that the values for the sensitivity do not change when increasing the width of the intervals. This sequence was constructed around the parameter sim;μ optimals identified in previous work (Menz & Kraft, 2013a). The modal or mean particle size ðηi Þ and the geometric sim;σ standard deviation ðηi Þ of the ensemble were calculated at each of these points. HDMR response surfaces were determined from this data for the two PSD descriptors. 4. Results and discussion The HDMR approach to parameter analysis yields the sensitivity of each experimental case to each parameter (first-order) or parameter pair combination (second-order). Only the first-order sensitivities will be discussed here, as it was observed that the only significant second-order interactions were composed of those which were significant in the first-order. 4.1. Breakdown of sensitivities In Fig. 3, the sensitivity of the mode or mean of a PSD is displayed. As there were 19 experimental cases, each with 7 firstorder sensitivities, a representative experiment was chosen from each experiment group to simplify visualisation of data. The absolute sensitivities of the PSD mode or mean are first compared across the parameters for each of the 6 groups in Fig. 3(a). It is clearly evident that the particle size is most sensitive to ASR;SiH4 , that is, the surface reaction process. It has the largest absolute sensitivity for each of the representative experimental cases. This confirms the importance of accurately describing surface reaction process in models for silicon nanoparticle formation. The second most sensitive parameter is A1;LP , the silane decomposition pre-exponential. This is not surprising, as it has often been targeted for optimisation in previous studies (Körmer et al., 2010; Menz et al., 2012). In Fig. 3(b), the first-order sensitivities are compared within experiments. Again, it is evident that the surface reaction pre-exponential is the dominant parameter for all cases. However, the relative importance A1;LP and to a lesser extent, A3;LP and AH2 are seen here. It is interesting to observe that AH2 accounts for approximately 5% of the sensitivity for the cases of

W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199


Fig. 3. Comparison of the sensitivity of the PSD mode/mean (μ case) across parameter types and different representative experiments. (a) The absolute sensitivities of the PSD mode compared across parameter types. (b) The normalised sensitivities of the PSD mode compared across experiments.

Frenklach et al. (1996) and Onischuk et al. (2000), as both cases were experimentally reported to yield particles with high concentrations of hydrogen. The sensitivities are broken-down differently when considering sensitivity of the PSD's geometric standard deviation, shown in Fig. 4. Again, the model is most sensitive to ASR;SiH4 , A1;LP and less so A3;LP , AH2 , depicted in Fig. 4(a). However, the internal breakdown of sensitivities is quite different. It is illustrated in Fig. 3(b) that A1;LP plays an important role for all experimental cases. Since A1;LP has a large impact on the inception process of particles, this observation reinforces the


W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199

Fig. 4. Comparison of the sensitivity of the PSD geometric standard deviation (σ case) across parameter types and different representative experiments. (a) The absolute sensitivities of the PSD geometric standard deviation compared across parameter types. (b) The normalised sensitivities of the PSD geometric standard deviation compared across experiments.

conclusions of Nguyen & Flagan (1991) and Kómer et al. (2010) that controlling the nucleation rate is critical in controlling the dispersion of particles. Further, the absolute sensitivities with respect to geometric standard deviation are generally greatest for the cases of Flint et al. (1986), Nguyen & Flagan (1991) and Onischuk et al. (2000) (Fig. 4(a)). These three cases correspond to those where aggregate particles are formed, where the coagulation rate is likely to be high. In such systems, it is common to see a spike in geometric standard deviation prior to it reaching its final value—either controlled by a finite sintering rate or a primary coalescence (Menz & Kraft, 2013b). It is therefore suggested that by varying the parameters describing inception,

W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199


Fig. 5. Global sensitivities (Sn) for A1;LP , A3;LP , ASR;SiH4 , AH2 shown as a function of process conditions. Each circle represents an experimental case (Table 3) and its diameter is proportional to the absolute sensitivity.


W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199

the PSD's geometric standard deviation has been observed at different points along this process, leading to an increased sensitivity of the distribution width for aggregate-forming cases. Finally, the case with the least sensitive geometric standard deviation is that of Wu et al. (1987). In this case, spherical particles are produced in a very hot reactor. It is likely that due to fast sintering (Menz & Kraft, 2013b) and strong particle nucleation (Menz & Kraft, 2013a), the system reaches the self-preserving PSD for the majority of parameter combinations. Thus, while the modal point of the PSD is very sensitive to surface growth (Fig. 3(a)), the comparative width of the distribution is relatively static (Fig. 3(a)). 4.2. Sensitivities as a function of process conditions In the previous section, the sensitivity the model at specific experimental conditions was addressed. The HDMR results can also be used to investigate how sensitivity varies with the process conditions across all experiments. In order to do this, the absolute sensitivities were plotted on a bubble plot in Fig. 5, where the size of the circle is proportional to the sensitivity. The left column of Fig. 5 shows the sensitivity of the particle size descriptor, while the sensitivity of the geometric standard deviation is given in the right column. The choice of temperature and initial silane pressure ðP ySiH4 Þ was made to clearly separate the experimental cases, although in practice the experiments are actually functions of more variables—such as residence time and reactor configuration—than these two. In any case, the clear importance of the surface reaction process at all process conditions is again demonstrated here. It appears that A3;LP is generally more important at lower pressures. This could indicate that the reaction pathway for forming Si2H6 becomes more important as silane pressure decreases. This observation would be consistent with homogeneous nucleation theory, as the critical cluster size increases as the partial pressure of silicon decreases. That is, larger gas-phase clusters (e.g. Si2H6 þ) must be formed before particles are stable enough to grow. Hence, the parameter describing one of the major pathways for this process has a stronger effect on PSD characteristic at lower silane pressure. There are also some trends within experimental cases. For example, it is evident in the cases of Kö rmer et al. (2010) that the sensitivity with respect to geometric standard deviation is roughly inversely proportional to silane pressure, while direct proportionality is observed for PSD mode sensitivity. This phenomenon is consistent with the experiments of Körmer et al. (2010), where increasing the pressure of silane would cause runaway nucleation, yielding broader PSDs and larger particles. Thus, as the pressure increases and the PSD broadens, the system nears the self-preserving distribution and the PSD's geometric standard deviation becomes insensitive to input parameters. 5. Conclusions This paper has presented an adaptation of a High Dimensional Model Representation (HDMR) to a detailed model for silicon nanoparticle synthesis. The HDMR yielded information about the sensitivity of the model outputs to various input parameters across a range of experiments conditions and configurations. The output responses investigated were the particle size and the distribution width, the latter represented as geometric standard deviation. It was found that the particle size was most sensitive to the parameter controlling the surface reaction rate. The same was observed for the geometric standard deviation, however, it was also sensitive to the initial decay rate of silane. The sensitivity analysis data could also be linked to phenomena reported experimentally. For example, the model was shown to be sensitive to the hydrogen release rate in cases where high concentrations of hydrogen were experimentally reported. The knowledge of the model's sensitivity can now be applied to embark on an experimental discrimination pathway. This can potentially indicate which experiments pull parameter optimals in which directions. Alternatively, it can show which experiments—or the model's representation of them—are incorrect. Although these aspects remain to be completed, it has been demonstrated that use of a HDMR as part of parameter estimation procedure can reveal new layers of insight from a complex model.

Acknowledgements W.J.M. acknowledges financial support provided by the Cambridge Australia Trust. M.K. gratefully acknowledges the DFG Mercator programme and the support of CENIDE at the University of Duisburg Essen. This paper was presented and benefited from discussions at PBM2013, Bangalore, India. References Azadi, P., Brownbridge, G.P.E., Mosbach, S., Smallbone, A., Bhave, A., Inderwildi, O., & Kraft, M. (2014). The carbon footprint and non-renewable energy demand of algae-derived biodiesel. Applied Energy, 113, 1632–1644. Balthasar, M., & Kraft, M. (2003). A stochastic approach to calculate the particle size distribution function of soot particles in laminar premixed flames. Combustion and Flame, 133(3), 289–298. Braumann, A., Kraft, M., & Mort, P.R. (2010). Parameter estimation in a multidimensional granulation model. Powder Technology, 197(3), 196–210. Braumann, A., Man, P.L.W., & Kraft, M. (2011). The inverse problem in granulation modeling—Two different statistical approaches. AIChE Journal, 57(11), 3105–3121.

W.J. Menz et al. / Journal of Aerosol Science 76 (2014) 188–199


Cannon, W.R., Danforth, S.C., Flint, J.H., Haggerty, J.S., & Marra, R.A. (1982). Sinterable ceramic powders from laser-driven reactions: I, Process description and modeling. Journal of the American Ceramic Society, 65(7), 324–330. Celnik, M., Patterson, R.I.A., Kraft, M., & Wagner, W. (2007). Coupling a stochastic soot population balance to gas-phase chemistry using operator splitting. Combustion and Flame, 148(3), 158–176. Chen, D., Zainuddin, Z., Yapp, E., Akroyd, J., Mosbach, S., & Kraft, M. (2013). A fully coupled simulation of PAH and soot growth with a population balance model. Proceedings of the Combustion Institute, 34, 1827–1835. Eibeck, A., & Wagner, W. (2001). An efficient stochastic algorithm for studying coagulation dynamics and gelation phenomena. SIAM Journal on Scientific Computing, 22(3), 802–821. Flint, J.H., Marra, R.A., & Haggerty, J.S. (1986). Powder temperature, size, and number density in laser-driven reactions. Aerosol Science and Technology, 5(2), 249–260. Frenklach, M., Ting, L., Wang, H., & Rabinowitz, M.J. (1996). Silicon particle formation in pyrolysis of silane and disilane. Israel Journal of Chemistry, 36(3), 293–303. Goodson, M., & Kraft, M. (2002). An efficient stochastic algorithm for simulating nano-particle dynamics. Journal of Computational Physics, 183(1), 210–232. Ho, P., Coltrin, M.E., & Breiland, W.G. (1994). Laser-induced fluorescence measurements and kinetic analysis of Si atom formation in a rotating disk chemical vapor deposition reactor. The Journal of Physical Chemistry, 98(40), 10138–10147. Houf, W.G., Grcar, J.F., & Breiland, W.G. (1993). A model for low pressure chemical vapor deposition in a hot-wall tubular reactor. Materials Science and Engineering: B, 17(1–3), 163–171. Kastner, C.A., Braumann, A., Man, P.L.W., Mosbach, S., Brownbridge, G.P.E., Akroyd, J., Kraft, M., & Himawan, C. (2013). Bayesian parameter estimation for a jet-milling model using metropolis-hastings and wang-landau sampling. Chemical Engineering Science, 89, 244–257. Knipping, J., Wiggers, H., Rellinghaus, B., Roth, P., Konjhodzic, D., & Meier, C. (2004). Synthesis of high purity silicon nanoparticles in a low pressure microwave reactor. Journal of Nanoscience and Nanotechnology, 4(8), 1039–1044. Körmer, R., Jank, M.P.M., Ryssel, H., Schmid, H.J., & Peukert, W. (2010a). Aerosol synthesis of silicon nanoparticles with narrow size distribution—Part 1: Experimental investigations. Journal of Aerosol Science, 41(11), 998–1007. Körmer, R., Schmid, H.J., & Peukert, W. (2010b). Aerosol synthesis of silicon nanoparticles with narrow size distribution—Part 2: Theoretical analysis of the formation mechanism. Journal of Aerosol Science, 41(11), 1008–1019. Kraft, M., & Mosbach, S. (2010). The future of computational modelling in reaction engineering. Philosophical Transactions of the Royal Society A, 368(1924), 3633–3644. Lavvas, P., Sander, M., Kraft, M., & Imanaka, H. (2011). Surface chemistry and particle shape: Processes for the evolution of aerosols in titan's atmosphere. The Astrophysical Journal, 728(2), 80. Li, G., Wang, S.-W., & Rabitz, H. (2002). Practical approaches to construct RS-HDMR component functions. The Journal of Physical Chemistry A, 106(37), 8721–8733. Man, P.L.W., Braumann, A., & Kraft, M. (2010). Resolving conflicting parameter estimates in multivariate population balance models. Chemical Engineering Science, 65(13), 4038–4045. Mangolini, L. (2013). Synthesis, properties, and applications of silicon nanocrystals. Journal of Vacuum Science & Technology B, 31(2), 020801. Menz, W.J., & Kraft, M. (2013a). A new model for silicon nanoparticle synthesis. Combustion and Flame, 160, 947–958. Menz, W.J., & Kraft, M. (2013b). The suitability of particle models in capturing aggregate structure and polydispersity. Aerosol Science & Technology, 47, 734–745. Menz, W.J., Akroyd, J., & Kraft, M. (2014). Stochastic solution of population balance equations for reactor networks. Journal of Computational Physics, 256, 615–629. Menz, W.J., Patterson, R.I.A., Wagner, W., & Kraft, M. (2013). Application of stochastic weighted algorithms to a multidimensional silica particle model. Journal of Computational Physics, 248, 221–234. Menz, W.J., Shekar, S., Brownbridge, G., Mosbach, S., Körmer, R., Peukert, W., & Kraft, M. (2012). Synthesis of silicon nanoparticles with a narrow size distribution: A theoretical study. Journal of Aerosol Science, 44, 46–61. Morgan, N.M., Wells, C.G., Goodson, M.J., Kraft, M., & Wagner, W. (2006). A new numerical approach for the simulation of the growth of inorganic nanoparticles. Journal of Computational Physics, 211(2), 638–658. Mosbach, S., Braumann, A., Man, P.L.W., Kastner, C.A., Brownbridge, G.P.E., & Kraft, M. (2012). Iterative improvement of Bayesian parameter estimates for an engine model by means of experimental design. Combustion and Flame, 159(3), 1303–1313. Mosbach, S., Hong, J.H., Brownbridge, G.P.E., Kraft, M., Gudiyella, S., & Brezinsky, K. (2014). Bayesian error propagation for a kinetic model of npropylbenzene oxidation in a shock tube. International Journal of Chemical Kinetics, 46(7), 389–404. Murthy, T., Miyamoto, N., Shimbo, M., & Nishizawa, J. (1976). Gas-phase nucleation during the thermal decomposition of silane in hydrogen. Journal of Crystal Growth, 33(1), 1–7. Nguyen, H.V., & Flagan, R.C. (1991). Particle formation and growth in single-stage aerosol reactors. Langmuir, 7(8), 1807–1814. Onischuk, A.A., Levykin, A.I., Strunin, V.P., Sabelfeld, K.K., & Panfilov, V.N. (2000). Aggregate formation under homogeneous silane thermal decomposition. Journal of Aerosol Science, 31(11), 1263–1281. Patterson, R.I.A., Singh, J., Balthasar, M., Kraft, M., & Norris, J.R. (2006). The linear process deferment algorithm: A new technique for solving population balance equations. SIAM Journal on Scientific Computing, 28(1), 303. Petersen, E.L., & Crofton, M.W. (2003). Measurements of high-temperature silane pyrolysis using SiH4 IR emission and SiH2 laser absorption. The Journal of Physical Chemistry A, 107(50), 10988–10995. Rabitz, H., & Aliş, Ö.F. (1999). General foundations of high-dimensional model representations. Journal of Mathematical Chemistry, 25(2–3), 197–233. Sander, M., West, R.H., Celnik, M.S., & Kraft, M. (2009). A detailed model for the sintering of polydispersed nanoparticle agglomerates. Aerosol Science and Technology, 43(10), 978–989. Shekar, S., Menz, W.J., Smith, A.J., Kraft, M., & Wagner, W. (2012a). On a multivariate population balance model to describe the structure and composition of silica nanoparticles. Computers & Chemical Engineering, 43, 130–147. Shekar, S., Smith, A.J., Menz, W.J., Sander, M., & Kraft, M. (2012b). A multidimensional population balance model to describe the aerosol synthesis of silica nanoparticles. Journal of Aerosol Science, 44, 83–98. Sinniah, K., Sherman, M.G., Lewis, L.B., Weinberg, W.H., Yates, J.T., Jr, & Janda, K.C. (1990). Hydrogen desorption from the monohydride phase on Si 100. The Journal of Chemical Physics, 92, 5700. Swihart, M.T., & Girshick, S.L. (1999). Thermochemistry and kinetics of silicon hydride cluster formation during thermal decomposition of silane. The Journal of Physical Chemistry B, 103(1), 64–76. Tomlin, A.S. (2013). The role of sensitivity and uncertainty analysis in combustion modelling. Proceedings of the Combustion Institute, 34(1), 159–176. Vikhansky, A., & Kraft, M. (2004). A Monte Carlo methods for identification and sensitivity analysis of coagulation processes. Journal of Computational Physics, 200(1), 50–59. Vikhansky, A., & Kraft, M. (2006). Two methods for sensitivity analysis of coagulation processes in population balances by a Monte Carlo method. Chemical Engineering Science, 61(15), 4966–4972. Wu, J.J., Nguyen, H.V., & Flagan, R.C. (1987). A method for the synthesis of submicron particles. Langmuir, 3(2), 266–271. Ziehn, T., & Tomlin, A.S. (2008). Global sensitivity analysis of a 3D street canyon model–Part I: The development of high dimensional model representations. Atmospheric Environment, 42(8), 1857–1873.