Computational Modeling of Silicon Nanoparticle Synthesis: I. A General Two-Dimensional Model

Abstract

A two-dimensional model has been developed for silicon nanoparticle synthesis by silane thermal decomposition driven by laser heating in a tubular reactor. This fully coupled model includes fluid dynamics, laser heating, gas phase and surface phase chemical reactions, and aerosol dynamics which includes particle transport and evolution by convection, diffusion, thermophoresis, nucleation, surface growth, and coagulation processes. A moment method, based upon a lognormal particle size distribution, and a sectional method are used to model the aerosol dynamics. The simulation results obtained by the two methods are compared. The sectional method is capable of capturing the bimodal behavior that occurs locally during the process, while the moment method is computationally more efficient. The effect of operating parameters, such as precursor concentration, gas phase composition, inlet gas velocity and laser power input, on the characteristics of the particles produced are investigated. Higher temperature generates more large particles with higher precursor conversion. Shorter residence time, from high inlet velocity, produces more small particles at the cost of lower precursor conversion. Increasing H₂ concentration suppresses particle formation by reducing the rates of gas phase and surface reactions, leading to fewer and smaller particles. In addition, the relative importance of the interconnected mechanisms involved in the particle formation is considered. The results make clear that spatial variations in reaction conditions are the primary source of size polydispersity and generation of non lognormal overall size distributions in a laser-driven process like that considered here.

INTRODUCTION

Synthesis of silicon nanoparticles is of great interest because of their unique optical and electronic properties. Various methods, such as silane thermal decomposition (Li et al. 2003; Ostraat et al. 2001) and laser ablation (Seto et al. 2001), have been used to produce silicon nanoparticles of high quality. Fundamental understanding of the various interconnected mechanisms involved in the particle formation, such as gas phase and gas-surface phase chemical kinetics and particle size evolution through nucleation, growth, coagulation and coalescence, would be of great value in designing and optimizing processes for producing silicon nanoparticles. On the other hand, the same understanding can also contribute to contamination control in the semiconductor industry, where particles formed by homogeneous nucleation within the processing environment are rapidly becoming the most important source of yield loss as integrated circuits become smaller. Moreover, silicon particle formation by thermal decomposition of silane is a prototypical example of a system where particle nucleation and growth are believed to occur via a complex network of chemical reactions, as opposed to nucleation from, and condensation of, a supersaturated vapor.

The evolution of nanoscale particles dispersed in a gas is governed by the aerosol general dynamics equation (GDE). The GDE (Friedlander 2000) is an integro-partial-differential equation whose analytical solution is unavailable except for extremely idealized conditions. The lack of analytical solution necessitates the use of numerical techniques involving different degrees of approximation at different computational costs. Among them, the moment method and the sectional method are frequently used in the field of aerosol dynamics.

The method of moments (Friedlander 2000) describes the aerosol dynamics through the evolution of a set of moments of the particle size distribution. The k^th moment of a particle size distribution is defined as:

where n(v, t) is the particle size distribution function. The evolution of the k^th moment of particle size distribution is obtained by multiplying both sides of GDE by v^k and integrating over all particle sizes (Randolph and Larson 1971). The method of moments is computationally efficient while providing the most important information about the process, such as particle concentration, average size and polydispersity (Barrett and Webb 1998; Frenklach 2002; Frenklach and Harris 1987; Lee and Liu 1980; Otto et al. 1999; Park et al. 1999a; Park et al. 1999b; Seigneur et al. 1986; Settumba and Garrick 2003; Talukdar and Swihart 2004; Zhang et al. 1999). However, the method of moments usually requires a priori assumption of the shape of the particle size distribution. The lognormal function is often elected not only because many experimental results indicate that aerosol size distributions are often approximately lognormal, but also because the physical origin of lognormal particle size distributions in open, stationary systems has been provided (Kiss et al. 1999). As a result, measured aerosol size distributions are often characterized using the parameters of a lognormal function. McGraw et al. compared the aerosol properties and evolution for explicit dissimilar test distributions having identical moments (McGraw et al. 1998). It was found that despite their dissimilarity, those distributions exhibit virtually identical physical and optical properties and dynamics. Another disadvantage of the moment method is that it requires closure of the moment equations, namely all the terms in the moment equations have to be expressed in terms of the moments themselves. Such closure has been achieved for a broad range of particle sizes over continuum, gas-slip, and free molecular regimes for coagulation analytically, taking into account the effects of both particles size and the size spread (Lee and Chen 1984; Lee et al. 1984). A simple model for simultaneous aerosol nucleation, condensation, and coagulation was proposed by Pratsinis (1988). This model covered the entire particle size spectrum approximating the size distribution by a unimodal lognormal function and has high potential to be used for the description of aerosol dynamics in highly complex aerosol systems for which detailed knowledge of the particle size distribution may not be required. There are also numerous works using the method of moments without pre-assumption of the form of the particle size distribution (Frenklach 2002; McGraw 1997; McGraw and Wright 2003). However, those methods either increase the computational cost or have the drawback of being implicit with respect to the moments and particle size distribution.

The sectional method approximates the continuous size distribution by a finite number of sections or bins, within which one numerically conserved aerosol property is held constant. The use of the sectional method in some applications is desirable because it offers generality and flexibility in describing the evolution of an aerosol, but it is computationally more expensive than the moment method. The computational expense of the sectional method mainly comes from the complicated integrals in the coagulation term in the GDE when size-dependent inter-bin and intra-bin coagulation coefficients are used (Gelbard et al. 1980; Landgrebe and Pratsinis 1990). A simpler approach of using a constant collision kernel within each bin was used to simplify the calculation and decrease the computational cost. (Hounslow et al. 1988; Kumar and Ramkrishna 1996a; Lister et al. 1995; Talukdar and Swihart 2004). Potential difficulties with the sectional method include numerical diffusion and stability problems when fixed sections are used and particle growth is important. A moving sectional method can overcome these problems by varying the section size for each interval to adapt the size distribution during the aerosol evolution process (Kumar and Ramkrishna 1996b). An effective moving sectional model has been used to simulate the flame synthesis of titania nanoparticles, accounting for gas phase chemical reactions, coagulation, surface growth, and sintering (Tsantilis et al. 2002). However, the moving sectional method is difficult to implement in a multidimensional flow environment like that which is the focus of the present study.

Here, we focused on developing a two dimensional model of a cylindrical reactor for synthesizing silicon nanoparticles using silane as precursor. The reactor wall is transparent to a CO₂ laser as shown in Figure 1. Silane and carrier gases (He and H₂) enter the reactor along the z direction. Silane absorbs the laser energy while the carrier gases do not. As silane is heated by the laser, it decomposes into more reactive species that further react to form larger silicon hydride clusters. When the cluster size exceeds a certain critical value, particle nuclei are formed. The silicon particles in the solid phase can then grow by surface reactions with gas phase silicon hydride clusters and by coagulation with each other. The introduction of the laser beam perpendicular to the flow breaks the axial symmetry of the system to some degree. However, the reacting flow (e.g., the coupled fluid dynamics and chemical reactions) and particle dynamics are still close to axisymmetric. In the following sections, a two dimensional model is developed with fully coupled fluid dynamics, laser heating, silane decomposition chemistry and aerosol dynamics described by either a moment method or sectional method. This is, to the best of our knowledge, the first two-dimensional model of a laser-driven aerosol reactor that includes all relevant modes of aerosol evolution: nucleation, growth, thermophoresis, diffusion, convection, and coagulation. Modeling of such a laser-driven system is especially challenging because it requires consideration of a two-dimensional domain with large temperature and concentration gradients. Moreover, detailed silane decomposition chemistry was used to link gas phase kinetics to aerosol dynamics. This inclusion of detailed nucleation and growth chemistry is needed to fully describe particle formation in this system, in contrast to many flame-based processes for producing metal oxides where a single reaction, or small set of reactions, generates a metal-oxide monomer that then irreversibly coagulates with particles and other monomers.

FIG. 1 Schematic of the tubular reactor (not to scale).

THEORY

In the two-dimensional model, we assume that the boundary layer approximations are valid (Coltrin et al. 1989; Coltrin et al. 1984; 1986). That is, the governing equations can be reduced to a system of parabolic partial differential equations describing the conservation of mass, momentum, energy, chemical species, and aerosol particles. The approximation relies on the existence of a principal flow direction in which convective transport is dominant and diffusive effects are negligible, which corresponds to the z direction here. In this approach the pressure is impressed on the flow by the boundary conditions and is uniform in the cross stream direction. The imposition of pressure replaces the cross-flow momentum equation entirely. The equations that are used to describe the fluid dynamics are as follows:

Mass conservation:

Axial momentum conservation:

Pressure:

State:

where, z and r are the axial and the radial directions of the cylindrical coordinate system. u and v are the axial and radial velocities, respectively. Because of the low operating pressure used, 0.55 atm, the ideal gas law is used here with denoting the average molecular weight of the reacting flow, R being the gas constant, and T being the operating temperature. The density ρ refers to the gas phase. Due to the large amount of carrier gas used in synthesizing silicon nanoparticles in the tubular reactor, the mass fraction of the particles is insignificant in the reacting flow and the density of the aerosol is assumed to be the same as the density of the gas.

The driving force for silane pyrolysis in this system is the continuous CO₂ laser beam of 0.63 cm diameter, centered at 0.4 cm above the inlet. As mentioned above, silane absorbs the laser energy as the reacting flow crosses the laser in the reaction zone. The energy conservation equation is thus:

where C_p is the molar specific heat of the gas mixture, h_k is the specific enthalpy of gas phase species k, and λ is the thermal conductivity of the gas mixture. On the left-hand side (LHS) are the convective terms. The first term on the right-hand side (RHS) is for conduction; the second term arises from species radial diffusion; the third term represents chemical reactions taking place in the gas phase (WDOT _k ) and on the surface of particles (SDOT _k *SMOM) between gas phase and surface phase species. SMOM is the total surface area of aerosol per unit volume.

Note that the laser energy has a Gaussian distribution which is axisymmetric along the laser direction x as described by Equation (7):

where I_o is a coefficient that incorporates both laser intensity and the absorption coefficient of the silane. However, the two dimensional model of the tubular reactor is axisymmetric along the main flow direction z. In order to maintain the symmetry of the two dimensional model, the laser intensity at any location with the same radial coordinate at specific axial position is taken be the average of the energy intensities at those locations at specific z, as described by the Equation (8) for the last term in Equation (6):

where C _SiH4 is the molar concentration of silane, d is the laser beam diameter (0.63 cm), z_o is the z coordinate of the laser center (0.4 cm), and Bessel_Io is the modified Bessel function of the first kind.

Absorption of the laser energy by the silicon particles is neglected in this work, because silicon absorbs very weakly at the laser wavelength of 10.6 μm, and because the particle size is extremely small relative to the laser wavelength. The extinction coefficient (imaginary part of the refractive index) of silicon at this wavelength is about 2 × 10^–4. The particle size is always at least 2 orders of magnitude smaller than the laser wavelength, so they are far into the Rayleigh scattering regime. For larger particles, or less transparent particles, laser heating of the particles could be important and would need to be addressed. Here, it is expected to be negligible compared to absorption by the SiH₄ precursor.

Silane decomposition is the starting point for making silicon nanoparticles in the tubular reaction. Numerous studies of silane decomposition kinetics have been reported (Becerra and Walsh 1992; Ho et al. 1994; Moffat et al. 1992a, b; Purnell and Walsh 1984; Wong et al. 2004). This work used an extensive chemical kinetic mechanism for silicon hydride cluster formation during silane pyrolysis (Girshick et al. 2000; Nijhawan et al. 2003; Swihart and Girshick 1999). This reaction mechanism includes 133 gas phase species with a critical size of 10 silicon atoms. It consists of reversible reactions among silicon hydrides containing up to 10 silicon atoms and irreversible formation of silicon hydrides containing 11 or more silicon atoms. It provided a means of calculating a particle nucleation rate that can be used as the nucleation source term in the aerosol dynamics models that predict particle formation, growth, and transport. The governing equation for each gas phase species is as follows:

where Y_k is the mass fraction of gas phase species k. j _k,r is the diffusive flux of species k in the radial direction. W_k is the molecular weight of species k, and WDOT _k and SDOT _k · SMOM are the production rates of species k by gas phase reactions and surface phase reactions (on the surface of particles).

There are two surface species, SiH(s) and Si(s), and one bulk phase species Si(b), with three types of surface reactions, with rate parameters from (Ho et al. 1994)

Saturated silanes (in which each Si atom is bound to four other atoms) react only with “bare” surface sites, and with a temperature dependent rate constant, while unsaturated silenes and silylenes react with any surface site, with unity sticking coefficient. In the current system, SiH ₄ is the most important species for surface reaction because of its dominant concentration among other gas phase silicon hydride clusters. To simplify the calculation, based on the assumption that the surface reactions are fast compared with gas phase reactions, the surface species Si(s) and SiH(s) are assumed to be in instantaneous equilibrium with respect to the dominant gas phase species SiH ₄.

The method of moments with a unimodal lognormal distribution is first used in this study to describe the aerosol dynamics. We assumed instantaneous coalescence between two coagulating particles, i.e., particles are treated as spherical. The conservation equation for the first three moments of the particle size distribution is as follows:

where, v_p is the particle volume and k is the index for the k^th moment. The first two terms on the LHS of Equation (11) represent the convective transport of the particle moments. The remaining two terms on the LHS are the particle diffusion and thermophoresis in the radial direction with j_vp being the diffusive flux and v_th being the thermophoretic velocity. An accurate expression for the particle diffusion coefficient is given by Friedlander (2000), but it does not facilitate closure. Instead, an approximate expression spanning the transition regime is used (Lee and Liu 1980) to achieve the closure of the particle diffusion term:

One of the features of the tubular laser reactor is rapid heating and cooling which results in huge temperature gradients in the reaction zone. Therefore, thermophoresis must be considered in this work. The thermophoretic velocity v_th is independent of particle sizes for small particles (d <λ) and is given by Waldmann and Schmitt (1966) as

For larger particles, v_th is dependent on the thermal conductivity of the particle material. However, this effect is negligible in the particle size range considered here. Thus,

The source term on the RHS of Equation (11) is related to particle nucleation as a result of silane decomposition, particle surface growth through chemical reactions, and coagulation between particles:

In the nucleation term, Equation (17), I is the sum of the production rates of silicon hydride clusters bigger than the critical size according to Swihart (Swihart and Girshick 1999). v* is the volume of the particle of critical nucleus size:

In the free molecular regime, the particle growth rate is expressed as Equation (18), with G_linear being the linear growth rate defined as the production rate of bulk silicon per unit surface area by surface chemical reactions:

Similarly, in the continuum regime:

which reflects the potential for diffusion limitations on the growth rate. The growth term in the particle moment equation is taken to be the harmonic mean of values in the free molecular regime and the continuous regime:

The coagulation term in the moment equation is also calculated from the harmonic mean of the values in the free molecular regime and continuum regime following Pratisnis' (1988) formulations in each of the two regimes.

The initial conditions at the inlet of the tubular reactor are taken to be constant temperature and fully developed laminar flow with a parabolic velocity profile. Non-slip boundary conditions are used at the reactor wall.

The above governing equations for mass, momentum, energy, gas phase species and particle moment evolution are a system of coupled partial differential equations (PDEs). The method of lines is used to discretize the radial derivatives to reduce the PDEs to a system of ordinary differential equations (ODE), where the central difference approach is used for the first and the second order radial derivatives, as described by Coltrin et al. (1984). DASSL (Brenan et al. 1989), which is a solver for large scale differential algebraic equations is used to solve the above ODEs simultaneously. The Chemkin family of reacting flow codes (Coltrin et al. 1996; Kee et al. 1986; Kee et al. 1996) was used for computing reaction rates and thermodynamic and transport properties.

RESULTS

Moment Method

The optimal number of nodes in the radial direction that is necessary to achieve accurate solutions at reasonable cost was first determined. Figure 2 shows the convergence of the simulation with increased number of nodes. The labels on the horizontal axis, for example , represent the maximum or average relative difference in results between simulations with the indicated number of nodes in the domain. Only the zeroth moment, temperature and velocity are plotted here. It can be seen from Figure 2a that the maximum difference in the zeroth moment (total particle number concentration) is less than 5% between simulations with 121 nodes and 101 nodes. Changes in other solution variables, such as temperature and velocity, are much smaller. Of course, at most locations, the differences are much smaller than the maximum values, as indicated by Figure 2b. Further increasing the number of nodes used does not increase the accuracy significantly, and therefore the increase in computational cost for additional nodes was not justified. We thus concluded that the node independency in the radial direction was sufficiently achieved with 101 nodes and this is the value used throughout the remainder of this work. It should be noted that while, in some contexts 5% variation in a solution variable as observed for M₀ would be quite significant, in this context it represents quite tight convergence. M₀ changes by many orders of magnitude across very small distances in the reactor (vide infra) and the remaining differences are confined to those regions of very steep gradients.

FIG. 2 Node independency study for discretization in the radial direction, showing the maximum (a) and average (b) relative changes in solution variables with changing number of nodes. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10 mol% SiH₄, 80 mol% H₂ and 10 mol% He; I_o =1.22 × 10¹⁵ erg/(mol s).

Figure 3 shows the temperature profile for a cross section of the tubular reactor. With the energy provided by the laser, the reacting flow is heated from 300 K at the inlet (z= 0) to the peak temperature of 1086 K for this case within a distance of 0.5 ∼ 0.6 cm, which corresponds to a residence time of less than 10 milliseconds. The huge temperature gradients around the center of the region are caused by the localized energy distribution of the laser. Figure 4 shows the silane mole fraction profile over the domain. At the center of the region, silane decomposes into more reactive species and further forms larger silicon hydride clusters and solid particles, where silane concentration decreases. When the temperature decreases downstream after passing the laser, silane concentration recovers, through the inward diffusion of silane from the edge.

FIG. 3 Temperature profile of the two-dimensional tubular region with inlet at z= 0. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10 mol% SiH₄, 80 mol% H₂ and 10 mol% He; I_o = 1.22 × 10¹⁵ ergs/(mol s).

FIG. 4 Precursor concentration profile of the two-dimensional tubular region with inlet at z= 0. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10 mol% SiH₄, 80 mol% H₂, and 10 mol% He; I_o =1.22 × 10¹⁵ ergs/(mol s). (Figure provided in color online.)

Figure 5a shows the particle number concentration profile, namely the zeroth moment. Particle number concentration peaks around z≈ 0.7 cm, compared to about 0.6 cm for the temperature peak. At high temperature, particle nucleation is strong, but as the particle number concentration increases, coagulation becomes important. As the decrease in particle number concentration by coagulation overcomes the increase in particle number concentration by nucleation, the particle number concentration reaches its peak and starts to decrease. In contrast, the particle total volume concentration increases throughout the process. This is because the processes of particle formation and growth by nucleation and surface reactions are irreversible. Note that there is a valley in particle number concentration along the centerline downstream of the peak at z≈ 0.7 cm in Figure 5a. Even though temperature is highest along the centerline, the rapid depletion of the precursor around the location where the laser passes through lowers the nucleation rate there. Figure 5b is the particle number concentration calculated under the same operating conditions, except that the aerosol dynamics is described by the sectional method which will be discussed in the next section. Figures 7b and 8b are plotted for the same reason.

FIG. 5 Particle number concentration profile of the two-dimensional tubular region with inlet at z= 0. (a) is calculated with moment method, and (b) is calculated with sectional method. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10 mol% SiH₄, 80 mol% H₂, and 10 mol% He; I_o =1.22 × 10¹⁵ ergs/(mol s). (Figure provided in color online.)

FIG. 7 Particle geometric mean diameter profile within the two-dimensional tubular region with inlet at z= 0. (a) is calculated with moment method and (b) is calculated with sectional method. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10 mol% SiH₄, 80 mol% H₂ and 10 mol% He; I_o =1.22 × 10¹⁵ ergs/(mol s). (Figures provided in color online.)

FIG. 8 Particle geometric standard deviation profile of the two-dimensional tubular region with inlet at z= 0. (a) is calculated with the moment method and (b) is calculated with the sectional method. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10 mol% SiH₄, 80 mol% H₂, and 10 mol% He; I_o =1.22 × 10¹⁵ ergs/(mol s). (Figures provided in color online.)

Figure 6 shows the contours of nucleation rate at the center of the tubular reactor. It can be seen that near the reactor inlet, the centerline has the highest nucleation rate compared with other radial positions at the same z coordinate. Further downstream, the nucleation rate becomes lower at the centerline compared with the off-axis locations with the same z coordinate. This results in lower particle number concentration at the centerline after the peak at about z≈ 0.7 cm as shown in Figure 5a.

FIG. 6 Contours of particle nucleation rate at the center of the tubular region. [number/(cm³ s)]. a = 3.3 × 10¹⁴, b = 6.6 × 10¹⁴, c = 9.9 × 10¹⁴, d = 1.3 × 10¹⁵, e = 1.65 × 10¹⁵, f = 1.98 × 10¹⁵, g = 2.31 × 10¹⁵, h = 2.64 × 10¹⁵.

The particles produced at conditions listed in Figures 3–5a have a mixing cup average geometric mean diameter of 7.2 nm with a geometric standard deviation of 2.2 at the outlet, which is 5.5 cm downstream from the inlet. These overall size statistics are based on taking a mixing-cup average of the size distributions, and then computing the geometric mean diameter and standard deviation of the resulting (arithmetically) averaged size distribution. The spatial distributions of the geometric mean diameter and the geometric standard deviation in the tubular region are shown in Figures 7a and 8a, respectively.

Figure 9 shows the particle size distributions along the centerline of the tubular reactor at different axial positions. The particle size corresponding to the peak of each distribution is the geometric mean diameter of the particles at that location and the area under each curve is the particle number concentration. Similar to the huge temperature gradients close to the center of the tubular reactor, there is also a sudden burst in particle concentration near the laser center. For example, as shown in Figure 9, the difference in particle number concentration between z= 0.5 cm and z= 0.6 cm is about 8 ∼ 9 orders of magnitude within a residence time of about 1 millisecond.

FIG. 9 Particle size distributions along the centerline of the two-dimensional tubular region with inlet at z= 0. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10 mol% SiH₄, 80 mol% H₂, and 10 mol% He; I_o =1.22 × 10¹⁵ ergs/(mol s).

The generation and evolution of the silicon nanoparticles are affected by multiple interconnected physico-chemical processes including convection, diffusion and thermophoresis, particle nucleation by gas phase chemical reactions, particle growth by surface chemical reactions between gas phase species and particle surfaces, and by coagulation of solid particles. Figure 10 shows the relative importance of each process to the change in particle number concentration. Because there are essentially no silicon particles produced before z≈ 0.4 cm, the graph is plotted starting from z= 0.4 cm, where particle nucleation starts to become significant. The vertical axis of the graphs indicates the contributions from different processes mentioned above to the change in particle number concentration. Because surface growth does not change the particle number concentration, the growth process is not plotted here. Figure 10a is for the centerline of the tubular reactor and Figures 10b, 10c, and 10d are for the locations at 0.1 cm, 0.2 cm, and 0.3 cm away from the centerline. Close to the centerline (a, b, and c) and near the inlet, particle nucleation is the most important as temperature ramps up rapidly. When the particle number concentration is high enough, coagulation starts to dominate as nucleation is quenched by lower temperatures after the reacting flow passes the laser. At locations far away from the centerline, such as Figure 10d, the particles are mainly transported from the inner core of the region by convection and diffusion.

FIG. 10 Contributions of different processes to particle number concentration. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10 mol% SiH₄, 80 mol% H₂, and 10 mol% He; I_o =1.22 × 10¹⁵ ergs/(mol s).

The synthesis of silicon nanoparticles in the tubular reactor is driven by the CO₂ laser. The maximum temperature of the reacting flow at the center of the tubular region increases almost linearly with increasing laser power. Higher temperature favors particle nucleation and surface growth. As a result, the number and sizes of the particles produced at the outlet (mixing cup values) showed increasing trends with increased laser power input in Figure 11, where I_o is the constant described in Equation (7). The corresponding peak temperature increases almost linearly from about 1005 K to 1115 K as I_o increases from 1.0 × 10¹⁵ to 1.3 × 10¹⁵ erg mol^–1 s^–1.

FIG. 11 Effect of laser power on the number and average size of the final particle product. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10 mol% SiH₄, 80 mol% H₂, and 10 mol% He.

The velocity of the reacting flow can affect the production rate of the particles at the outlet. At higher inlet velocity, in order to achieve the same maximum temperature at the center of the region, higher laser power must be used, because there will be more gas to be heated up per unit time. With the same maximum temperature, reacting flows with higher inlet velocities produce more particles at the outlet because of higher precursor flowrates as shown in Figure 12. However, the particles are smaller because the residence time is shorter at higher velocities. Thus, processes that make particles larger (surface growth and coagulation) have less time to contribute. Hydrogen is the byproduct of silane decomposition. Increasing H₂ concentration suppresses silane decomposition, which results in both fewer and smaller particles produced as shown in Figure 13, in which the H₂ to He ratio of the carrier gas is varied as the silane concentration is held constant.

FIG. 12 Effect of inlet velocity on the number and average size of the final particle product. Initial conditions: T = 300 K; 10 mol% SiH₄, 80 mol% H₂ and 10 mol% He; T_max=1080 K.

FIG. 13 Effect of initial H₂ concentration on the number and average size of the final particle product. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10 mol% SiH₄ with He and H₂ as carrier gases; I_o = 1.1 × 10¹⁵ ergs/(mol s). The He mole fraction is 0.9 minus the H₂ mole fraction shown.

Validation of the Lognormal Assumption with a Sectional Method

The moment method used here assumes a unimodal lognormal particle size distribution. In order to verify this assumption, the moment method in the two dimensional model was replaced by a sectional model to describe aerosol dynamics with the mass, momentum, energy, and species conservation equations unchanged as in the previous section.

The sectional method approximates the continuous size distribution by a finite number of sections or bins, within which one numerically conserved aerosol property is held constant. The sectional method can capture the details of particle evolution where the assumption of a lognormal size distribution in the moment method breaks down. The approach used here is similar to the one reported by Talukdar and Swihart (Talukdar and Swihart 2004). The difference is, in this work, for numerical convenience, the mass of particles in each bin (per mass of gas, and therefore formally a mass fraction), Y _np,i , with respect to the aerosol were considered as the solution variables and the governing equations are for mass conservation of particles in each bin. The particle size distribution can be obtained as:

where v _mid,i is the average volume of particles in bin i, ρ is the density of the aerosol which is approximated by the density of the gas phase, ρ _bulk is the density of bulk silicon, and φ _i is a basis function that is equal to one within bin i and zero elsewhere. The governing equation of aerosol dynamics for bin i is:

The terms on the LHS of the equation are for particle transport through convection, diffusion, and thermophoresis, with D _np,i and v _th,i being the diffusivity and thermophoretic velocity of the particles in bin i, respectively. On the RHS are the creation and depletion of particles in bin i by nucleation, surface growth and coagulation. GrowthLaw _i is the growth law covering both the free molecular and continuum regimes similar to the one used in the moment method. The coagulation coefficient is calculated based on Fuchs' interpolation formula (Seinfeld 1986). In Equation (22), the growth term is written in the form of a two-point approximation in particle size. However, in order to reduce numerical diffusion (in particle size) and to increase stability of the sectional method, a power essentially non-oscillatory (ENO) method is used to calculate the Y _np,i at the boundaries of bin i (Serna and Marquina 2004).

The accuracy of the sectional method depends on the number of bins used. The more bins used, the more accurate the solutions are, at the price of increased computational cost. Figure 14 presents the particle size distributions along the centerline of the tubular reactor obtained with the sectional method using different numbers of bins. At z= 0.6 cm, a significant number of particles have been produced by nucleation. The nuclei grow by surface reaction and coagulation, which results in the expansion of the distribution towards the right side. At z= 0.7 cm, nucleation is still strong, and the plateau in the distribution for particles less than 10 nm is the result of two overlapping modes of a particle size distribution which are relatively close to one another. The mode with small particle sizes is caused by nucleation, and the other mode for larger particles is produced by coagulation and surface growth. With further growth of the particles in the second mode, the two modes separate completely as shown in Figure 14c. There, it can also be seen that nuclei are diminishing because of the quenched nucleation at lower temperature and the continuous growth of particles. At the outlet at z= 5.5 cm, the second mode becomes dominant, with a tail of small particles at significantly lower concentration. Calculations with different numbers of bins showed similar trends. The results between simulations with 50 and 60 bins were almost identical. As a result, 50 bins were used in the sectional method for all remaining simulations and comparisons with the moment method. With 50 bins, the computational cost of simulations using the sectional method was about twice that of simulations using the moment method.

FIG. 14 Particle size distributions along the centerline of the tubular region. 20, 30, 40, 50, and 60 bins are used in the sectional method. The distributions calculated from the moment method are plotted for comparison. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10% (mass) SiH₄, 80% (mass) H₂ and 10% He (mass); I_o =1.22 × 10¹⁵ ergs/(mol s). As the number of bins increased from 20 to 60, the spacing factor for the bins decreased from 2.46 to 1.35.

The particle size distributions along the centerline predicted by the moment method are also plotted in Figure 14. Near the beginning of the process (z= 0.6 cm), the moment method and the sectional method predict similar distributions, except for very small particles. At this time, nucleation is dominant, and surface growth and coagulation are starting to show impact on the distribution. It is expected that, at the beginning of the process, particle nuclei have the highest concentration, as predicted by the sectional method. The moment method fails to capture this feature because of the assumption of a lognormal distribution. Further down the reactor (z= 0.7 cm) the difference between the two methods become more obvious: the sectional method is able to capture the high concentration of nuclei at that position, while the moment method under-predicts the concentration of particles with small sizes. Further downstream (z= 1.0 cm), the difference between the two methods becomes less significant again. At the outlet of the tubular region, (z= 5.5 cm), the sectional method predicts a similar particle size distribution to the moment method, except for the small tail of very small particles, with very low number concentration.

The distribution profiles of particle number concentration, geometric mean diameter, and geometric standard deviation in the tubular reactor calculated with the sectional method are plotted in Figures 5b, 7b, and 8b, respectively. The profiles of particle number concentration and geometric average diameter predicted by the moment method and the sectional method are nearly identical. However, there are significant differences in the distributions of the geometric standard deviation between the two methods. The moment method tends to underestimate the value of the geometric standard deviation because of the lognormal distribution assumption in particle sizes.

Figure 15a shows the particle number and volume flow rates in the reactor. The horizontal axis is the flow direction. The results from the sectional method and the moment method are very close to each other, with a maximum difference of 5% in number flowrate and 6% in volume flowrate throughout the process.

FIG. 15 (a) Particle number and volume flowrates predicted by the moment method and sectional method. (b) Arithmetic and geometric average diameters calculated by moment method and sectional method. (c) Geometric standard deviation calculated by moment method and sectional method. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10% (mass) SiH₄, 80% (mass) H₂, and 10% He (mass); I_o = 1.22 × 10¹⁵ ergs/(mol s).

Figure 15b shows the particle diameters averaged over the cross section perpendicular to the flow direction. d_MOM and d_SM are the arithmetic average diameters of the particles calculated from the moment method and the sectional method respectively. Similarly, dg_MOM and dg_SM are the geometric average diameters of the particles from the two methods. The arithmetic average diameters predicted by the two different models are almost identical, which is not surprising since the particle number and volume flowrates are in good agreement between the two models. On the other hand, the geometric average diameters predicted by both methods are smaller than their corresponding arithmetic average values. This can be explained by the definition of geometric average diameter, which gives more weight to smaller particles. Note that in Figure 15b, the geometric average diameters from the sectional method decrease slightly from z= 2.1 cm. Figure 15a has shown that the flowrate of the number of particles through this region (z> 2.1 cm) decreases while the particle volume flowrate increases. A physically meaningful average diameter should increase in this situation. However, the definition of the geometric average diameter allows it to decrease even though the particle volume flowrate increases and the particle number flowrate decreases. This illustrates the challenges inherent in describing a size distribution with parameters characteristic of a lognormal distribution, when the true distribution is far from lognormal. Therefore, two different kinds of average diameters, arithmetic and geometric, may have to be considered simultaneously to provide a physically meaningful view of the average size of particles. Figure 15c shows the geometric standard deviation of the mixing-cup-averaged size distribution, which is much larger than the local values of the geometric standard deviations plotted in Figure 8, because spatial non-uniformity contributes more to the polydispersity of the final size distribution than local effects like nucleation and coagulation.

For the moment method, the particle size distribution at any point of the tubular region is assumed to be lognormal, but the distribution of all the particles over the cross section perpendicular to the flow direction is not constrained to be lognormal. Figure 16 shows the size distributions of particles collected at the outlet calculated by both the moment method and sectional method. It can be seen the distributions differ the most for small particles. This is the cause of the difference in geometric average diameter predicted by those two methods.

FIG. 16 Particle size distributions of the particles at the outlet calculated by the moment method and sectional method. Initial conditions: T = 300 K; parabolic velocity profile with a maximum value of 80 cm/s; 10% (mass) SiH₄, 80% (mass) H₂, and 10% He (mass); I_o =1.22 × 10¹⁵ ergs/(mol s).

Even though the moment method could not capture the bimodal particle size distribution and tends to slightly over predict particle number and volume concentrations, the moment method provides results at the outlet of the tubular reactor that are very similar to those from the sectional method. Because of the relatively low computational cost, the moment method might be more suitable for simulations with more complicated reactor geometry as well as for situations where additional effects, such as finite sintering time and the resulting distribution in both size and surface area, are to be included.

CONCLUSION

We have presented a two dimensional model for generation of silicon nanoparticles in a tubular reactor with fully coupled fluid dynamics, laser heating, detailed silane decomposition chemistry and aerosol dynamics described by a moment method as well as by a sectional method for comparison in a limited number of cases. We found that even though the moment method cannot capture the bimodal behavior that occurs locally during the process, its predictions for the character of the final product at the outlet of the reactor are very close to those of the sectional method, with much lower computational cost. Higher temperature increases the production rate of the solid particles and the precursor conversion at the outlet, and results in larger particles. Shorter residence time from high inlet velocity produces more small particles at the cost of lower precursor conversion. Increasing H₂ concentration suppresses particle formation by reducing the rates of gas phase and surface reactions, leading to fewer and smaller particles.

Acknowledgments

We thank the U.S. National Science Foundation (CTS-0500249) for financial support.