Computational Fluid Dynamics of a Cylindrical Nucleation Flow

Article pubs.acs.org/JPCA

Computational Fluid Dynamics of a Cylindrical Nucleation Flow Reactor with Detailed Cluster Thermodynamics Baradan Panta, Walker A. Glasoe, Juliana H. Zollner, Kimberly K. Carlson, and David R. Hanson* Chemistry Department, Augsburg College, 2211 Riverside Avenue South, Minneapolis, Minnesota 55454, United States S Supporting Information *

ABSTRACT: Particle formation and growth with H2SO4 molecules in an axially symmetric flow reactor was simulated with computational fluid dynamics. A warm (∼310 K) gas containing H2SO4 flows into a cooled section (296 K) that induces particle formation. The fluid dynamics gives flow fields, temperatures, and reactant and cluster distributions. Particle formation and growth are simulated with detailed H2SO4 cluster kinetics with chemistry based on measured small cluster thermodynamics and on bulk thermodynamics for large clusters. Results show that particle number densities have power law dependencies on sulfuric acid of ∼7, in accord with the thermodynamics of the cluster chemistry. The region where particle formation rates are largest has a temperature that is within 3 K of the wall. Additional simulations show that the H2SO4 concentration in this region is 5 to 10 times greater than the measured H2SO4: this information allows for direct comparisons of experiment and theory. Experiments where ammonia was added as a third nucleating species were simulated with a three-dimensional model. Ammonia was dispersed quickly and particle formation during this mixing was seen to be low. Downstream of the initial mixing region, however, ammonia greatly affected particle formation.

1. INTRODUCTION

within it. In addition, the model yields power dependencies on H2SO4 for the nucleation rate. Additional CFD simulations are presented: models of the H2SO4 detection region and a model of ammonia dispersal in the top half of the flow reactor. Detailed modeling of the H2SO4 detection region is important for comparison of the Zollner et al. results to theoretical and other experimental nucleation rates. Modeling the dispersal of ammonia is important for estimating its abundance where particles are formed. Also, simplified DCT models that include ammonia in the clusters were used to assess the amount of particle formation induced by ammonia.

Binary nucleation involving sulfuric acid and water is important for understanding particle formation in cold portions of the atmosphere, particularly the upper troposphere.1−6 Therefore, nucleation within the H2SO4/H2O system has been studied extensively in the laboratory,1−3,7−10 often with the use of flow reactors. Results2,3,7−10 vary widely, which hinders the understanding of this process and its application to the atmosphere. Experimental results are interpreted with flow reactor analysis models, such as well-mixed batch models, simple plug-flow analysis,11 fully developed laminar flow analysis,12 etc. Recently, Herrmann et al.13 presented a computational fluid dynamics (CFD) model simulation of their flow reactor experiments,10 and that effort led to an improved interpretation of these results. Here, CFD models (using FLUENT14) have been developed to understand the flow, temperature, and reactant and cluster distributions of the cylindrical flow reactor presented by Zollner et al.2 Since the main parts of the reactor2 are axially symmetric, the primary CFD model developed here is two-dimensional and axially symmetric. Particle production was simulated with detailed cluster thermodynamics (DCT), which is based on known and interpolated thermodynamics of H2SO4/H2O mixtures. The CFD model with DCT allows for the identification of the region, where the majority of particles are formed, which is called the nucleation zone (NZ), and the H2SO4 concentration, temperature, and typical residence times © 2012 American Chemical Society

2. MODEL DEVELOPMENT 2.1. Flow Reactor Geometry and Experimental Conditions. The flow, heat transfer, and H2SO4 clustering reactions inside the flow reactor were characterized using the two-dimensional (2D) capabilities of FLUENT in the primary model described here. The 2D model of the flow reactor2 is shown in Figure 1. It has a 5 cm ID and a total length of 157.5 cm; the flow inlet is at the top, and gas flows generally downward as it travels from warm (313 K) to cold regions (296 K). The model has four sections with walls held at different temperatures: the mixing region (MR) at 307 to 313 K, a Received: March 13, 2012 Revised: August 24, 2012 Published: September 17, 2012 10122

dx.doi.org/10.1021/jp302444y | J. Phys. Chem. A 2012, 116, 10122−10134

The Journal of Physical Chemistry A

Article

mostly N2 with up to 0.07 sLpm H2O vapor. The flow reactor has a total pressure of 0.97 atm. A trace amount of H2SO4 is present in the inlet flow, at a mixing ratio typically between 400 and 2000 pptv (pmol/mol.) The mixture passes first through the MR, then the transition region, with a length of 18 cm and wall temperature of 300 K, and then enters the nucleation region, which has a wall temperature of 296 K and a length of 100 cm. At the end of the reactor (wall at 300 K), a side port at X = ∼1.45 m allows for a particle counter inlet tube to sample the center of the flow (not included in the simulation.) At the end of this section, the model has the flow funnel down to a 1.2 cm ID exit tube ∼4 cm long. The actual H2SO4 detection region in the experiment,2 however, has a 7 cm long tube of 0.9 cm ID that leads to a ∼3 cm diameter sphere where a chemical ionization mass spectrometer (MS) monitors H2SO415 or Nbases.16 This latter geometry was used in additional simulations discussed below. 2.2. Cell Size and Thermal Properties. The primary 2D model has ∼32 000 cells, quadrilateral mesh faces, and maximum face diagonals (i.e., ‘face size’ in ANSYS) of 2 mm. The cells are mostly rectangular except near the exit funnel where some are skewed. Element thickness was chosen to be smallest near the wall where gradients in H2SO4 are largest. Elements next to walls are 0.1 mm thick by ∼1.5 mm long and each succeeding cell inward increases in thickness by 10% (i.e., geometric inflation of 1.1 is applied) for 28 cells; the full radius of the flow reactor is covered by an additional 20 square cells. Thermal properties of nitrogen and water vapor and their variation with temperature are detailed in Table 1.17 The effect Table 1. Thermal Parameters As a Function of Temperature and Water Content for N2/H2O Mixtures at a Total Pressure of 0.97 atma parameter

T (K)

Cp

290 300 600 284.05 323.15 290 320 296 313

Figure 1. Nucleation flow reactor of Zollner et al.2 as simulated here (top half and bottom half are side by side.) Inset at the bottom is a depiction of a separate three-dimensional (3D) simulation of the sulfuric acid detection region.

μ Kt Dcd

transition region at 300 K, a nucleation region that contains the NZ (wall at 296 K), and a particle sampling region at 300 K. In the experiment, two flows of nitrogen, one laden with sulfuric acid vapor and another with water vapor, enter the MR via small showerheads that enhance mixing. The warm temperature in this region suppresses particle formation. The showerheads were not modeled here, and the inlet has a 5 cm ID with gases entering well-mixed at a temperature of 313 K and a velocity profile close to plug flow. Note that gas entering the MR in the experiment was at room temperature, which was quite variable. Therefore, the inlet gas temperature was varied from 300 to 321 K in some runs (see Supporting Information for results.) The ammonia addition port in gray (X = 34 cm) and the mass spectrometer sampling region (dashed box in Figure 1) are modeled in separate simulations, described below. Sulfuric acid vapor was assumed to be lost on each collision with a surface and has a diffusion coefficient in the N2−H2O mixture, Dc, of ∼0.08 cm2/s that changes with relative humidity (RH) and temperature due to clustering with H2O.15 The total flow rate is 6 sLpm (STP, 273 K, 1 atm, L min−1), which is

value at 0.005 atmb value at 0.011 atmc 1048.17 1049.19 1079.89 1.703 × 10−5 1.878 × 10−5 0.0251887 0.0273720 8.064 × 10−6 9.492 × 10−6

1053.2 1054.3 1085.8 1.698 × 10−5 1.873 × 10−5 0.025148 0.027322 7.688 × 10−6 9.014 × 10−6

unit J/kg K

kg/m s W/m K m2/s

Specific heat Cp, viscosity μ, thermal conductivity Kt, diffusion coefficient, Dc. The temperature and water vapor dependent monomer diffusion coefficient includes hydration.15 bpH2O = 0.005 atm, 18% RH at 296 K. cpH2O = 0.011 atm, 40% RH at 296 K. dValues from Dc(T)/(m2/s) = A1 + A2T where Ai are functions of pH2O: for pH2O = 0.0044 atm, A1 = −1.68 × 10−5, A2 = 8.40 × 10−8 K−1; at pH2O = 0.0111 atm, A1 = −1.54 × 10−5, A2 = 7.80 × 10−8 K−1. a

of water vapor on the density and thermal properties of the N2−H2O mixture was calculated assuming ideal mixing. Thermal properties of sulfuric acid and its clusters are not important as they are present at ppbv (nmol/mol) levels or less. Upon a change in conditions, generally 800 iterations were found to yield converged cluster concentrations. Runs were made with two larger mesh sizes: one with ∼50% larger mesh edges that results in about 2/3 the number of cells (medium) and another mesh with cell sides ∼125% larger that resulted in about 65% less cells than the normal run (large). There were small, insignificant differences in H2SO4 concen10123



Article

2.3.3. Three-Dimensional Model with NH3 Addition. The ammonia addition experiments of Zollner et al.2 were simulated with a 3D model of a 88 cm length of the flow reactor: from the top of the MR to about the middle of the nucleation region. A side arm inlet (1/8″ OD Teflon tube with a 0.17 cm ID) introduced NH3 into the main flow, positioned at the wall, 32 cm from the top, with its flow directed inward. There were ∼138 000 volume elements with a maximum face size of 1.2 cm. A variable flow of nitrogen (0.025 to 0.25 sLpm) carried ammonia into the flow reactor such that a level of 16 pptv NH3 would be attained in the ∼6 sLpm main flow if no losses occurred. The diffusion coefficient for ammonia in N2 at 296 K is 0.22 atm cm2/s.18 A similar 3D simulation of ammonia dispersal into an N2 flow in a larger cylindrical reactor was presented by Hanson and Eisele.19 2.4. Cluster Kinetics and Chemistry. The detailed cluster thermodynamic (DCT) scheme used to simulate particle nucleation and growth consists of a series of forward and backward reactions involving sulfuric acid molecules and its cluster products (i.e., as in Lovejoy et al.33 and Friedlander et al.20) Particle formation starts with the clustering of two sulfuric acid monomers to form a sulfuric acid dimer. The dimer reacts with a monomer to form a sulfuric acid trimer and so on. Forward reactions were terminated at either the sulfuric acid dodecamer or the tetracosamer (24-mer). Clusters decompose via an accompanying reverse reaction: a dimer breaks down to two monomers, a trimer to a dimer and a monomer, a tetramer to a trimer and a monomer, and so on. Since clustering is terminated at a certain size, a buildup occurs, which is taken to be the concentration of new particles. The majority of the results reported here include reactions up to the dodecamer, with some runs performed up to the tetracosamer. Water molecules were not explicitly included in the DCT kinetics: this is the “quasi-equilibrium”20,21 approach. It is based on the assumption that each cluster rapidly equilibrates with the surrounding water vapor, which greatly simplifies the cluster calculations. Changes in RH affect the water content of H2SO4 clusters, which affects their size and their loss rate of H2SO4 molecules. The latter effect is dominant, and thus, the changes in the thermodynamics due to changes in RH were incorporated into the decomposition reactions. This is important as the model needs to simulate RH changes as temperature changes at a constant partial pressure of H2O. For example, the water content of the decamer changes with RH, and the thermodynamics were adjusted accordingly, which determined the rate of decamer decomposition. The change in size of a cluster with RH affects the kinetics, but the effect is small and was not included in the simulation: for a typical temperature change of 313 to 296 K, the RH changes from 6 to 16%, the collision rate with the monomer would increase ∼6%, and the decamer diffusion coefficient would decrease ∼10%. The thermodynamics, ΔGn0, for the dimer and trimer (n = 2 and 3, respectively) are taken from the measurements of Hanson and Lovejoy22 and for the nonamer and above are assumed to be equal to the bulk standard free energy of vaporization of the relevant sulfuric acid solution (calculated from Vermeulen et al.23 and Clegg et al.;24 both based on Giauque et al.25) The individual ΔH0 and ΔS0 between the trimer and the nonamer were linearly interpolated. A linear interpolation is simple and yields reasonable results; alternate interpolation schemes were explored, and these results are included in the Supporting Information.

trations ([monomer]) and axial velocity between the normal and medium mesh. For example, on-axis values of the [monomer] showed deviations of less than 1% from for X < 1.47 m and ∼5% different at X = 1.55 m (near the funnel at the end of the reactor.) The average deviation of on-axis axial velocities between the normal and medium meshes was 3 × 10−4 m/s (except at the funnel.) Differences were larger for the largest mesh: over the first 1.47 m of the reactor, [monomer] was within 5% of the normal run and average axial velocity deviation was 6 × 10−4 m/s. 2.2.1. Spatial Discretization. The first-order upwind procedure in FLUENT’s solver for scalars was used for the majority of the simulations presented here. The convective term is solved using concentrations at the faces of cells, which are interpolated from the cell center values.14 A set of simulations was also performed using second-order upwind discretization, and there was a small but significant effect on the monomer and cluster distributions: [monomer] increased about 8% in the NZ and [dodecamer] increased by about 50%. Further effects on monomer distributions are discussed in the Supporting Information. 2.3. Additional Simulations. Three additional geometries were developed for CFD: (i) the end portion of the primary 2D model was extended to include the MS detection region (dashed box in Figure 1), (ii) the MS detection region was also simulated with detailed 3D models, and (iii) a 3D model of the top half of the flow reactor was developed to simulate ammonia addition experiments. 2.3.1. Exit Tube in the Extended Two-Dimensional Model. Extension of the primary 2D model provides an estimate of [monomer] that is present in the MS detection region compared to that in the NZ, which is important for direct comparison to theory and experiment. A preliminary result was shown in Zollner et al., and it is also discussed here. Cluster chemistry was not simulated because clusters have extremely low abundances and do not impact the monomer or the ions. The bottom portion of the simulated flow reactor was modified: the ID of the small exit tube was decreased to 0.9 cm and lengthened to 13 cm, and it had a 3 cm diameter bulge centered 8.5 cm along its length. Note that the ion source and inlet were not included in the 3 cm bulge as they are 3D features (however, see section 2.3.2). Total flow was 6 sLpm or 4.5 sLpm; the latter was used to simulate monomer loss at X > 1.48 m because 1.5 sLpm is sampled by the particle detector at this position. There were 26 000 quadrilateral elements, maximum face size was 2 mm, and there were 18 inflation layers. 2.3.2. Three-Dimensional Model of MS Detection Region. A detailed 3D model of the MS detection region was developed to better model [H2SO4] in the exit tube and detection region. This model included a portion (∼20 cm) of the 2.5 cm radius flow reactor, the exit tube (ID of 0.9 cm), which is the inlet to the MS region, and an outlet tube of 0.6 cm ID. A cross-section of the geometry is shown in a contour plot in the Supporting Information. Flow is withdrawn by the particle detector at a rate of 1.5 sLpm at a position of 1.48 m (the flow rate downstream is reduced to 4.5 sLpm). The initial H2SO4 and velocity distributions were taken from the 2D simulation (section 2.3.1.) There were 6.7 × 105 tetrahedron elements, and face size maximum was 2 × 10−3 m. This and another, more detailed 3D model of the MS region are discussed in the Supporting Information. 10124



Article

Table 2. Values of ΔS0 (cal/mol·K) and ΔH0 (kcal/mol) As a Function of Cluster Size at Select Conditions

Standard thermodynamics and kinetic theory (i.e., as in Lovejoy et al.33) were used to obtain rate coefficients from the equilibrium constant Kc = RTKp ⎛ −ΔG 0 * ⎞ k1, n − 1 n ⎟⎟ = Kc*(n) = RTK p*(n) = RT exp⎜⎜ kn ⎝ RT ⎠

(1)

(H 2SO4 )n − 1* + H 2SO4 * ↔ (H 2SO4 )n *

(2)

∼6% RH, 69 wt %

pn = p(H 2SO4 )n * =

∑ p(H 2SO4 )n (H 2O)m m

K p*(n) =

ΔH0

ΔS0

ΔH0

ΔS0

2 3a 4 5 6 7 8 9b N>9

−18.3 −21.6 −22.4 −23.3 −24.1 −24.9 −25.8 −26.6 −26.6

−39.5 −45.5 −44.1 −42.6 −41.2 −39.8 −38.3 −36.9 −36.9

−18.4 −21.8 −23 −24 −25 −26 −27 −28 −28

−39.5 −45.5 −44 −42.5 −41 −39.5 −38 −36.5 −36.5

−18.5 −22.5 −23.9 −25.4 −26.8 −28.3 −29.7 −31.2 −31.2

−39.5 −45.5 −44.2 −42.8 −41.4 −40.1 −38.8 −37.4 −37.4

Values at ∼250 K from Hanson and Lovejoy22 with a small RH dependence for ΔH0. bBulk values assumed for nonamer and larger; 290−310 K vapor pressure data from Vermeulen et al.23 and Clegg et al.24 Both are based on Giauque et al.25 data, and enthalpy and entropy changes are within 0.8 kcal/mol and 0.8 cal/(mol K), respectively. Thermodynamics for n = 4 to 8 are linear interpolations. The effect of the hydration of the monomer on the bulk thermodynamics was not considered.

(3)

(4)

where pn is the pressure of the nth cluster and all its hydrates, and p1 is the pressure of the neat monomer plus its hydrates. The forward rate coefficients k1,n‑1 were taken from kinetic theory27 assuming unit reaction probability and the bulk density for 60 wt % sulfuric acid.17 For simplicity, the T1/2 temperature dependence and the change in cluster size with composition (due to RH changes, discussed above) were not included in k1,n‑1. The reverse rate coefficients kn are given by kn = A f T β exp( −Ea /RT )

ΔS0

a

pn pn − 1 p1

∼40% RH, 48 wt %

ΔH0

a

Written as functions of cluster number n, Kp is the standard state (atm−1) equilibrium constant and Kc is the molar (M−1) value with the gas constant R = 0.082 M−1 atm K−1. The forward rate coefficient for reaction 2 is denoted k1,n−1 in eq 1, and the decomposition reaction for the nth cluster is denoted kn. The * indicates that the water equilibrium assumption is in effect, and presented here for completeness, the relation to the actual hydrate distribution is:

∼16% RH, 60 wt %

n

approximation was used with the diffusion coefficients of clusters in N2 independent of temperature, while a temperature dependence of T1.75 was used for the monomer as it experiences a wide range of temperatures and because the hydration of the monomer changes significantly in the experiment.15 The latter two were simulated with a linear dependence of Dc on temperature, as shown in Table 1. For all conditions simulated here, the dominant sink for the monomer is wall loss because cluster concentrations are small. For most conditions simulated here, the primary sink for clusters up to the accumulating cluster is growth via monomer. Therefore, intercluster reactions (i.e., coagulation) were not included. The largest contributor to this would be reactions of clusters with the [dimer], which, for most simulations, is less than 0.001 times the [monomer]; reactions of the monomer dominate the cluster kinetics. 2.4.1. Monomer Hydration. The bulk thermodynamics give the vapor pressure of the unhydrated monomer (p1,0) over the bulk solution.26 Accounting for monomer hydration would have a significant effect on the thermodynamics of monomer evaporation rates from the large (n ≫ 3) clusters. Note that the effects of monomer, dimer, and trimer hydration were included in the thermodynamic measurements of these clusters.22 Incorporating the effect of monomer hydration could be achieved by making the large clusters less stable: i.e., to correct for assuming that p1 = p1,0 in eqs 3 and 4. A typical hydration at 16% RH gives an equilibrium hydration level that 2 out of 3 monomers have one or more H2O ligands: p1/p1,0 = 3.15 This would translate into an effective change of +0.6 kcal/mol in ΔG0. There is a difference of ∼0.6 kcal/mol in ΔG0 between the two23,24 bulk treatments. Note also that vapor pressure measurements of neat liquid sulfuric acid34,35 lead to a difference of a factor of 2.3 at 296 K, which is about the size of the effect of hydration. Monomer hydration was not taken into account because its effect is roughly the size of these uncertainties.

(5)

where Ea is the activation energy, which, at constant, RH can be taken equal to −ΔH0, and the product of the pre-exponential factors Af and Tβ can be obtained from eq 1, the k1,n‑,1 and the ΔS0 values. The experimental conditions focused on here are 16 and 40% RH at 296 K, which are partial pressures of H2O of 0.0045 and 0.011 atm, respectively. Since the gas temperature is as high as 313 K, thermodynamics were also needed for clusters at 6% RH conditions. These three RH conditions of 6.1, 16 and 40% give equilibrium compositions of bulk H2SO4 solutions of 69, 60, and 48 wt % H2SO4, respectively. Listed in Table 2 are the ΔH0 and ΔS0 for the clusters in kcal/mol and cal/(mol K) for 6.1, 16 and 40% RH for temperatures near 296 K. Since the simulation and experiment are not conducted at constant RH, Ea larger than |ΔH0| for clusters n ≥ 4 were used. Also, since Af is limited to a value of 1038 M−1 s−1 K−β, a β = 2 was used. Table 3 lists the cluster/particle diameter (Dp), their diffusion coefficients Dc, the forward and reverse reaction rates for the 16% RH thermodynamics at 296 K, and the values of Ea and Af for β = 2. Note that these values simulate the change in cluster thermodynamics due to their change in composition with RH in the experiment. For example, the Ea and Af were adjusted for each cluster so that the reverse rates closely matched those given by the thermodynamics in Table 2 for 16% RH at 296 K and for 6% RH at 313 K. The diffusion coefficient and diameter for each cluster are based on 60 wt % H2SO4 bulk densities,17 and the diffusion coefficients in N2 for dimer and greater are based on kinetic theory.27 Since water vapor is at most 1.2% of the mixture, full multicomponent diffusion was not simulated. The dilute

3. RESULTS AND DISCUSSION 3.1. Temperature Profile and Flow Visualization. Temperature changes along the flow reactor have a significant 10125



Article

Table 3. Diameter, Diffusion Coefficients, Forward Rate Coefficients, Inverse of 296 K Reverse Reaction Rate (Lifetime of Cluster), and Ea and Af for krev for Hydrated H2SO4 Clusters at pH2O = 0.0044 atm (Equivalent to 16% RH at 296 K) n 1 2 3 4 5 6 7 8 9 10 11 12

Dp (nm)a 0.690 0.869 0.995 1.095 1.179 1.253 1.319 1.379 1.434 1.485 1.533 1.578

Dc(cm2/s)b 0.08−0.095 0.0582 0.0445 0.0367 0.0316 0.0280 0.0253 0.0231 0.0214 0.0199 0.0187 0.0177

kfor (M−1s−1)c

1/krev(s)d

Eae

× × × × × × × × × × ×

7.19 × 10−6 1.71 × 10−4 1.81 × 10−3 1.94 × 10−2 2.10 × 10−1 2.30E+00 2.53 × 101 2.80 × 102 2.68 × 102 2.56 × 102 2.46 × 102

76.6 92.0 120.7 149.4 178.0 206.7 235.4 255 255 255 255

Af e

f

2.56 2.83 3.11 3.38 3.63 3.87 4.10 4.31 4.52 4.71 4.90

1011 1011 1011 1011 1011 1011 1011 1011 1011 1011 1011

5.16 1.17 1.28 1.40 1.50 1.60 1.69 4.58 4.79 5.00 5.19

× × × × × × × × × × ×

1013 1015 1019 1023 1027 1031 1035 1037 1037 1037 1037

Assuming 3.5 water for each H2SO4, i.e., about 60 wt % bulk density (CRC). bDiffusing in 0.97 atm total pressure (N2 with ∼0.5% H2O vapor.) Rate coefficient to form nth cluster, from monomer plus (n − 1)th cluster. Rate from kinetic theory as presented in McMurry.27 dLifetime of cluster against loss of H2SO4 at 296 K. eβ = 2 to simulate change in thermodynamics with RH. Ea in kJ/mol; Af in M−1 s−1 K−2. fVaries with temperature due to changes in hydration;15 see Table 1. a c

Figure 2. Contour plots from 2D simulations with vertical coordinate, radius, multiplied by 6.7. Flow inlet at the left, gravity pulls to the right, the reactor axis is on the bottom of the contour, and the wall is at the top. Particle sampling occurs just before the funnel at the right, and the flow exiting the small diameter port enters the mass spectrometer detection region. (a) Temperature and (b) axial velocity contour plots.

Figure 3. Pathlines originating at 1 mm intervals radially (Y = 0 is the flow reactor axis) at X = 0.44 m showing the dominance of recirculation in this region and at X = 0.7 m showing that particles formed in this region flow directly to the particle counter at X = 1.45 m.

effect on flow patterns and on particle formation. Contours of temperature (a) and axial velocity (b) are shown in Figure 2

(radial coordinate scaled by a factor of 6.7 compared to the axial coordinate.) Gas from the warm mixing region cools 10126



Article

Figure 4. Set of contour plots for 1200 pptv initial H2SO4. The Y axis (radial coordinate) has been amplified by a factor of 5; 2.5 cm is the radial extent of the reactor. Axial position (X axis) is demarcated every 20 cm. The plots are in mass fraction from top to bottom: monomer, hexamer, decamer, undecamer, and dodecamer. Mole fraction can be obtained by dividing mass fraction values by ∼3.5, 21, 35, 38.5, and 42 for the monomer through dodecamer, respectively.

Within the negative velocity zone, the flow is primarily recirculating and none of the streamlines exited this region in ∼90 s of simulation time, while streamlines originating in the NZ flow toward the exit of the flow reactor and arrive at the particle detector in about 12 s. More flow streamlines and particle paths generated with FLUENT’s postprocessing capabilities are shown in the Supporting Information. 3.2. Nucleation Zone and Particle Formation. The DCT within the primary 2D simulation provides particle formation rates and particle/cluster abundances for comparison with experimental results. The simulation also provides information on how many regions other than the NZ contribute to particle abundances. Discussion of these starts with cluster abundances, and cluster mass fractions are presented in Figure 4 as 2D contour plots for the monomer, hexamer, decamer, undecamer, and dodecamer. The radial coordinate has been scaled by a factor of 5 compared to the axial coordinate for presentation. As it is carried down the reactor, monomer abundance steadily decreases with axial distance X due to diffusion and loss on the wall with a 90% or greater loss by the end of the reactor. As temperature decreases, the hexamer, decamer, and undecamer become more stable, and their abundances peak in the region just downstream of the buoyant zone. These clusters are lost on the wall (slowly) and grow to larger clusters (rapidly) as they travel further down the reactor. The dodecamer abundance ([dodecamer]) is the largest of any cluster because it is the accumulation cluster: its abundance is

before significant particle nucleation occurs, and this cooling takes place near the wall. Gas density increases, and the flow is large near the wall at X ≈ 25 cm and from X = 60 to 120 cm, the flow converges toward the reactor axis. Gas temperature stays above ∼300 K for some distance into the nucleation region. H2SO4 diffuses and is lost on the walls, and therefore, [monomer] decreases toward the walls and down the reactor. As a result, [monomer] is largest on the axis, and nucleation rates peak at the center of the reactor at about X = 70 cm, ∼30 cm into the 296 K nucleation region. Above this region, a warm, buoyant section develops centered on the reactor axis. Gas in the buoyant zone generally rises along the axis and warms to ∼310 K before turning radially, toward the wall. A rough estimate of the amount of mass flow upward is about 0.2 sLpm based on a velocity of ∼1 cm/s over a roughly 1 cm radius area, about 3% of the total flow rate. Much of this flow comprises gas that is recirculating, shown with streamlines in Figure 3 and discussed below. It is apparent in Figure 2b that negative velocities are somewhat larger in the transition region than in the nucleation region; probably due to a larger temperature difference in the walls of the mixing and transition regions than that between the transition and nucleation regions. The temperature differences of the mixing, transition, and nucleation region walls affect the extent of the buoyant zone; see the Supporting Information for simulations with variable mixing region temperatures. Streamlines in Figure 3a originate in the negative velocity zone at X = 0.44 m and (b) originate in the NZ at X = 0.7 m. 10127



Article

Figure 5. (a) Nucleation rate and axial velocity along the flow reactor axis, Y = r = 0. Left axis for microscopic nucleation rates: the net flux between n = 7 and 8 (J8) and between n = 11 and 12 (J12). Axial velocity (solid curve) plotted with the right axis. (b) Number density of dodecamer (solid curve, right axis), its decomposition rate (dash-dot), and J12 (dash) along the flow reactor axis. The nucleation zone is depicted with the dotted box. (c) Monomer mixing ratio (solid) and temperature (dash) along the flow reactor axis.

taken to be that of the particles, Np. After [dodecamer] peaks, it drops slowly via diffusion and loss on the wall, and it is spread out axially compared to the others. Some dilution of [dodecamer] occurs as X increases: centerline flow with high [dodecamer] mixes with flow containing lower [dodecamer] that originated off axis. Microscopic nucleation rates were obtained from the model by calculating the net flux between adjoining clusters: i.e., Jn = formation rate of n from n − 1 plus monomer minus the decomposition rate of n.20 Presented in Figure 5a are two

different Jn as well as the axial velocity at Y = 0 as a function of X, the distance down the reactor (left axis, nucleation rate; right axis, axial velocity). J8 represents the flux between the heptamer and octamer: J8 = [heptamer][monomer]k1,7 − [octamer]k8, where [nmer] is the concentration of the nth cluster. Similarly, J12 is the flux between the undecamer and dodecamer. Over the length of the flow reactor, J8 and J12 are in good agreement. Nucleation rates peak at X = 70 cm with a value of ∼0.35 cm−3 s−1 for these conditions. 10128



Article

particle detector. Therefore, particle tracks were generated to estimate the flux of particles that can escape the BZ. Inert, 6 nm diameter particles allowed to undergo Brownian motion but not evaporate were released in the BZ and in the NZ. Ten particles were released at points evenly spaced along Y = 0 to ∼1 cm radial distance from the axis, and they were tracked for up to 300 s. Plots of their X and Y positions versus time are shown in the Supporting Information. All particles released within the BZ were carried down to the particle detector where they had an average radial position of 0.07 cm (maximum of 0.13 cm). Particles released in the NZ reached the detector much sooner and had an average radial position of 0.4 cm (max position of 0.8 cm, particles were confined to a region of ∼1.6 cm diameter). Roughly, the flux of particles from a given source region is proportional to the relative area that they cover in the target region. With about 1/2 the Np in the buoyant zone as in the NZ (see Figure 5b), the spatially averaged contribution of BZ particles is roughly 1.5% (i.e., 1/2 times the square of the ratio of average radii, 0.07 cm/0.4 cm) of the amount of NZ particles. This probably does not fully represent the relative flux of particles from the two zones as the BZ particle’s median arrival time is about 15 times that of the NZ released particles. Assuming the flux of particles is also inversely proportional to median arrival time suggests that this 1.5% potential contribution should be further divided by 15. This analysis suggests there is a potential 0.1% contribution of BZ particles to the total that are measured. 3.3. Power Dependence on Sulfuric Acid. The [dodecamer] = Np was calculated for initial sulfuric acid concentrations ranging from 400 to 1600 pptv. Figure 6 shows the resulting macroscopic nucleation rate Jp as a function of the H2SO4 concentration (monomer) at the peak value of the microscopic nucleation rate J12 (X = 0.7 m, Y = 0 m); [monomer] ranged from 4.5 × 109 cm−3 to 1.8 × 1010 cm−3

A nucleation zone, NZ, can be nominally defined as the region where J is within 50% of peak J, and it has an on-axis length of ∼48 cm. However, the 50% cutoff upstream of peak nucleation falls inside the negative velocity region. Although it has significant nucleation occurring in its lower portion, particles formed in the negative velocity zone do not contribute significantly to those observed downstream (see section 3.2.1). Therefore, the upstream side of the NZ is taken to be 57 cm, the point on axis where the velocity is zero. The NZ with a length of ∼40 cm is shown in Figure 5b, which is a plot of the Y = 0 (on axis) values of [dodecamer], its decomposition rate, and J12. Although the [dodecamer] (solid line) is substantial in this negative velocity region, its decomposition is, too: [dodecamer] decreases as the flow warms (see temperature profile in Figure 5c). Also, when J12 (dashed line) is positive, recirculating flow dilutes [dodecamer]. At the warmest part of this region, the high dodecamer decomposition rate results in a bump in the [undecamer] as seen in Figure 4 (note that this bump in [undecamer] is an artifact due to the truncation of clusters at the dodecamer.) Clusters evaporate at the top of the buoyant zone and do not contribute significantly to the observed particles downstream (see more discussion of this below.) Hence, the starting point for particle production in the experiment begins near 57 cm, or perhaps even 63 cm when axial velocity becomes greater than ∼1 cm/s. Following the flow downward through the NZ, [dodecamer] is boosted via particle formation until about 90 cm, where it declines as dilution and diffusion losses overtake nucleation. The large decrease that occurs after 1.5 m is due to disturbances in the flow entering the warmed sample region and could in part be artifact due to meshing or other difficulties associated with the converging geometry at the end of the reactor. Figure 5c is a plot of temperature and [monomer] at Y = 0 as a function of X. In the NZ, the temperature ranges from 300.5 to ∼297 K, and the monomer mole fraction drops from 0.6 to 0.4 ppbv from an initial value of 1.2 ppbv. At the peak of J12, X = 70 cm, T = 299 K, and [monomer] = 0.5 ppbv. A macroscopic nucleation rate Jp can be defined as Np/tnuc where Np is particle number density and tnuc is an average residence time in the NZ. Np is equal to [dodecamer] near the center of the reactor at X = 1.45 m, and tnuc is estimated to be 10 s. This time was estimated by considering the on-axis NZ parameters: it is ∼35 cm long and has an average axial velocity of ∼3 cm/s. Note that [dodecamer] in low velocity regions such as the very top 5 cm of the NZ might have a very long residence time. But this region’s contribution to Np is small at Y = 0, X = 1.45 because the flux of gas through it is a very small portion of the total flow. With the [dodecamer] at about 2 cm−3 from Figure 4, Jp = Np/10 = 0.2 cm−3 s−1. Note that the macroscopic Jp obtained in this manner is very close to the average value of the microscopic J12 in the NZ shown Figure 5b. 3.2.1. Particles in the Buoyant Zone. The DCT predicts significant nucleation rates within the buoyant zone (BZ) or negative velocity region where temperatures and H2SO4 are elevated. While the simulation shows that newly nucleated particles in the BZ evaporate as they rise and warm, the DCT cannot give direct evaporation information on clusters containing more than 12 (or 24) H2SO4 molecules. Small particles in the recirculating pathlines in Figure 3a can, via diffusion or inertia, enter streamlines that are sampled by the

Figure 6. Nucleation rate (macroscopic, Jp) as a function of [monomer] at the peak of J12: Y = 0 m, X = 0.7 m. Experimental data2 for 14 and 20% RH are scattered and fall within the outlines shown. The simulation is the solid squares; other predictions (triangles29 and diamonds30) are also shown. Power dependencies of J on [H2SO4] are indicated next to the respective data. 10129



Article

Figure 7. Results with the DCT with the accumulating cluster at n = 24. Mass fraction contour plots as in Figure 4 for monomer, hexamer, and decamer; mole fractions can be obtained from mass fraction for the last two clusters by dividing 80.5 for the tricosamer and 84 for the tetracosamer.

= 12 with k1,n‑1 increasing, and Dc decreasing, with cluster size beyond this. In Figure 7 are contour plots of the mass fraction of monomer, hexamer, decamer, tricosamer, and tetracosamer with initial [H2SO4] of 1200 pptv. The mole fraction of the tetracosamer in the particle sampling region was slightly larger than that obtained for the dodecamer DCT, about 20% greater. This is likely due to diffusion coefficients decreasing with cluster size, and more of the accumulating cluster survives. Because the DCT truncates the cluster size, diffusional loss is enhanced in the model and thus Np is likely to be underpredicted. Nonetheless, the simulations show that particles readily grow to 2 nm in diameter; a size at which the Zollner et al. particle detectors are fairly efficient. The ratio of [tetracosamer]/[tricosamer] is 310, near the particle sampling region (X = 1.5 m, Y = 0 m) for 1200 pptv initial H2SO4. For the primary simulation with the dodecamer as the accumulating cluster, the [dodecamer]/[undecamer] ratio is 315. These high ratios indicate plenty of potential for particle growth beyond the accumulating cluster for these conditions. An additional test for particle size was performed with model runs using thermodynamics for 40% RH at 296 K. Clusters up to n = 24 were included. At a very low initial sulfuric acid concentration, 200 pptv, the [tetracosamer] was about 15 times larger than the [tricosamer] at X = 1.45 m indicating that there is still potential for growth beyond the tetracosamer even under low growth rate conditions at low [H2SO4] (about 2/3 of the lowest used in the Zollner et al. experiment.) Included in the Supporting Information are the 40% RH cluster thermodynamics, contour plots of selected clusters, and comparisons to the 16% RH results.

(180 to 720 pptv). For the lowest [monomer], the nucleation rate is about 10−4 cm−3 s−1, and as [monomer] increases, Np increases until Jp reaches ∼1 cm−3 s1 for a [monomer] of 1.8 × 1010 cm−3. The power dependence of Jp on sulfuric acid is 6.8 thus the DCT indicates that the critical cluster1 contains 6 to 7 sulfuric acid molecules. The nucleation theorem1,28 states that the size of the critical cluster is near the maximum in the cumulative free energy change (ΔGcum*), which occurs at n = 6 using the thermodynamics in column two of Table 2 for 299 K and the pH2SO4 listed above (a plot of a few ΔGcum* versus n is shown in the Supporting Information, Figure A4b). The DCT gives nucleation results that are in accord with the nucleation theorem. The Zollner et al.2 measurements at 14 and 20% RH are quite variable and are encompassed by the polygons in Figure 6. The simulation shows a decent congruence with these data especially the power dependence: a typical power dependence of 6 for these results is also shown. Note that the peak of the nucleation rate in Figure 5c occurs at a temperature of about 298−299 K. Discussed below are possible reasons for the deviation of the DCT predictions and the experimental data. Predictions from theoretical treatments29,30 of H2SO4/H2O homogeneous nucleation at 298 K and 16% RH are also shown in the figure, and they contrast strongly with the DCT results and the experimental2 data. The predictions are based on the classical liquid drop model and both have power dependencies of ∼16 for these conditions. 3.4. Particle Size. When the DCT model was extended to the tetracosamer, the number of reactions roughly doubled to 46, and the dodecamer (Dp ≈ 1.6 nm) was replaced by the tetracosamer (Dp ≈ 2 nm) as the accumulating cluster. The same thermodynamic and kinetic parameters were used up to n 10130



Article

3.4.1. Cluster Size Distributions. A reviewer pointed out that a detailed cluster treatment should in principle yield cluster size distributions at steady state. However, clusters were terminated at a certain size (the accumulation cluster) due to computational constraints: FLUENT caps the number of species at 50. A few runs were made at low [H2SO4] to see if the concentration of the accumulating cluster (i.e., tetracosamer) might become representative of steady state (i.e., not too much greater than concentrations of smaller ones) but cluster concentrations became so small that the CFD solver diverged. The effect of truncating the clusters at a certain size can be demonstrated with the simulation data presented above. Cluster distributions from the 2D primary model at two different points (X = 0.7 and 1.45 m, both Y = 0 m) are shown in Figure 8 for the dodecamer (filled symbols) and the

significantly affected by other losses, such as diffusional loss to the wall. 3.5. Ammonia Addition. The 3D model simulation of the top half of the flow reactor shows the dispersal of ammonia vapor into the reactor and consequent reactions with H2SO4 molecules. Trace amounts of NH3 in N2 were entrained into the main flow through a small inlet tube. The inlet tube of 1.7 mm ID was positioned with its outlet flush to the flow reactor wall, pointing toward the axis. The distribution of NH3 in the main flow depends on the flow rate of nitrogen in the inlet tube. Mixing via diffusion is fast, and wall loss is substantial. Dc for NH3 is 0.22 cm2/s at 296 K.18 A set of NH3 contour plots for a range of N2 flows are shown in the Supporting Information, and they show that, for N2 flows of 0.03 to 0.10 sLpm, ammonia concentrations are centered on the axis, which allows for an effective interaction between NH3 and H2SO4. For an N2 flow of 0.05 sLpm, on-axis NH3 mixing ratios are 15 to 20 pptv about 20 cm downstream of the addition point, approximately the value if evenly distributed and without loss. However, there is ammonia loss which is reflected in decreased [NH3] near the wall and down the reactor. Furthermore, when a strong interaction between NH3 and H2SO4 was assumed, there was additional loss of gas-phase NH3. The base addition protocol of Zollner et al. (i.e., an inlet N2 flow rate of ∼0.05 sLpm) as modeled here shows that base abundances in the NZ are within ∼30% of the no-loss expected value. To evaluate potential particle production in the region where ammonia initially enters the main flow, a DCT scheme involving clusters with one NH3 and up to 12 H2SO4 molecules was implemented. The forward rates for ammonia addition to the neat H2SO4 clusters were taken to be gas-kinetic. The loss rate of NH3 from H2SO4·NH3 and (H2SO4)2·NH3 was 20 and 0.002 s−1, respectively. Ammonia was assumed to be stable in the n = 3 clusters and higher: loss rate of ammonia from these clusters was 0 s−1. The loss rates of H2SO4 from the ammoniated clusters were taken to be 0.001 times the rate of the corresponding neat cluster. The abundances of the ammonia containing clusters peak in the center of the flow reactor and 20−30 cm downstream of the addition point. Cluster abundances are not elevated in or near where the ammonia addition plume first contacts the main flow. The high ammonia levels in the plume do not significantly enhance cluster formation because of depletion of H2SO4 where the flows mix. Contour plots and details of two NH3DCT thermodynamic schemes are presented in the Supporting Information. Cluster distributions on axis about 8 cm downstream of the NH3 addition point are shown in Figure 9. The nonammoniated clusters for n ≥ 3 (open symbols) are far less abundant than the clusters containing one NH3 molecule (filled symbols); indicating a potentially large effect on nucleation due to ammonia. The enhancement for the n = 12 cluster is about a factor of 108 (i.e., [(H2SO4)12·NH3]/[(H2SO4)12]), which is higher than the factor of 5 × 105 (corrected for RH difference, 27 vs 16%) reported by Zollner et al. for comparable conditions (∼20 pptv NH3 and 3 × 109 cm−3 H2SO4 at X = 0.7 m.) The simulated Np of ∼5 × 104 cm−3 compares more favorably with the measured number density of 1 × 103 cm−3. The over prediction by the model indicates that the stabilities of the ammoniated clusters are probably overestimated in the crude DCT used here.

Figure 8. Semilog plot of the on-axis (Y = 0 m) simulated cluster distributions at two axial positions: X = 0.7 (blue) and 1.45 m (red.) Results are shown using the 2D primary simulation with two DCT models: either truncated at the dodecamer (filled triangles and diamonds) or at the tetracosamer (squares, circles).

tetracosamer (open symbols) as accumulating clusters. The cluster abundances for n = 2 to 10 differ by less than 0.1% between the two DCT models; the undecamer is ∼10% different, and the dodecamers differ by a factor of about ten. The very small affect for clusters that are 2 or more molecules away from the accumulating cluster indicates that microscopic nucleation rates are reliable for these sizes (even J12 is not greatly affected in the NZ, see Figure 5a.) Note also that the agreement of the macroscopic nucleation rates for these two DCTs indicate that effects due to truncation are small or negligible. Finally, the shape of the distributions at X = 0.7 and X = 1.45 m differ, indicating a change in the loss process for the clusters. The general decrease in cluster concentration with size at X = 0.7 m reflects that growth by the addition of monomer is the dominant loss process. However, the shape of the distribution at X = 1.45 m indicates that cluster concentrations are 10131



Article

of axial distance at the funnel constriction, from X = 1.51 to 1.53 m, well above the ∼1.5% loss per cm in the main portion of the flow reactor. A few different meshes of the exit tube assembly were generated. The mesh was adapted (reduced) in and near the funnel such that the face size was 1/4 of that previously and loss decreased to ∼2% per cm at 1.52 m but [monomer] in this region (and downstream) did not converge (the solver even diverged in some runs). A triangular mesh was implemented, and monomer distribution in the exit tube was better behaved, but the axial loss was still large, 7% per cm at 1.52 m. Note that the 2D primary model with an exit tube of 0.6 cm diameter exhibited much more gradual and steady decreases with axial distance. The second order upwind solver scheme was selected for some runs, and it had a significant effect on simulated monomer distributions in this region (see the Supporting Information.) The geometry of the exit tube assembly when the flow constricts from a tube radius of 2.5 to 0.45 cm requires additional meshing efforts in 2D as well as exploration of different solver schemes. Also, the more reasonable performance of a triangular mesh indicates that tetrahedral meshes in 3D may help better catch mass transport in this portion of the flow reactor. These simulations are presented in the Supporting Information along with the procedure for calculating f NZ. The estimates of f NZ from the various simulations range from 4.3 to 10.7. The variability due to the different meshes and solver schemes is particularly large near the exit tube in both the 2D and 3D simulations. Despite the attempts to build a reliable simulation of this region of the flow reactor, the monomer distribution and thus the uncertainty in f NZ remain large. Therefore, the present work does not improve upon the f NZ reported by Zollner et al.: an f NZ of 7 with an uncertainty of factor of 1.4 (i.e., ×1.4 and ÷1.4). A factor of 7 was applied to the measured H2SO4 so that the experimental data can be plotted in Figure 6. Note that this does not include the uncertainty in the measured value of [H2SO4], likely to be on the order of 25−50% (see Appendix A2 of ref 19).

Figure 9. Cluster distributions (X = 0.4 m, Y = 0 m) given by the 3D simulations of NH3 addition to the flow reactor (this DCT is detailed in the Supporting Information.) Three distributions as a function of nH2SO4 molecules are shown: filled squares are those containing one NH3 molecule; open squares and + symbols have no NH3. The + symbols are a control run in the absence of NH3, and neat clusters are slightly higher than the neat clusters when NH3 was present.

The + symbols are simulated neat clusters in the absence of ammonia: they are nearly identical to the neat clusters in the presence of ammonia. Note that the neat clusters in the 3D simulation are far different than those of the primary 2D model. For the same initial monomer, the primary model gives a [monomer] at X = 0.7 that is about a factor of 2.5 times larger, and [dodecamer] is about 3000 times larger. This is probably due to numerical diffusion in the 3D simulations that leads to a large loss of monomer. Note that this deficiency in the 3D model for nonammoniated clusters does not destroy the comparison of simulated to measured ammoniated clusters because the [monomer] at X = 0.7 m is close to that estimated for the Zollner et al. ammonia experiments. 3.6. Inclusion of MS Detection Region. There is large H2SO4 loss between the NZ and the MS detection region. The CFD simulations can be combined to yield a factor (f NZ) that can be applied to the measured H2SO4 to yield the value in the NZ. This is similar to the wall loss factor previously introduced9,10 with the main difference being the recognition that peak H2SO4 in the MS region differs significantly from the measured value.2,7 Ions drift perpendicularly across the flow in the MS region, reacting with a range of [H2SO4] between the ion source and the MS inlet, and product ions are representative of the average [H2SO4] in this region. Thus, f NZ is the ratio of the on-axis H2SO4 at peak nucleation to the measured H2SO4. Note that [H2SO4] varies significantly (e.g., 400 to 600 pptv) along the NZ, which is not taken into account here. For the 2D extended simulation, large axial gradients in monomer distributions in the exit tube assembly were evident (see Supporting Information) at points where the geometry changed drastically. For example, there was a 22% loss in 2 cm

4. COMPARISON TO PREVIOUS WORK A 2D approach was used by Herrmann et al.13 to model a nucleation experiment10 that is similar to Zollner er al.2 For example, the Herrmann reactor has a 6 cm ID and is about 2 m in length, and experimental conditions are similar: a warm flow of air at 6 to 10 sLpm is cooled to ∼298 K, total pressure is 1 atm, and H2SO4 is present in the ppbv range. Herrman et al., however, simulated particle formation and growth with a fine particle model (FPM)31 add-on code to FLUENT. Another important difference with the present study is that Hermann et al. limited their model to the final section of the experiment, which was held at 298 K, the nucleation region. They assumed that the gas from the ∼320 K transition section in the experiment10 had uniformly cooled to 302 K at the entrance to the nucleation region; they also assumed a plug flow velocity profile at that point. This is not in agreement with the present simulations, where warm temperatures persist for some distance into the next section of the flow reactor, and velocity profiles upon entering a cooler section are far from plug flow. Note that the 302 K assumption was based on a temperature measurement, but this measurement was admittedly inaccurate.10 Nonetheless, the boundaries of the NZ and the profiles of particles are qualitatively similar between the two studies reflecting the similarity of the experimental conditions. There 10132



Article

this discrepancy: (i) the presence of an impurity in the experimental work, which could be below detection,2 (ii) uncertainties in the measured [H2SO4] of ∼30% and in the simulated ratio (f NZ) of NZ to measured [H2SO4] of ∼40%, and (iii) the temperature of the gas and the inside wall of the mixing region was measured to be ∼6 K less than the 313 K set point. This latter effect was shown to be significant (see the Supporting Information): a lower temperature leads to an enhancement in particle number of about a factor of 4 for a 6 K cooling. Future development of the simulations of the mass spectrometer detection region can help to reduce the ∼1.4 factor uncertainty in f NZ. Variations in experimental conditions such as the sampling, ion source flow rates, and the size of the ion source opening can also provide information that will help reduce uncertainty in f NZ. The DCT can be altered to explore the sensitivity of Np to the chosen thermodynamics, i.e., setting the ΔG0 equal to the bulk value at clusters much larger than n = 9. For the addition of base molecules, it can be further developed to explore the addition of more than one base molecule and also to investigate species more basic than ammonia, e.g., amines.

are significant differences in the monomer distributions between these studies, which affects the interpretation of the respective experimental results. The Herrmann et al. simulation better agreed with measurements (particle sizes and H2SO4 wall loss calculations) when a value of 0.06 cm2/s was used for the H2SO4 diffusion coefficient; considerably less than the measured values15 and the values used here (Table 1). Consequently, their simulated wall losses were much less than those presented here. This led to a smaller calculated value for [H2SO4] in the NZ than our simulation would give because the measured [H2SO4] at the end of the reactor is the anchor point, similar to the procedure used here. Lower [H2SO4] results in smaller calculated growth rates for the particles: the Hermann et al. simulates final particle sizes well below 5 nm in diameter. It is possible that their simulation conditions compensated for an assumed low value for the diffusion coefficient. Also, later measurements in the same apparatus of the H2SO4 loss factor32 do not agree with the earlier values.10 Another explanation is that potential impurities in the experiment10,32 led to higher rates of nucleation than can be explained by H2SO4−H2O homogeneous nucleation rate calculations. Note that the Brus et al.10,32 experimental nucleation rates do not agree with those of Zollner et al. Finally, quantitative comparison of the DCT versus the FPM was not done here because the FPM is not freely available. However, the FPM is based on the nucleation treatment of Vehkamaki et al.;29 thus, quantitative agreement between FPM and the current DCT are not expected.

■

ASSOCIATED CONTENT

S Supporting Information *

Model and experimental results for temperature variations of inlet gas and mixing region wall. Also included are model results for alternate DCT schemes, 40% RH conditions with the tetracosamer DCT, particle tracks, monomer distributions from the various models, and 3D NH3 addition and clustering kinetics. Finally, future development of the model is discussed. This material is available free of charge via the Internet at http://pubs.acs.org.

5. IMPLICATIONS FOR EXPERIMENTAL RESULTS AND FUTURE CFD WORK The CFD simulations provide the conditions in the NZ and the approximate nucleation times for recent experiments.2 Particle growth information from the simulation indicate that they were large enough to be efficiently detected. A 3D simulation of the addition of NH3 via a small side port shows the details of the initial mixing of the base molecule into the flow and the consequent cluster formation. Thermodynamics of H2SO4 clusters within the water equilibrium assumption were constrained by previous experimental data for small clusters and bulk liquids. With detailed cluster thermodynamics determining the kinetics, the CFD simulations give results that are in satisfactory agreement with experiment for absolute particle abundances and in good agreement with the experimental power dependencies on [H2SO4]. For the present DCT, truncating the cluster calculations at the dodecamer does not lead to artifacts in cluster distributions for the decamer and smaller. If truncation in the DCT occurs sufficiently far from the critical cluster, the DCT is sufficiently accurate. Nucleation rates calculated with the classical liquid drop model have power dependencies on sulfuric acid of ∼16, far from the observed and the DCT based calculations of 6 to 7. These latter two results suggest there are about 6−7 H2SO4 molecules in the critical cluster for these conditions. For the RH conditions investigated, the power dependency of J on RH is about 6, which agrees with the experimental power dependency.2 Note that power dependencies on RH are not equal to the number of water molecules in the critical cluster, please see Zollner et al.2 for further discussion. Absolute rates from the DCT are somewhat lower than the experimental rates, about a factor of 15 different at an [H2SO4] of 7 × 109 cm−3. There are a number of potential reasons for

■

AUTHOR INFORMATION

Corresponding Author

*Phone: (612)330-1620. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by NSF Grants 0943721, 1068201, and 1051396. Augsburg College’s URGO program is acknowledged for funding travel and summer support. We gratefully acknowledge D. Oberreit and D. Porter for help with ANSYS and FLUENT and the University of Minnesota Supercomputing Institute for support. We also are grateful for the helpful comments of an anonymous reviewer.

■

REFERENCES

(1) Laaksonen, A.; Talanquer, V.; Oxtoby, D. Annu. Rev. Phys. Chem. 1995, 46, 489−524. (2) Zollner, J. H.; Glasoe, W. A.; Panta, B.; Carlson, K. K.; McMurry, P. H.; Hanson, D. R. Atmos. Chem. Phys. 2012, 12, 4399−4411. (3) Kirkby, J.; et al. Nature 2011, 476, 429−433. (4) Kazil, J.; Lovejoy, E. R.; Jensen, R. J.; Hanson, D. R. Atmos. Chem. Phys. 2007, 7. (5) Clarke, A. D.; Kapustin, V. N.; Eisele, F. L.; Weber, R. J.; McMurry, P. H. Geophys. Res. Lett. 1999, 26 (16), 2425−2428. (6) Brock, C. A.; Hamill, P.; Wilson, J. C.; Jonsson, H. H.; Chan, K. R. Science 1995, 270, 1650−1653. (7) Ball, S. M.; Hanson, D. R.; Eisele, F. I.; McMurry, P. H. J. Geophys. Res. 1999, 104, 23 709−23718. (8) Berndt, T.; Bőge, O.; Stratmann, F.; Heintzenberg, J.; Kulmala, M. Science 2005, 307, 698−700.

10133



Article

(9) Benson, D. R.; Young, L. H.; Kameel, F. R.; Lee, S.-H. Geophys. Res. Lett. 2008, 35, L11801. (10) Brus, D.; Hyvrinen, A.-P.; Viisanen, Y.; Kulmala, M.; Lihavainen, H. Atmos. Chem. Phys. 2010, 10, 2631−2641. (11) Howard, C. J. J. Phys. Chem. 1979, 83, 3−9. (12) Brown, R. L. J. Res. Natl. Bur. Stand. 1978, 83, 1−8. (13) Herrmann, E.; Brus, D.; Hyvarinen, A.-P.; Stratmann, F.; Wilck, M.; Lihavainen, H.; Kulmala, M. J. Phys. Chem. A 2010, 114, 8033− 8042. (14) FLUENT 6.2 User’s Guide; Fluent Inc.: Lebanon NY, 2005. (15) Hanson, D. R.; Eisele, F. L. J. Phys. Chem. A 2000, 104, 1715− 1719. (16) Hanson, D. R.; McMurry, P. H.; Jiang, J.; Tanner, D.; Huey, L. G. Environ. Sci. Technol. 2011, 45 (20), 8881−8888. (17) Weast, R. C., Ed. CRC Handbook of Chemistry and Physics; CRC Press Inc.: Boca Raton, FL, 1979. (18) Hanson, D. R.; Kosciuch, E. J. Phys. Chem. A 2003, 107, 2199− 2208. (19) Hanson, D. R.; Eisele, F. L. J. Geophys. Res. D 2002, 107, 4158. (20) Friedlander, S. K. J. Colloid Interface Sci. 1978, 67, 388. (21) Mirabel, P.; Clavelin, J. L. J. Aerosol Sci. 1978, 9, 219. (22) Hanson, D. R.; Lovejoy, E. R. J. Phys. Chem. A 2006, 110, 9525− 9528. (23) Green, D.; Perry, R., Eds. Perry’s Chemical Engineers’ Handbook; McGraw-Hill Professional Publishing: Blacklick, OH, 2007. (24) Clegg, S. L.; Brimblecombe, P.; Wexler, A. S. J. Phys. Chem. A 1998, 102, 2137−2154. Extended AIM Aerosol Thermodynamics Model. http://www.aim.env.uea.ac.uk/aim/aim.php. (25) Giauque, W. F.; Hornung, E. W.; Kunzler, J. E.; Rubin, T. T. J. Am. Chem. Soc. 1960, 82, 62−70. (26) Bein, K. J.; Wexler, A. S. J. Chem. Phys. 2007, 127, 124316. (27) McMurry, P. H. J. Colloid Interface Sci. 1980, 78, 513−527. (28) McGraw, R.; Wu, D. T. J. Chem. Phys. 2003, 118, 9337. (29) Vehkamäki, H.; Kulmala, M.; Napari, I.; Lehtinen, K. E. J.; Timmreck, C.; Noppel, M.; Laaksonen, A. J. Geophys. Res. 2002, 107 (D22), 4622. (30) Yu, F. J. Geophys. Res. D 2008, 113, D24201. (31) Herrmann, E.; Lihavainen, H.; Hyvärinen, A.-P.; Riipinen, I.; Wilck, M.; Stratmann, F.; Kulmala, M. J. Phys. Chem. A 2006, 110, 12448−12455. (32) Brus, D.; Neitola, K.; Hyvärinen, A.-P.; Petaja, T.; Vanhanen, J.; Sipila, M.; Paasonen, P.; Kulmala, M.; Lihavainen, H. Atmos. Chem. Phys. 2011, 11, 5277−5287. (33) Lovejoy, E. R.; Curtius, J.; Froyd, K. D. J. Geophys. Res. D 2004, 109, 08204. (34) Ayers, G. P.; Gillet, R. W.; Gras, J. L. Geophys. Res. Lett. 1980, 7, 433. (35) Roedel, W. J. Aerosol Sci. 1979, 10, 375−386.

10134


Computational Fluid Dynamics of a Cylindrical Nucleation Flow

Recommend Documents