Stochastic numerics for the gas-phase synthesis of nanoparticles - - PowerPoint PPT Presentation

Stochastic numerics for the gas-phase synthesis of nanoparticles

Shraddha Shekar1, Alastair J. Smith1, Markus Sander1, Markus Kraft1 and Wolfgang Wagner2

1Department of Chemical Engineering and Biotechnology

University of Cambridge

2Weierstrass Institute for Applied Analysis and Stochastics, Berlin

March 2011

CoMo

GROUP

March 2011 1 / 45

Outline

1 Introduction

Motivation

2 Model

Type space Particle processes Algorithm

3 Numerical study 4 Conclusion

CoMo

GROUP

March 2011 2 / 45

Motivation

What are nanoparticles? Why are they important?

Particles sized between 1-100 nm. Both inorganic and organic nanoparticles find wide applications in various fields.

Why model nanoparticle systems?

To optimise industrial operations and to obtain products of highly specific properties for sensitive applications. To understand the molecular level properties that are difficult to be

• bserved experimentally.

CoMo

GROUP

March 2011 3 / 45

Motivation II

Salient features of the current model: Fully-coupled multidimensional stochastic population balance model. Describing various properties of nanoparticles at an unprecedented level of detail. Tracking properties not only at macroscopic level but also at a molecular level.

CoMo

GROUP

March 2011 4 / 45

Applications of silica nanoparticles

Silica nanoparticles are amorphous and have Si:O = 1:2. Their applications include: Catalysis Bio-medical applications Support material for functional nanoparticles Fillers/Binders Optics

CoMo

GROUP

March 2011 5 / 45

Physical system

To describe the system at a macroscopic level it is essential to understand it at a molecular level.

Precursor (TEOS) Flame reactor P ≥ 1 atm T ≈ 1100 - 1500 K Silica nanoparticles Industrial Scale Molecular Scale

CoMo

GROUP

March 2011 6 / 45

Type space I

Each particle is represented as: Pq = Pq(p1, . . . , pn(Pq), C). Particle Pq consists of n(Pq) primary particles pi with i ∈ {1, . . . , n(Pq)} and q ∈ {1, . . . , N}, where N is the total number

• f particles in the system.

Si O O O O Si O OH O Si Si Si

Pq = Pq(p1,...,pn(Pq),C)

O Si O O

pi = pi(ηSi,ηO,ηOH)

HO OH HO

Figure: Type Space.

CoMo

GROUP

March 2011 7 / 45

Type space II

Each primary particle pi is represented as: pi = pi(ηSi, ηO, ηOH) where ηx (ηx ∈ Z with ηx ≥ 0) is the number of chemical units of type x ∈ {Si, O, OH}.

CoMo

GROUP

March 2011 8 / 45

Type space III

C is a lower diagonal matrix of dimension n(Pq) × n(Pq) storing the common surface between two primary particles: C(Pq) =          · · · · · · C21 ... · · · . . . ... ... · · · . . . Ci1 · · · Cij ... . . . . . . · · · . . . · · · . . .          . The element Cij of matrix C has the following property: Cij =

• 0, if pi and pj are non-neighbouring ,

Cij > 0, if pi and pj are neighbouring.

CoMo

GROUP

March 2011 9 / 45

Particle processes

Particles are transformed by the following processes: Inception Surface reaction Coagulation Sintering Intra-particle reaction

CoMo

GROUP

March 2011 10 / 45

Particle processes

Particles are transformed by the following processes: Inception Surface reaction Coagulation Sintering Intra-particle reaction

CoMo

GROUP

March 2011 11 / 45

Inception

Two molecules in gas phase collide to form a particle consisting of one primary.

HO Si OH OH OH HO Si OH OH HO

+

HO Si OH OH OH HO Si OH OH HO HO Si OH O OH Si OH OH HO

• 2 H2O

[monomers] [primary particle]

Figure: Inception of primary particles from gas-phase monomers.

An inception event increases the number of particles in the system molecule + molecule → PN(p1, C), C = 0. Initial state of primary p1 given by: p1 = p1(ηSi = 2, ηO = 1, ηOH = 6),

CoMo

GROUP

March 2011 12 / 45

Inception rate

Inception rate for each particle (Pq) calculated using the free molecular kernel: Rinc(Pq) = 1 2K fmNA2C 2

g ,

NA is Avogadro’s constant, Cg is the gas-phase concentration of the incepting species (Si(OH)4), K fm = 4

• πkBT

mg (d2

g ),

kB is the Boltzmann constant, T is the system temperature, mg and dg are the mass and diameter respectively of the gas-phase molecule Si(OH)4.

CoMo

GROUP

March 2011 13 / 45

Particle processes

Particles are transformed by the following processes: Inception Surface reaction Coagulation Sintering Intra-particle reaction

CoMo

GROUP

March 2011 14 / 45

Surface reaction

Dehydration reaction between gas-phase monomer and particle surface:

O Si O O Si HO O

+Si(OH)4

• H2O

O Si O O Si O O Si OH OH HO

Figure: Surface reaction between a particle and a gas-phase molecule.

Surface reaction transforms particle as: Pq + molecule → Pq(p1, ., pi ′, .., pn(Pq), C′), p′

i → pi(ηSi + 1, ηO + 1, ηOH + 2).

CoMo

GROUP

March 2011 15 / 45

Particle rounding due to surface reaction

Surface reaction also alters the common surface (C → C′). Net common surface area of pi changes due to volume addition: ∆s(pi) = (v(pi ′) − v(pi)) 2σ dp(pi), where σ is the surface smoothing factor (0 ≤ σ ≤ 2). C′ is given by: C ′

ij =

• 0, if pi and pj are non-neighbouring ,

Cij + ∆s(pi), if pi and pj are neighbouring.

CoMo

GROUP

March 2011 16 / 45

Surface reaction rate

Surface reaction rate calculated using equation of Arrhenius form: Rsurf(Pq) = Asurf exp

• − Ea

RT

• ηOH(Pq)NACg,

Asurf is pre-exponential factor (obtained from collision theory), Ea is activation energy, ηOH(Pq) is the total number of –OH sites on particle Pq.

CoMo

GROUP

March 2011 17 / 45

Particle processes

Particles are transformed by the following processes: Inception Surface reaction Coagulation Sintering Intra-particle reaction

CoMo

GROUP

March 2011 18 / 45

Coagulation

Two particles collide and stick to each other:

+

Pq Pr Ps

Figure: Coagulation between two particles.

Coagulation of particles Pq and Pr forms new particle Ps as: Pq + Pr → Ps(p1, ..., pn(Pq), p(n(Pq)+1), ..., pn(Pq)+n(Pr), C).

CoMo

GROUP

March 2011 19 / 45

Coagulation II

Primary pi from particle Pq and primary pj from Pr are assumed to be in point contact. The matrix C(Ps) is calculated as: C(Ps) =          . . . C(Pq) · · · Cij · · · . . . . . . . . . Cji . . . C(Pr) . . .          and has dimension n(Ps) × n(Ps), where n(Ps) = n(Pq) + n(Pr).

CoMo

GROUP

March 2011 20 / 45

Coagulation rate

Coagulation rate between Pq and Pr calculated using transition coagulation kernel: K tr(Pq, Pr) = K sf(Pq, Pr)K fm(Pq, Pr) K sf(Pq, Pr) + K fm(Pq, Pr), where the slip-flow kernel is: K sf(Pq, Pr) = 2kBT 3µ 1 + 1.257Kn(Pq) dc(Pq) + 1 + 1.257Kn(Pr) dc(Pr)

• × (dc(Pq) + dc(Pr)) ,

and the free molecular collision kernel is: K fm(Pq, Pr) = 2.2

• πkBT

2

• 1

m(Pq) + 1 m(Pr) 1

2

(dc(Pq) + dc(Pr))2

CoMo

GROUP

March 2011 21 / 45

Particle processes

Particles are transformed by the following processes: Inception Surface reaction Coagulation Sintering Intra-particle reaction

CoMo

GROUP

March 2011 22 / 45

Sintering

Sintering described using viscous-flow model.

pi pj pi pj pk

No Sintering Partial Sintering Complete Sintering

Figure: Evolution of sintering process with time.

Sintering between pi and pj of a single particle Pq calculated on a primary particle-level.

CoMo

GROUP

March 2011 23 / 45

Sintering level

Sintering level defined to represent degree of sintering between pi and pj: s(pi, pj) =

Ssph(pi,pj) Cij

− 2

1 3

1 − 2

1 3

. Ssph(pi, pj) is the surface area of a sphere with the same volume as the two primaries. Pq conditionally changes depending on the sintering level s(pi, pj). Two types are defined depending upon a threshold (95%): Partial sintering s(pi, pj) < 0.95 Complete sintering s(pi, pj) ≥ 0.95

CoMo

GROUP

March 2011 24 / 45

Partial sintering

Surface areas of primaries are reduced by a finite amount. –OH sites at contact surface react to form Si–O–Si bonds:

Si O O O OH Si Si Si Si O O O HO Si Si Si OH HO

• 2H2O

OH OH OH OH OH HO OH OH OH OH OH OH OH OH OH OH OH OH OH OH OH OH OH Si O O O O Si Si Si Si O O O Si Si Si O OH OH OH OH OH HO OH OH OH OH OH OH OH OH OH OH OH OH OH OH OH OH OH

pi pj pi pj

Reactions at particle neck

Figure: Dehydration reaction due to sintering.

CoMo

GROUP

March 2011 25 / 45

Partial sintering II

Surface density of –OH sites assumed constant throughout sintering. The change in the internal variables of primaries pi and pj given by: ∆ηOH(pi) = ∆ηOH(pj) = ρs(Pq)∆Cij/2, ∆ηO(pi) = ∆ηO(pj) = −0.5 × ∆ηOH(pi), ∆ηSi(pi) = ∆ηSi(pj) = 0.

CoMo

GROUP

March 2011 26 / 45

Partial sintering III

Particle continuously transforms due to partial sintering as: Pq(p1, . . . , pn(Pq), C) → Pq(p1, .., p′

i, p′ j, .., pn(Pq), C′)

where, p′

i = pi(ηSi, ηO − ∆ηO(pi), ηOH − ∆ηOH(pi)),

p′

j = pj(ηSi, ηO − ∆ηO(pj), ηOH − ∆ηOH(pj)).

Element of C′ given by: C ′

ij = Cij −

∆t τ(pi, pj) (Cij − Ssph(pi, pj)) , where ∆t is a time interval.

CoMo

GROUP

March 2011 27 / 45

Complete sintering

pi and pj replaced by new primary p′′

k.

Particle transforms due to complete sintering as: Pq(p1, . . . , pn(Pq), C) → Pq(p1, .., p′′

k, .., pn(Pq), C′′),

where new primary: p′′

k = p′′ k (ηSi(pi) + ηSi(pj), ηO(pi) + ηO(pj), ηOH(pi) + ηOH(pj)) .

CoMo

GROUP

March 2011 28 / 45

Complete sintering II

C is changed by removing columns and rows i and j and adding new column and row k: C′′ =         · · · · · · . . . ... . . . . . . . . . C ′′

k1

· · · · · · . . . . . . . . . ... . . . C(n(Pq)−1)1 · · · C ′′

(n(Pq)−1)k

· · ·         .

CoMo

GROUP

March 2011 29 / 45

Sintering rate

Sintering rate between pi and pj equivalent to rate of change of their common surface ∆Cij in time ∆t: ∆Cij ∆t = − 1 τ(pi, pj)(Cij − Ssph(pi, pj)), Ssph(pi, pj) is the surface area of a sphere with the same volume as the two primaries.

CoMo

GROUP

March 2011 30 / 45

Characteristic sintering time

Characteristic sintering time of pi and pj is: τ(pi, pj) = As × dp(pi, pj) × exp Es T

• 1 −

dp,crit dp(pi, pj)

• ,

where dp(pi, pj) is the minimum diameter of pi and pj, and As, Es and dp,crit are sintering parameters.

CoMo

GROUP

March 2011 31 / 45

Particle processes

Particles are transformed by the following processes: Inception Surface reaction Coagulation Sintering Intra-particle reaction

CoMo

GROUP

March 2011 32 / 45

Intra-particle reaction

Reaction of two adjacent –OH sites on one particle:

O Si OH O O Si HO O O Si O Si O O O Si O O Si

• H2O

Figure: Intra-particle reaction.

Intra-particle reaction transforms particle as: Pq(p1, .., pi, .., pn(Pq), C) → Pq(p1, .., p′

i, .., pn(Pq), C),

where p′

i → pi(ηSi, ηO + 1, ηOH − 2).

CoMo

GROUP

March 2011 33 / 45

Intra-particle reaction rate

Deduce rate of intra-particle reaction from surface reaction rate and avergae sintering rate such that Si/O ratio of 1/2 is attained: Rint(Pq) = Asurf exp

• − Ea

RT

• ηOH(Pq)NACg

−ρs(Pq) 2  

n(Pq)

• i,j=1

Cij − Ssph(pi, pj) τ(pi, pj)   , where ρs(Pq) = ηOH(Pq)/S(Pq) is the surface density of active sites.

CoMo

GROUP

March 2011 34 / 45

Algorithm

Input: State of the system Q0 at initial time t0 and final time tf. Output: State of the system Q at final time tf. t ← − t0,Q ← − Q0; while t < tf do Calculate an exponentially distributed waiting time τ with parameter; Rtot(Q) = Rinc(Q) + Rcoag(Q) + Rsurf(Q) + Rint(Q). Choose a process m according to the probability P(m) = Rm(Q) Rtot(Q), where Rm is the rate of the process m ∈ {inc, coag, surf, int}; Perform process m; Update sintering level of all particles; Increment t ← − t + τ; end

CoMo

GROUP

March 2011 35 / 45

Numerical study

For a given property of the system ξ calculated using Nsp computational particles and L number of independent runs, the empirical mean is: µ(Nsp,L)

1

(t) = 1 L

L

• l=1

ξ(Nsp,l)(t). The empirical variance is: µ(Nsp,L)

2

(t) = 1 L

L

• l=1

ξ(Nsp,l)(t)2 − µ(Nsp,L)

1

(t)

2

.

CoMo

GROUP

March 2011 36 / 45

Confidence interval

The confidence interval IP within which there is a probability P of finding the true solution is then given by: IP =

• µ(Nsp,L)

1

(t) − cP, µ(Nsp,L)

1

(t) + cP

• .

cP = aP

• µ(Nsp,L)

2

(t) L . We use aP = 3.29 which corresponds to P = 0.999(99.9%)

CoMo

GROUP

March 2011 37 / 45

Error

The error e is: e(Nsp,L)(t) =

• µ(Nsp,L)

1

(t) − ζ∞(t)

• ,

ζ∞(t) is an approximation for the true solution which is obtained from a "high-precision calculation" with a very large number of particles. The average error over the entire simulation time is: ¯ e(Nsp, L) = 1 M

M

• j=1

e(Nsp,L)(tj), where the M time steps tj are equidistant.

CoMo

GROUP

March 2011 38 / 45

M0

The zeroth moment is the particle number density: M0(t) = N(t) Vsmpl

1.5 2 2.5 3 3.5 4 4.5 5 7.5 8 8.5 9 9.5

log(Nsp) log(¯ e) (cm−3)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 10

8

10

9

10

10

10

11

10

12

Time (s) M0 (cm−3)

M0 Nsp=64, L=256 Nsp=512, L=128 Nsp=16384, L=4 High Precision Solution

Figure: Convergence of zeroth moment (Nsp × L = 65536). Solid line indicates slope of -1.

.

CoMo

GROUP

March 2011 39 / 45

Volume

The average particle volume: V (t) = 1 N(t)ΣN(t)

q=1 V (Pq(t))

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1 2 3 4 5 6 7 x 10

−9

Time (s) Volume (cm3)

1.5 2 2.5 3 3.5 4 4.5 5 −18.8 −18.6 −18.4 −18.2 −18 −17.8 −17.6

log(Nsp) log(¯ e) (cm3)

Average Volume Nsp=64, L=256 Nsp=512, L=128 Nsp=16384, L=4 High Precision Solution

Figure: Convergence of average volume (Nsp × L = 65536). Solid line indicates slope of -1.

.

CoMo

GROUP

March 2011 40 / 45

Collision diameter

The average collision diameter of a particle: Dc(t) = 1 N(t)ΣN(t)

q=1 dc(Pq(t))

0.2 0.4 0.6 0.8 10 20 30 40 50 60 70 80

Time (s) Collision Diameter (nm)

1.5 2 2.5 3 3.5 4 4.5 5 −0.4 −0.2 0.2 0.4 0.6 0.8 1

log(Nsp) log(¯ e) (nm)

Average Collision Diameter Nsp=64, L=256 Nsp=512, L=128 Nsp=16384, L=4 High Precision Solution

Figure: Convergence of average average collision diameter (Nsp × L = 65536). Solid line indicates slope of -1.

.

CoMo

GROUP

March 2011 41 / 45

Computational time

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 10 20 30 40 50 60 70 80

Time (s) Collision Diameter (nm)

10

1

10

2

10

3

10

4

10

5

10 10

1

10

2

10

3

Nsp Computational Time (s)

Nsp=64, L=256 Nsp=512, L=128 Nsp=163824, L=4 High Precision Solution Nsp=64 Nsp=512 Nsp=163824 High Precision Solution

Figure: Computational time for different Nsp.

CoMo

GROUP

March 2011 42 / 45

Computer-generated TEM

Figure: TEM images generated by projecting particles onto a plane. Experimental values from Seto et al. (1995).

CoMo

GROUP

March 2011 43 / 45

Conclusion

Description of a detailed population balance model. Numerical studies performed. Demonstrated feasibility of using first-principles to model complex nanoparticle synthesis processes.

CoMo

GROUP

March 2011 44 / 45

Thank You!

CoMo

GROUP

March 2011 45 / 45