Dual Population Balance Monte Carlo Simulation of Particle Synthesis by Flame Spray Pyrolysis

Abstract

The Dual Population Balance Monte Carlo Method (DPBMC) takes into account the full size spectrum of the droplet and particle phase. Droplet and particle size distributions are rendered by weighted simulation particles. This allows for an accurate description of particle nucleation and coagulation and droplet combustion, simultaneously. Internal droplet properties such as temperature and concentrations fields are used to define criteria for the onset of droplet breakage in the framework of weighted Monte Carlo droplets. We discuss the importance of droplet polydispersity on particle formation in metal oxide particle synthesis, which is shown to strongly affect particle formation and growth. The method is applied to particle synthesis from metal nitrate precursor solutions with flame spray pyrolysis (FSP) and compared to experiments from literature.

1. Introduction

Flame spray pyrolysis (FSP) is widely used for synthesis of oxide particles. It offers in contrast to conventional flame synthesis more possibilities regarding material choice and design [1]. The main reason is the availability of liquid precursors for a broad spectrum of materials. In addition, the synthesis of composite materials is possible by mixing different precursor solutions. In recent years, the research in this field gained additional traction motivated by economic interests. At the laboratory scale, a production rate of over 100 g/h of nanoparticles can be achieved [2]. This makes the method very promising for large scale production. However, liquid precursors pose their own difficulties. Often, they can not be vaporized and have to be atomized before combustion. The addition of a liquid phase besides the gas phase and particle phase causes increased complexity in both experiment and modeling. Therefore, the FSP was largely regarded as a ‘black box’. Different parameters were correlated with material properties such as particle size [3], shell thickness [4] or composition [5]. A lack of in situ measurement techniques amplifies the difficulties in understanding the droplet-to-particle conversion process. Burning droplets may destroy transmission electron microscopy grids. This makes sampling of particles in the presence of droplets impractical. There is also no technique to measure droplet composition inside the flame. In conclusion, experimental information about the early particle formation stage is very difficult to obtain. The combination of a large variety of possible experimental parameters and lack of measurement techniques requires the development and experimental validation of dedicated models.

Population balance (PB) approaches showed promising results in predicting particle properties in FSP even though relatively simple models were used. A monodisperse droplet size distribution (DSD) was applied to predict the size evolution of particles in spray synthesis. Deviations between particle size measurements and modeling were usually attributed to the simplifications used (e.g., instant droplet evaporation or monodisperse primary and agglomerate particles [6]). Recent suggestions to account for polydispersity of the droplet phase impose several limiting conditions. Heine and Pratsinis [3] superimpose monodisperse droplet calculations to evaluate the influence of a polydisperse DSD on the resulting particle size distribution (PSD). They show that large droplets tend to produce fewer but larger particles. Widiyastuti et al. [4] account for droplet polydispersity by applying 18 size classes for droplets, but assume the total number of droplets to be constant during the whole synthesis process. They found their method to be accurate in predicting particle morphology. However, in a flame spray synthesis reactor, regions with different droplet concentrations and mean diameters exist. According to the $d^2$—law [7], small droplets have a shorter lifetime than large droplets. Hence, precursor release and particle formation along the reactor axis depend on the size and number of droplets present. Early particle formation may occur when there is still a considerable amount of droplets present. Therefore, the aim is to describe precursor release and early particle formation in the presence of burning droplets by applying a PB approach both to the particle and droplet phase.

The Monte Carlo method is able to capture discrete events very well due to its own discrete nature. This includes removal of fully evaporated droplets and droplet breakage. The latter promotes gas-phase particle nucleation and formation of smaller particles as reported by Rosebrock et al. [8]. They also observed cascade-like burning and break up of droplet fragments from the original droplet. Li et al. [9] as well as Rosebrock et al. [8] applied a multicomponent diffusion model to single droplet combustion for nanoparticle production and compared the results to experimental data. They found that temperature and concentration fields inside the droplet undergo changes as the volatile component evaporates and a viscous shell is formed.

The model employed in this work aims to be as computationally efficient as possible for computation of 100–1000 Monte Carlo droplets (MC-droplets) over up to 100 single MC-simulations, while still giving insight into whether the gas-to-particle or droplet-to-particle mechanism is preferred for the respective MC-droplet. In this way, the full size spectrum of the DSD can be taken into account. There are several methods available to simulate droplet and particle breakage in the framework of the Monte Carlo method if the appropriate breakage function and breakage rate are supplied [10,11]. For these kinds of problems, breakage functions and breakage rates are currently unknown. As a first step, a mechanism is suggested to predict a droplet fragmentation event.

Droplets with different sizes experience different velocities relative to the gas phase. Some droplets may meet the conditions for liquid-phase nucleation of particles, while others do not. The application of weighted particles facilitates following these features of individual droplets. Particles on the other side experience temperature and supersaturation fluctuations inside the flame. This leads to a rapidly varying Kelvin diameter and therefore alternating condensation and evaporation of particles, which is a challenge for a number of numerical schemes. The Monte Carlo method can simulate this behaviour well [12,13]. The drawback is a high computational cost of the Monte Carlo method, which can be reduced by parallel calculations on GPUs [14]. The model presented in this work can be implemented for parallel calculations on GPUs where each parallel thread processes all calculations necessary for a given particle or droplet, which helps in reducing the computation time.

2. Population Balance Modeling

The aim of this work is to describe the effect of droplet combustion kinetics on particle growth. Spray droplets are assumed to be evenly distributed in the simulation domain. Since the number concentration of droplets is small, interaction and coagulation between droplets are neglected. The vapour pressure of ${\mathrm{Fe}}_2{\mathrm{O}}_3$ at flame temperatures used here is nearly zero; it can be assumed that $v^{*}$ = $v_1$ and therefore evaporation is negligible. Consequently, single molecules are treated as the first step towards particles. This approach is described by Girshick et al. [15] for plasma synthesis of oxide particles. Particles grow due to Brownian coagulation in the free molecular regime. The general form of PBEs for particles and droplets is given in Equations (1) and (2):

$$\frac{\mathrm{d}{n}_{\mathrm{p}}(v)}{\mathrm{d}t}=\underset{0}{\overset{v}{\int }}\beta (v-v^{\prime },v^{\prime })n(v-v^{\prime })n(v^{\prime })\mathrm{d}{v}^{\prime }-n(v)\underset{0}{\overset{\infty }{\int }}\beta (v,v^{\prime })n(v^{\prime })\mathrm{d}{v}^{\prime }+J(t)\delta (v^{*}-v)$$

(1)

$$\frac{\mathrm{d}{n}_D(v)}{\mathrm{d}t}=-\frac{\partial [E(v)n_D(v)]}{\partial v}-\delta (v_1-v){\left.\frac{\partial [E(v)n_D(v)]}{\partial v}\right|}_{v=v_1}.$$

(2)

The Dirac-Delta function $\delta$ in Equations (1) and (2) describes discrete events for MC-particles and MC-droplets, whereas E gives the evaporation rate of the individual droplet. MC-particles are added in the case of nucleation and are removed in the case of complete evaporation or fragmentation. For this work, the gas-to-particle route has been used for the MC-Simulations as sketched in Figure 1. For nucleation and coagulation of the particle phase, a constant number Monte Carlo scheme using the concepts of ‘stochastic resolution’ and ‘low-weight-merging’ [16] was applied. One simulation yields then the average of 100 subsimulations. The particle PB model used here was validated using a discrete sectional model [16].

Additionally, instant precursor decomposition in the gas-phase is assumed. This equates then the monomer production rate from droplets to nanoparticle nucleation rate. The same approach is also used by Ulrich [17]. Since two Fe-atoms are needed to produce one ${\mathrm{Fe}}_2{\mathrm{O}}_3$-monomer, the nucleation rate is half of the precursor release rate. The nucleation rate couples the droplet to the particle phase according to Equation (3). After each timestep $\Delta t$, a MC-particle is created with weight ${\omega }_P$ equaling the monomer production rate $m_1$ in $\Delta t$ and having monomer volume $v_1$ and the nucleation rate is:

$${\mathrm{J}}^{*}=\frac{1}{2}\sum _{i=1}^N{m}_{1,i}{w}_{D,i}=\frac{{\omega }_P}{\Delta t}.$$

(3)

A population balance approach to the droplet phase is applied to improve the calculation of nucleated particles by including the full size spectrum of the DSD. In contrast to the particle phase, a non-constant number approach for the droplet phase is used. If a droplet evaporates below a preset threshold value, the respective MC-droplet is removed from the simulation. This is realistic since, during synthesis, the droplet concentration may decrease to due to combustion. Although a non-constant number approach results in increased statistical noise for the droplet phase, it prevents the simulation of temperature and concentration fields of low-weighted droplets and thus significantly reduces computation time. The influence of the number of initial MC-droplets on the Sauter mean diameter (SMD) and ${\mathrm{J}}^{*}$ is shown in Figure 2. As the number of MC-droplets decreases, the noise in SMD increases. This noise does not occur in ${\mathrm{J}}^{*}$ because droplets release monomers steadily as they evaporate, whereas the removal of a MC-droplet with the weight ${\omega }_{D,i}$ is a discrete event.

For more than 100 MC-droplets per single Monte Carlo simulation, virtually no changes are visible in the nucleation rate ${\mathrm{J}}^{*}$. This is thus the optimal number for maximum accuracy of ${\mathrm{J}}^{*}$ and consequently the particle phase with minimum computation time for the droplet phase. For discussion of DSD properties, 1000 MC-droplets will be used to minimize the evaporation noise.To reduce the computation time for the droplet phase further, a parallel scheme for droplet size change and monomer release is applied as illustrated in Figure 3.

As a first step, the size change and number of monomers released by every droplet is calculated by one GPU-thread each. Then, the Thrust library is used to perform parallel reduction on released monomers of each droplet. From the sum, the nucleation rate of particles can be calculated for every Monte Carlo simulation. The advantage of parallel computation on GPUs is that computation time does not increase proportionally to number of elements computed. Figure 4 shows the computation time on a Tesla P100 GPU associated with the number of initial MC-droplets per single Monte Carlo simulation. The comparison of the computation times for 100 and 2000 MC-droplets shows a doubling of computation time for a twenty times increase in computed elements. Thus, the advantage of parallelization compared to serial computing is a below proportional increase in computation time as the number of computed elements is increased.

Figure 4 also shows the effect of GPU memory bandwith limitation on the computation time. For 10,000 MC-droplets, a high increase in computation time occurs. Because more data from the main memory is requested than can be computed in parallel, additional data is then processed in serial. Consequently, the number of initial MC-droplets can be varied to achieve the desired accuracy for the droplet phase and computation time depending on the GPU hardware available.

3. Droplet Combustion

Rosebrock et al. [18] measured temperature gradients inside burning droplets using Rainbow Refractometry and compared them to different combustion models. They found that diffusion limited heat transport models to give good agreement with the measured data. One model of this type was developed by Law [19,20] and Sirignano [21] and is also applied in this work. Isobaric and spherically symmetric combustion of droplets in the presence of an oxidizer gas and an inert species is assumed. The gas-phase processes are assumed to be quasi-steady. Additionally, the flame sheet approximation is used, which approximates the flame front by an infinitesimally thin layer where stoichiometric combustion occurs. In case of ethanol combustion in hot air, this results in the following reaction mechanism (Equation (4)):

$${\mathrm{C}}_2{\mathrm{H}}_6\mathrm{O}+11.29·{\mathrm{N}}_2+3·{\mathrm{O}}_2\to 2·{\mathrm{CO}}_2+11.29·{\mathrm{N}}_2+3·{\mathrm{H}}_2\mathrm{O}$$

(4)

From Equation (4), the stoichiometric oxygen to fuel mass ratio $\nu =2.08$ and ambient mass fraction of water $Y_{H_{2}O,\infty }=0.1047$ can be calculated.

3.1. Heat and Mass Transport

To obtain the droplet evaporation rate $\widehat{m}$, the state of the droplet surface regarding temperature and composition at liquid and gas side has to be known. The heat and mass diffusion equations are given in Equations (5) and (6):

$$\frac{\partial T}{\partial t}=D_H\left(\frac{{\partial }^{2}T}{\partial r^2}+\frac{2}{r}\frac{\partial T}{\partial r}\right),$$

(5)

$$\frac{\partial Y_i}{\partial t}=D_M\left(\frac{{\partial }^2{Y}_i}{\partial r^2}+\frac{2}{r}\frac{\partial Y_i}{\partial r}\right).$$

(6)

Heat diffusion is solved for the complete liquid-phase while mass diffusion is solved for ethanol $Y_{C_2{H}_{6}O}$ and water $Y_{H_{2}O}$. Since the sums of mass fractions $Y_i$ and mole fractions $X_i$ equal 1, the mass fraction of the precursor $Y_{Prec}$ is then known. Their respective boundary conditions are given by Equations (7) and (8) as summarized by Sirignano [21].

$${\left.\frac{\partial T}{\partial t}\right|}_{r=0}=0;{\left.\frac{\partial T}{\partial t}\right|}_R=\frac{ln(1+B)}{R}\frac{{\lambda }_g}{{\lambda }_l}\left[\frac{T_{\infty }-T_R+\nu Q{Y}_{O,\infty }/c_p}{B}-\frac{L}{c_p}\right],$$

(7)

$${\left.\frac{\partial Y_i}{\partial t}\right|}_{r=0}=0;{\left.\frac{\partial Y_i}{\partial t}\right|}_R=\frac{{\rho }_g{D}_{M,g}}{{\rho }_l{D}_{M,l}}{\displaystyle \frac{ln(1+B)}{R}}(Y_{i,R}-{\epsilon }_i).$$

(8)

If the liquid side mass fractions and the surface temperature $T_R$ are known, Raoults Law and the Clausius–Clapeyron equation can be used to calculate the gas side mass and mole fractions of each species (Equation (9)):

$$X_{i,R^{+}}=X_{i,R^{-}}{\displaystyle \frac{1}{p_{\infty }}}exp\left[{\displaystyle \frac{L_i}{R{g}_i}}\left({\displaystyle \frac{1}{T_{b,i}}}-{\displaystyle \frac{1}{T_R}}\right)\right].$$

(9)

The fractional vaporization rates ${\epsilon }_i$ can be calculated from the surface gas side composition with Equation (10):

$${\epsilon }_i=Y_{i,R^{+}}+(1-Y_{F,R^{+}}){\displaystyle \frac{Y_{i,R^{+}}-Y_{i,f}}{Y_{F,R^{+}}-Y_{F,f}}},$$

(10)

$$B_i={\displaystyle \frac{Y_{i,R}-Y_{i,\infty }}{{\epsilon }_i-Y_{i,R}}}.$$

(11)

The Spalding mass transfer number $B_i$ given by Equation (11) is needed for the boundary conditions in Equations (7) and (8). Since water is the only non-combustible species in ambient atmosphere, it can be accounted for by setting $Y_{H_{2}O,f}$ and $Y_{F,f}$ to $Y_{H_{2}O,\infty }$. Since metal nitrates decompose before vaporization, they do not have a relevant vapor pressure. To account for this, Equations (10) and (11) are only applied if $T_R$ is equal to or is higher than the boiling temperature $T_{b,i}$ of the respective precursor species. In this case, $X_{i,R^{+}}$ is set to $X_{i,R^{-}}$. This models instant evaporation of all precursor molecules on the surface. Now, only the temperature function $\widehat{H}$ and $\widehat{m}$ are needed to obtain the droplet regression rate $\alpha$, given by Equations (12) and (13):

$$\widehat{H}=\frac{{\widehat{T}}_{\infty }-{\widehat{T}}_R+\tilde{L}-\tilde{Q}+(1+{\displaystyle \frac{Y_{O,\infty }}{\nu }})[\tilde{Q}-\tilde{L}(1-Y_{F,f})/(1-Y_{F,R^{+}})]}{(1+{\displaystyle \frac{Y_{O,\infty }}{\nu }})(1-Y_{F,f})/(1-Y_{F,R^{+}})-1},$$

(12)

$$\widehat{m}={\displaystyle \frac{{\epsilon }_i}{m_i}}=ln\left(1+\frac{{\widehat{T}}_{\infty }-{\widehat{T}}_R+{\displaystyle \frac{Y_{O,\infty }}{\nu }}\tilde{Q}}{\tilde{L}+\widehat{H}}\right).$$

(13)

To take the relative velocity of the droplet into account, an empirical correction to $\widehat{m}$ according to Ranz and Marshall [22] is applied. The corrected vaporization rate ${\widehat{m}}_c$ is then given by Equation (14):

$${\widehat{m}}_c=\widehat{m}\left[1+0.3P{r}^{1/3}{(2Re)}^{1/2}\right].$$

(14)

For each individual MC-droplet, $\alpha$ (Equation (15)) is calculated to obtain the size evolution.

$${\displaystyle \frac{\mathrm{d}{R}^2}{\mathrm{d}t}}=-{\widehat{m}}_c{\displaystyle \frac{2{\lambda }_g}{{\rho }_l{c}_p}}=\alpha .$$

(15)

Equations (5)–(8) together with Equation (15) describe a moving boundary problem. Here, it is solved with a forward-time and centered-space finite differences (FTCS-FD) approach in one dimension assuming radial symmetry. To preserve the spatial resolution of temperature and concentration fields, the spatial grid is reconstructed after every timestep with a fixed number of gridpoints. For FTCS-FD, the error $ϵ$ is $O({ϵ}^2)$ [23]. In order to validate the FD method, the number of gridpoints is increased until no changes in the nucleation rate are observable. For more than 200 gridpoints per droplet, essentially no changes in the nucleation rate occur as shown in Figure 5.

3.2. Average Liquid and Gas-Phase Properties

A quasi-steady state model requires an averaging of liquid and gas-phase properties at a reference temperature $T_{Ref}$. According to Law and Williams [24], $T_{Ref}$ can be reasonably approximated by the arithmetic mean of flame temperature $T_f$ and boiling temperature $T_b$. $T_{Ref}$ is thus given by Equation (16).

$$T_{Ref,g}=1/2(T_f-T_b)$$

(16)

Following Abramzon and Sirignano [25], the same can be done for the liquid-phase given by Equation (17), where $T_0$ is the initial droplet temperature:

$$T_{Ref,l}=1/2(T_0+T_b).$$

(17)

Material properties are then set to $T_{Ref}$. The boiling point of the mixture and heat capacity of the gas phase are calculated by weighting the mole fractions $X_i$ of all N existing species in the respective phase (Equations (18) and (19)):

$$T_{b,Mix}=\sum _i^N{X}_{i,l}{T}_{b,i},$$

(18)

$$c_{p,\infty }=\sum _i^N{X}_{i,g}{c}_{p,i}.$$

(19)

3.3. Droplet Events

For MC-droplets, the onsets of superheating, liquid-phase precipitation and precursor decomposition are used to define conditions for droplet fragmentation. For this, the following conditions are evaluated for each droplet after every timestep and for every gridpoint:

If the temperature $T(r,t)$ at r exceeds the mixture boiling temperature $T_{b,Mix}(r,t)$ at r, this then marks onset of internal boiling at r:

$$T(r,t)\ge T_{b,Mix}(r,t).$$

(20)

If the temperature $T(r,t)$ at r exceeds the decomposition temperature $T_{Dec}$ of the precursor at r, solid phase formation is assumed:

$$T(r,t)\ge T_{Dec}.$$

(21)

If the solubility limit is reached at $S(r,t)>1$, liquid phase precipitation of the precursor is assumed at r:

$$S(r,t)>1.$$

(22)

For ethanol, the boiling temperature used is 351 K. For metal nitrate precursors, their decomposition temperatures are approximated by their boiling temperatures, which then gives $T_{Dec}=398$ K for ${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$ and ${\mathrm{Zn}(\mathrm{NO}}_3{)}_2$ and $T_{Dec}=423$ K for ${\mathrm{Al}(\mathrm{NO}}_3{)}_3$.

4. Results and Discussion

4.1. Droplet Temperature and Concentration Fields

Rosebrock et al. [8] investigated single droplet combustion of aluminium nitrate (${\mathrm{Al}(\mathrm{NO}}_3{)}_3$), iron nitrate (${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$) and zinc nitrate (${\mathrm{Zn}(\mathrm{NO}}_3{)}_2$) solution droplets. They found for ${\mathrm{Zn}(\mathrm{NO}}_3{)}_2$ the droplet-to-particle route as the preferred mechanism, whereas ${\mathrm{Al}(\mathrm{NO}}_3{)}_3$ and ${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$ solution droplets undergo microexplosions and the resulting particles indicate gas-to-particle conversion. To evaluate the mass transfer model for particle synthesis from low cost metal nitrate precursors, the simulation results in this chapter are compared qualitatively to the experiments from Rosebrock et al. [8]. To our knowledge, there is no measurement method to resolve concentration fields of burning droplets during flame spray synthesis. Therefore, we evaluate if calculated temperature and concentration fields provide an explanation for the different combustion behaviour of the metal nitrate precursor solutions as observed by Rosebrock et al. In the following, temperature and concentration fields are simulated for an initially 100 $\mathsf{\mu }$$\mathrm{m}$ sized droplet containing $0.5$ $\mathrm{mol}$/$\mathrm{L}$ ${\mathrm{Al}(\mathrm{NO}}_3{)}_3$, ${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$ or ${\mathrm{Zn}(\mathrm{NO}}_3{)}_2$ ethanol solutions. Droplets are burning at 2300 K.

In Figure 6, the size evolution for all three materials is shown. Immediately, liquid-phase precipitation begins in ${\mathrm{Al}(\mathrm{NO}}_3{)}_3$ and ${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$ solutions, since they are oversaturated at the beginning. Then, all materials follow initially the nearly identical surface regression. This is a result of preferential evaporation of ethanol. Around ${\mathrm{D}}^2/{\mathrm{D}}_0^2=0.6$, all three materials experience superheating in the region below the surface as shown in Figure 7a). When the droplet surface is depleted of ethanol, droplets obtain $T_{Dec}$ of their respective precursor. This is accompanied with a change of regression rate in the region below ${\mathrm{D}}^2/{\mathrm{D}}_0^2=0.4$. The results also show that liquid-phase precipitation precedes precursor decomposition. In Figure 6b), the influence of droplet relative velocity $v_{rel}$ is shown. Due to the increased regression rate, $T_{Dec}$ is reached only at the very end of droplet lifetime. If, in the case of ${\mathrm{Zn}(\mathrm{NO}}_3{)}_2$, droplet-to-particle formation is preferred as reported by [8], the resulting oxide particles will be smaller.

Figure 7a,b show the temperature fields and concentration fields at the onset of superheating. Superheating occurs in proximity to the droplet surface, enabling homogeneous vapor nucleation inside the droplet for all three precursors. On the other side, only ${\mathrm{Al}(\mathrm{NO}}_3{)}_3$ and ${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$ solution droplets reach supersaturation at onset of superheating. The supersaturation is high at the droplet surface und low in the interior. This indicates that precipitates will form preferentially on the surface, meeting the requirement for shell formation, whereas precipitates inside the droplet enable heterogeneous vapor nucleation. Nucleated vapor bubbles then expand while being heated up, leading to a rupture of the surface layer and causing droplet fragmentation. The combination of a diffusional barrier of surface precipitates and heterogeneous and homogeneous vapor bubble nucleation could explain why micro explosions are observed in ${\mathrm{Al}(\mathrm{NO}}_3{)}_3$ and ${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$, while these are absent in the case of ${\mathrm{Zn}(\mathrm{NO}}_3{)}_2$ [8]. If vapor bubble nucleation occurs in the case of ${\mathrm{Zn}(\mathrm{NO}}_3{)}_2$, the bubbles can leave the droplet without causing disruption, since no pressure build up takes place due to a lack of a shell.

The interpretations above differ from those of Rosebrock et al. [8,26] and Li et al. [9] in the form that precursor decomposition is not required to form a viscous shell. Especially in the case of nitrates with comparatively low solubility in ethanol, such as ${\mathrm{Al}(\mathrm{NO}}_3{)}_3$ and ${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$, precipitates may act as nucleation sites for vapor bubbles and form a viscous shell.

To conclude the discussion above, droplet fragmentation event can be decided in the framework of DPBMC simulation when two consecutive events take place. First, liquid-phase precipitation or precursor decomposition needs to occur. Then, internal boiling needs to happen in the form of heterogeneous or homogeneous vapor nucleation. Widiyastuti et al. [27] modeled spray pyrolysis under low-pressure conditions and concluded, from the obtained very high evaporation rates, that droplets are ruptured nearly instantaneously, releasing all monomers at once. This might prove a feasible approximation in the case of violent droplet disruptions.

4.2. Droplet Population Evolution

To investigate the influence of the DSD on particles, nanoparticle synthesis in a flame spray burner is simulated. The ambient pressure and temperature are 1 bar and 2300 K with initial droplet parameters given in Appendix A. As a model system, ${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$ in ethanol is chosen. The relative velocity is approximated by Equation (23) for droplets between 5 $\mathsf{\mu }$$\mathrm{m}$ and 100 $\mathsf{\mu }$$\mathrm{m}$ and set to 0 $\mathrm{m}$/$\mathrm{s}$ for particles below 5 $\mathsf{\mu }$$\mathrm{m}$ and to 80 $\mathrm{m}$/$\mathrm{s}$ for particles bigger than 100 $\mathsf{\mu }$$\mathrm{m}$. This is close to Phase Doppler Anemometry (PDA) measurements [3]

$$v_{Rel}=\frac{1}{s}0.8\times {10}^6·\mathrm{D}.$$

(23)

To illustrate the influence of droplet polydispersity on particle formation, only the heat conduction model following from Equations (5) with (7) is applied, while a uniform precursor concentration inside the droplet is assumed. All DSDs are assumed to be initially lognormal distributed. The initial mean droplet diameter ${\mathrm{d}}_{\mathrm{g},{\mathrm{D}}_0}$ equals 10 $\mathsf{\mu }$$\mathrm{m}$ in all simulations, while the initial geometric standard deviation ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}$ of the DSD is set to 1.1, 1.3 and 1.7. Flame spray synthesis of nanoparticles follows three distinct stages. First, only droplets are present. Second droplets and particles coexist in the lower flame region. In the third stage, only particles are left. However, earlier publications in PB modeling of FSP applied restrictions on droplet behaviour. Mueller [6] assumed instant droplet evaporation, while Widiyastuti et al. [4] used constant droplet concentrations. With these approaches, these three distinct stages can not be discriminated. By applying the DPBMC technique, no further restrictions or simplifications regarding the droplet phase have to be imposed. The resulting overall particle and droplet dynamics are exemplary shown in Figure 8 for an initial ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}=1.3$.

At HAB = 0 mm, only droplets exist. The particle nucleation mode is formed when droplets first start to evaporate and produce monomers at HAB = 0.1 mm. Evaporation of small droplets leads to a decrease of droplet concentration and a broadening of the DSD between HAB of 0.1 mm to 5 mm (Figure 8a). At the same HAB, self-coagulation of the nucleation mode leads to the formation of first larger particles (Figure 8b). Continuous nucleation and coagulation continues until the majority of droplets has burned off at approximately HAB = 40 mm and thus halts nucleation. Shortly before this point, ${\mathsf{\sigma }}_{\mathrm{g},\mathrm{P}}$ reaches its maximum value. Continuous build up of the nucleation mode with simultaneous coagulation results in a bimodal PSD and consequently a high ${\mathsf{\sigma }}_{\mathrm{g},\mathrm{P}}$. From here, no droplets exist anymore and particle coagulation dominates the temporal evolution of the PSD. In the end, the PSD reaches the self-preserving form.

The influence of transient heating on a droplet population is shown in Figure 9. Small droplets heat up first and start evaporating. Following that, bigger droplets heat up and start to evaporate at around 50 mm above the burner. This results in a broadening of the DSD, while the droplet concentration continuously decreases. In the end, only very few droplets between 30 $\mathsf{\mu }$$\mathrm{m}$ and 100 $\mathsf{\mu }$$\mathrm{m}$ exist. Therefore, if the DSD consists mainly of large droplets, as is the case for ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}=1.7$, the longer heating time will delay droplet evaporation and particle formation.

The resulting mean droplet diameter ${\mathrm{d}}_{\mathrm{g},\mathrm{D}}$ and the normalized spray load ${\mathrm{M}}_{\mathrm{D}}/{\mathrm{M}}_{\mathrm{D},0}$ over HAB are shown in Figure 10a,c for all ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}$. Broader DSDs contain more large droplets, which have a longer lifetime than small droplets. This is also shown by Figure 10b as SMD stays unchanged and then increases as the weighted average shifts to larger droplets. In contrast to that, SMD decreases for ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}$ = 1.1 initially. This is the expected behaviour if only one size class for the droplet phase is simulated, as droplets can only decrease in size due to combustion. Heine and Pratsinis [3] show PDA measurements where SMD increases initially during ${\mathrm{ZrO}}_2$ synthesis by FSP. According to them, the initial increase of SMD is due to evaporation of smaller droplets, while droplet coagulation or droplet dispersion also might affect SMD. The simulation results in Figure 10c show that this behaviour can be explained by combustion of small droplets, although it is highly dependent on ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}$.

Finally, Figure 11a–d show the PSD properties over time for the each initial DSD. Particle concentration is initially the lowest for ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}$ = 1.7 due to the lowest nucleation rate in the case of ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}$ = 1.7, whereas, for ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}$ = 1.1 and 1.3, particles nucleate in a burst-like fashion as shown in Figure 11b. The spikes in ${\mathsf{\sigma }}_{\mathrm{g},\mathrm{P}}$ are caused by the formation of a nucleation mode, which are more pronounced for ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}$ = 1.3 and 1.7. A prolonged nucleation of monomers thus result in particles with higher polydispersity. After nucleation slows down, coagulation dominates the growth process until ${\mathsf{\sigma }}_{\mathrm{g},\mathrm{P}}\approx 1.46$ is reached. This is expected when the self-preserving form is reached in the free molecular regime. At around HAB = 800–1000 mm particle concentration, ${\mathrm{d}}_{\mathrm{g},\mathrm{P}}$ and ${\mathsf{\sigma }}_{\mathrm{g},\mathrm{P}}$ reach the same value for all three simulations as a direct consequence of attaining the self-preserving form.

Depending on the particle residence time in the reactor, the PSD may be quenched by exhaust gas before coagulation depletes the nucleation mode, resulting in a bimodal PSD. Otherwise, the PSD becomes independent of the initial DSD. Similar results were found by Heine and Pratsinis [3] for monodisperse droplet calculations. Large droplets of 40 $\mathsf{\mu }$$\mathrm{m}$ release monomers longer over time because of their higher lifetime compared to smaller droplets and consequently lead to broader PSDs. They calculated for 40 $\mathsf{\mu }$$\mathrm{m}$ droplets a ${\mathsf{\sigma }}_{\mathrm{g},\mathrm{P}}$ over 2.2 in the simulated time. They also found that self-preserving form is reached relatively quickly once all droplets evaporate. It can be concluded that prolonged monomer release promotes the formation of bimodal PSDs. To describe the monomer release and nucleation rate correctly, the droplet polydispersity has to be taken into account.

For all simulations, isothermal flame temperature is assumed. To test this assumption, simulations with flame temperatures of 1700 K, 2000 K and 2300 K are calculated and evaluated with respect to ${\mathrm{J}}^{*}$ and single droplet results regarding surface temperature $T_R$ and surface regression rate $\alpha$ for an initially 100 $\mathsf{\mu }$$\mathrm{m}$ droplet. In all three cases, the droplet reaches its maximum temperature at around HAB = 300 mm (Figure 12b)), which is slightly below ethanol boiling temperature. Here, the droplet also attains its maximal $\alpha$ in each case. For the flame temperature 1700 K and 2300 K, the maximum $\alpha$ equals approximately $1.75\times {10}^{-6}$ ${\mathrm{m}}^2/\mathrm{s}$ and $2.1\times {10}^{-6}$ ${\mathrm{m}}^2/\mathrm{s}$ respectively, which is fairly similar considering 600 K difference in flame temperature. It can be concluded that ${\mathrm{J}}^{*}$ is not sensitive regarding the flame temperature.

However, the flame temperature has a major influence on primary particle size as shown by Heine and Pratsinis [3]. When particles experience a temperature drop after leaving the flame region, primary particle growth stops while agglomerate growth continues. In the case of ZnO synthesis, the flame may exceed the decomposition temperature of ZnO (2250 K). In this case, ZnO particles dissociate into Zn and oxygen and may form ZnO again. Rosebrock et al. [8] state this to be the reason for different morphology of particles produced in single droplet experiments compared to particle synthesis by FSP. As far as only ${\mathrm{J}}^{*}$ is concerned, the isothermal assumption can be applied leading to a negligible error in ${\mathrm{J}}^{*}$.

5. Conclusions

The Monte Carlo method presented in this work allows for simulating the temporal evolution of the particle size distribution coupled to the temporally evolving droplet size distribution. Droplets and nanoparticles are rendered by weighted MC-particles and MC-droplets, respectively, which allow for taking the full size spectrum of the droplet size distribution into account. In contrast to nanoparticles, droplets experience a velocity differential relative to the gas phase. This relative velocity leads to an increased droplet evaporation rate. Using MC-droplets, the individual relative velocity can be applied to each fraction of the droplet size distribution. In this way, measured droplet velocity distributions can be applied to the droplet size distribution. Additionally, MC-droplets can be used to track several physical properties of the respective droplet at once. The internal temperature and concentration fields can be used to predict whether a droplet might experience a fragmentation event. The temperature-and concentration field calculations show that internal boiling occurs below the droplet surface for ${\mathrm{Al}(\mathrm{NO}}_3{)}_3$, ${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$ and ${\mathrm{Zn}(\mathrm{NO}}_3{)}_2$ solution droplets. However, at the onset of internal boiling, only ${\mathrm{Al}(\mathrm{NO}}_3{)}_3$ and ${\mathrm{Fe}(\mathrm{NO}}_3{)}_3$ solution droplets show a possible formation of a viscous shell by surface and liquid phase precipitation of the precursor. This might explain why, for ${\mathrm{Zn}(\mathrm{NO}}_3{)}_2$ solution droplets, no disruptive burning is observed in experiments by Rosebrock et al. [8]. As a result, a droplet fragmentation event is defined as surface precipitation or decomposition of precursor followed by internal boiling below the droplet surface. For implementation of breakage algorithms in the framework of Monte Carlo modeling, more research is needed to find suitable breakage rates and functions. This will likely further improve the calculations of mass transfer from droplets into the gas phase.

Calculating temperature and concentration fields can help to decide whether the gas-to-particle or droplet-to-particle route is preferred. Additionally, by tracking the surface temperature of each droplet, their evaporation rate and thus precursor release rate can be calculated more precisely, which is shown to vary greatly along the reactor axis depending on initial droplet polydispersity. The usage of weighted MC-particles for the droplet phase allows for easily treating removal of evaporated droplets, by setting their weight to zero. This way, no simplifying assumptions to the droplet phase need to be applied (e.g., instant droplet evaporation or constant droplet evaporation) as used in earlier publications regarding flame spray synthesis of nanoparticles. As a result, three distinct stages of flame spray synthesis of nanoparticles can be distinguished. The first stage consists of a droplet only region at approximately HAB = 1 mm. After the first particles nucleate, a mixed zone can be observed where droplets and particles coexist. In the third stage at HAB = 40 mm, nearly all droplets have burned off and particles remain.

The initial droplet polydispersity also shows that it has a large impact on the particle nucleation rate and polydispersity. Higher initial droplet polydispersities promote the formation of bimodal PSDs. In the case of ${\mathsf{\sigma }}_{\mathrm{g},{\mathrm{D}}_0}$ = 1.7, an initial increase of the Sauter mean diameter of the droplet phase is observed. This is a result of small droplets evaporating and shifting the weighted average to larger droplet diameters. This behaviour is frequently observed in PDA measurements of the droplet phase in flame spray synthesis and can only be simulated if the droplet size distribution is assumed to be polydisperse. Overall, it is shown that droplet phase calculations benefit largely from parallel computations on the GPU. The computation time for 100 MC-droplets per single Monte Carlo simulation compared to 2000 MC-droplets per simulation roughly doubles, while nearly eliminating evaporation noise for the droplet phase. If the number of MC-droplets is optimized regarding the particle nucleation rate and computation time, 100 MC-droplets per single Monte Carlo simulation are shown to be sufficient.