1 Introduction

Aerosols, fine solid particles or liquid droplets suspended in gas, are ubiquitous both in industrial products and in nature, for example, soot in flames, cement dust, smog, smoke, and clouds. Aerosol dynamics generally involve three fundamental processes: nucleation, coagulation, and condensation. The evolution of aerosol particles is highly sensitive to the history of vapor concentration and temperature that these particles undergo in turbulence. Therefore, it is very difficult to study the intricate interaction between turbulent mixing and aerosol dynamics.

In their seminal work, Lesniewski and Friedlander^[1] measured the nucleation and growth of dibutyl phthalate (DBP) particle in a free turbulent jet. Since then, a few numerical simulations have been done to investigate aerosol particles formation and growth in turbulent jets. Garmory and Mastorakos^[2] and Veroli and Rigopoulos^[3] simulated DBP aerosol evolution in the turbulent jet flow the same as the experimental configuration of Ref. ^[1]. In the simulation of Garmory and Mastorakos^[2], Reynold's stress model for turbulence was used, combined with the stochastic fields transported probability density function (PDF) method to model aerosol evolution. Large eddy simulation was used by Zhou and Chan^[4] in the simulation of a similar configuration as that of Ref. ^[1]. Veroli and Rigopoulos^[3] developed the transported population balance equation-PDF approach, based on a $k$-$\epsilon$ model for turbulence. Garmory and Mastorakos^[2] and Veroli and Rigopoulos^[3] found discrepancies when comparing their simulation results with the measurements of, and attributed those discrepancies to the inability of the numerical methods to correctly model the coupling between turbulence and aerosol dynamics, and the uncertainty in the aerosol kinetics and the experimental procedure.

In a previous study, direct numerical simulation (DNS) combined with the quadrature method of moments (QMOM)^[5-6] was used to simulate DBP aerosol evolution in a turbulent mixing layer to investigate the intricate interaction between aerosol dynamics and turbulence. A Lagrangian particle scheme^[7] was used to deal with the moments set, which was used to describe the aerosol particles. Due to the limitation of the method of moments, the aerosol particle size distribution (PSD), which carries the most complete and useful information of the aerosol particle, was not available in the previous simulation. In this work, Monte Carlo (MC) simulation of aerosol evolution is performed along the Lagrangian trajectories obtained in the previous simulation so as to investigate the aerosol PSD. This research provides a successful prototype to simulate aerosol evolution in turbulent flows through combing MC simulation and DNS in a feasible and efficient way.

2 Methodology 2.1 Fluid dynamics

The DNS for the flow field was performed with the parallel flow solver "NGA"^[8], which uses a finite difference method on a spatially and temporally staggered grid. The semi-implicit fractional step method^[9] is used for temporal integration. The second-order centered scheme is used to discretize the velocity. The message passing interface is used for parallel simulation. The library HYPRE^[10] is used to solve the pressure Poisson equation, which is highly efficient on massively parallel machines.

The velocity profile at the inlet is to join two laminar boundary layers. At the outlet, the free convective outflow^[11] condition is used. The periodic boundary condition is used in the spanwise direction $z$. The free slip boundary condition is used in the crosswise direction $y$.

The two streams have velocities $u_1=15$~m/s and $u_2=5$~m/s. The random velocity fluctuation is $0.4(\omega_t-1) u$, superimposed on the velocity profile, where $\omega_t$ is a random variable from the standard uniform distribution. The computational domain is discretized with $768\times398\times256$ (80 million grid points). A posteriori analysis shows that the size of the grid is enough to resolve the Kolmogorov scale.

The temperature $T$ and the vapor mass fraction $Y$ are assumed to be passive scalars and are modeled by the convection-diffusion equation. The effect of aerosol dynamics on the temperature field is neglected due to low vapor mass fraction. However, the vapor concentration consumption is considered, i.e.,

where ${u}$ is the gas flow velocity, $D$ is the diffusion coefficient, and $S_{\mathrm{Aerosol}}$ is the source term, which comes from the gas-to-particle conversion. Both the Lewis number and the Schmidt number are assumed to be one, i.e., thermal diffusivity, mass diffusivity, and kinematic viscosity are all the same. At the inlet, a hyperbolic tangent profile is used for the temperature with $T_{\mathrm{max}}=400$ K and $T_{\mathrm{min}}=300$ K at the fast and slow streams, respectively. Similarly, the profile for the vapor mass concentration is imposed with $Y_{\mathrm{max}}=5\times 10^{-3}$ and $Y_{\mathrm{min}}=0$.

More than a million CPU hours was used for the simulation on the supercomputer, "Shaheen", at King Abdullah University of Science and Technology. A schematic view of the simulation configuration is shown in Fig. 1.

2.2 Aerosol dynamics

In this work, all aerosol particles (droplets) are taken to be spherical, characterized by the diameter $\xi$. The PSD $n({x}, t;\xi)$ satisfies the general dynamic equation (GDE)^[12]

where the three terms on the right-hand side denote nucleation, condensation, and coagulation, respectively. In the GDE, the diffusion term is neglected because of the large Schmidt number of aerosol particles. For example, under standard atmospheric conditions, spherical particles with a diameter of 10 nm have a Schmidt number equal to 290, and even higher Schmidt numbers for larger particles^[12]. Nucleation is modeled by the classic Becker-Döring theory. Condensation and coagulation are usually modeled differently in different Knudsen number regimes, which are defined according to the particle size compared to the gas mean free path, i.e., free molecular and continuum regimes. The harmonic mean of the models in the two regimes are used. Details of the physical models can be found in Friedlander^[12] and Zhou et al.^[13].

Instead of solving the PSD directly from Eq. (2) , only a few moments (four moments in this simulation) are solved in the QMOM, which is much more economic computationally than the former and is able to capture the main characteristics of the aerosol statistics. The moment (of $k$th ($k=0, 1, 2, \cdots$) order) is defined by

Therefore, $M_0$ is the particle number density (which is also denoted as $N_\mathrm d$ thereafter, with unit $1/\rm{m}^3$), $M_1$ is the diameter "density" ($\rm{m}/\rm{m}^3$), $\pi M_2$ is the surface area "density" ($\rm{m}^2/\rm{m}^3$), and $\frac{\pi}{6} M_3$ is the volume fraction (which is also denoted as $F_\mathrm v$ thereafter, with unit $\rm{m}^3/\rm{m}^3$), etc. The details of the QMOM can be found in Ref. ^[5]. Actually, there are various methods of moments, such as TEMOM^[14-16], DQMOM^[17], and MOMIC^[18]. It is found that the QMOM is the most suitable for this simulation after taking efficiency, precision, and robustness into account.

2.3 Lagrangian particle scheme

The Lagrangian particle scheme circumvents the moments realizability problems^[19] encountered in conventional advection schemes for aerosol dynamics. Besides, the low numerical diffusion property of the Lagrangian particles scheme makes it easy to track the sharp front of the motion of aerosol particles, and makes it an ideal tool to investigate the differential diffusion of aerosol particles and the gas.

In the Lagrangian particles scheme (here, particles refer to mathematical objects, other than physical particles), the movement of a large number of Lagrangian particles is tracked. Physical variables of interest (aerosol moments here) are tied to these Lagrangian particles and evolve along the Lagrangian trajectories. Under this frame work, the conservation equations (here, aerosol diffusion is neglected) are separated to two parts, convection and unsteady evolution. The convection is dealt with through an average scheme over particles, and the unsteady evolution of the variables along a trajectory is obtained by solving the corresponding control equations, which are only ordinary differential equations. The details of the scheme can be found in Ref. ^[7].

2.4 MC simulation

The aerosol evolution along a fluid parcel trajectory is simulated by an MC method developed by Zhou et al.^[20]. A number of notional simulation particles are used to simulate the aerosol particles evolution. The notional particles undergo various aerosol processes to model the aerosol dynamics. Notional particles are generated in nucleation process, and these particles grow via condensation process. Coagulation is simulated through the interaction of the notional particles by the stochastic algorithm^[21]. In principle, the coagulation algorithm can be described as:

Step 0 Give initial conditions (diameter, surface, volume, etc.) for $N$ particles, decide the simulator size $V=N/n_0$ from the initial number density $n_0$, and calculate the coagulation rate $C_{ij}=\beta(v_i, v_j)/V (i=1, 2, \cdots, N-1$ and $j=i+1, i+2, \cdots, N)$, where $\beta$ is the coagulation coefficient to describe the rate that two particles coagulate.

Step 1 Produce a random time $\tau$, where $\tau$ follows the Poisson distribution

Here, $C_0$ denotes to sum all coagulation rate, i.e., $C_0=\sum\limits^{N-1}_{i=1}\sum\limits^{i+1}_{j=1}C_{ij}$.

Step 2 Randomly pick up two particles to coagulate by the chance $P(i, j)=C_{ij}/C_0$.

Step 3 Update $N$ and $C_{ij}$, and then repeat Steps 1 and 2 till the end of the simulation time $t_{\rm stop}$.

More details of the MC algorithm can be found in Ref. ^[20]. It is worth pointing out that MC methods for the aerosol dynamics simulation are very popular, such as Refs. ^[22] and ^[24], just to name a few.

The MC simulation takes temperature and vapor concentration along Lagrangian trajectories as input, and it simulates the PSD evolution along the trajectories. MC is computationally far more expensive then the QMOM, that is why MC is not coupled with the fluid dynamics solver, but it is used during post-processing. However, MC is generally more flexible and accurate. Moreover, MC can provide the full PSD, while the QMOM can only predict a few lower order moments (number density, volume fraction, etc.) of the PSD.

3 Results and discussion 3.1 Aerosol evolution overview

The main findings in the previous relevant simulation^[13] are shortly reviewed to provide a convenient overview of the aerosol evolution in the turbulent mixing layer. For the number density, high value is found in the laminar region near the inlet, which is due to high nucleation rate and weak mixing there. In the turbulent region, high number density is typically found in the outer region of a vortex. Very differently from the number density, the volume fraction increases along the streamwise direction, because of the increase of the aerosol residence time. High volume fraction typically appears in the inner region of a vortex, because condensation growth has the highest rate at the intermediate temperature. Typically, particles on the hot rich vapor side are larger than those on the cold lean vapor side. A persistent narrow layer of nanometer-sized particles is found on the cold lean vapor side.

3.2 Aerosol evolution along sample trajectory

In the Lagrangian particle scheme, around 500 million particles are tracked through the flow domain. Figure 2 shows the temperature profile along a representative Lagrangian particle trajectory. The corresponding nucleation rate is also shown. Since vapor depletion is negligible^[13], the vapor concentration is almost linearly correlated with the temperature. Hence, the nucleation rate is almost solely determined by the temperature. Nucleation mostly takes place in a narrow temperature region (marked with the gray block in Fig. 2) . Large number of aerosol particles are generated when a fluid parcel passes this active nucleation region (in the temperature space).

The MC simulation along Lagrangian trajectories is used to verify the results obtained from the QMOM. Figure 3 compares the number densities and volume fractions obtained by the QMOM and the MC along the selected trajectory. The QMOM and the MC give the same result for the number density. The number density profile ($N_\mathrm d$) behaves like a staircase, and is almost constant between abrupt jumps, which corresponds to the bursts of aerosol particles from nucleation. The volume fraction given by the QMOM also agrees with the MC result very well. The volume fraction increases all the way due to condensation. Although the QMOM and the MC agree with each other extremely well here, this agreement cannot be guaranteed in general. In this simulation, the number density is dominated by the nucleation process, and coagulation is negligible owing to the low number density. Nucleation is only dependent on the temperature and vapor concentration and is independent of the existing aerosol particles. Therefore, the same number density profile is observed. Moreover, condensation growth is mainly in the free molecular regime, which is independent of the particle size. In other words, the evolution of the volume fraction is insensitive to the actual PSD. When coagulation is significant and condensation has stronger dependence on particles size, the results from the QMOM and the MC might differ by a greater extent.

One big advantage of the MC over the QMOM is that the MC can provide the PSD. Figure~4 shows the PSD at the end of the trajectory (where $d_\mathrm p$ is the particle diameter). The results from 300 and 500 thousand samples are compared to check the stochastic convergence. The PSD shows complex multi-modal distribution. The peaks in the PSD are closely related to the nucleation rate curve in Fig. 2. The rightmost peak in the PSD corresponds to the leftmost spiky peak in the nucleation rate curve,where the nucleated particles have the longest residence time to grow into the largest ones. Particles in the big peak in the PSD from 100 nm to 400 nm are generated in the nucleation burst during 0.01 s to 0.018 s in Fig. 2. The big trough in the PSD corresponds to the dip in the nucleation rate around 0.02 s. Particles under a few tens of nanometer in the PSD are generated in a narrow time period after 0.02 s. The relatively broad distribution of these nanometer-sized particles is due to the fast condensation growth rate for small particles. As this example PSD shows,the combination of nucleation and condensation can render very complicate distribution. The distribution may even be discontinuous when nucleation ceases to generate particles for an extensive period of time and the particles still grow due to condensation.

3.3 Mean particle size distribution

The instantaneous PSD is very difficult to obtain experimentally. The mean PSD is usually reported by continuous sampling at a fixed point in experiments. Here, the aerosol PSD is simulated along Lagrangian trajectories. To calculate the mean PSD at a spatial point, Lagrangian trajectories passing through the point (actually a small region) at the same time instant are collected. The PSDs along these trajectories are averaged with weighing on the corresponding number density to obtain the mean PSD, i.e.,

where $m$ is the number of trajectories, and $\Psi_{{\rm PSD}i}$ and ${N_\mathrm d}_i$ are the PSD and the number density along trajectory $i$, respectively.

Lagrangian trajectories are collected at 12 points, summarized in Table 1. These points are located on 4 cross planes, i.e., $x^*=125$, 249, 415, and 581. On each plane, 3 points are selected, with $y^*$ corresponding to $\overline{\phi}=0.3$, 0.5, and 0.7, where $\phi$ is the mixture fraction, which is the normalized temperature to the range from 0 to 1.

Figure 5 shows trajectories passing through the point at $x^*=581$ with $\overline{\phi}=0.5$. For clarity, only around one tenth of the total trajectories are shown in Fig. 5. It is clear that the point (corresponding to $\overline{\phi}=0.5$) is below the line of $y^*=0$, which means that the trajectories tend to "bend down" from the fast stream to the slow stream. Trajectories originated around $y^*=0$ fluctuate more strongly, which reflects that mixing is the strongest in the middle of the mixing layer. Trajectories affected by the large roller structure can also be easily identified, such as the upper most one.

Figure 6 shows the mean PSDs (normalized by the total number density) at the sampling points summarized in Table 1, which are calculated according to Eq. (5) . In addition to collecting all trajectories passing through the same spatial point (illustrated as in Fig. 5) , trajectories at 25 different time instants are collected as well. The time span is big enough to avoid significant statistical bias when evaluating the mean PSD. Total around 10 thousand trajectories are used to compute the mean PSD at one point. Overall, the PSDs are similar with each other, with a high peak at small diameter position, which corresponds to the huge number of nucleated particles, and with a sharp decrease at large diameter position, because the condensation growth rate decreases sharply when particle size becomes bigger, and hence the biggest particles are of similar size. The peak on the right-hand side is also related to the typical function of the condensation growth rate, which makes particles to accumulate to the right on the diameter space. Along the cross-wise direction (for $\overline{\phi}$ from 0.3 to 0.7) , the PSD at higher mixture fraction ($\phi$) has a lower dip around $d_\mathrm p=3$ nm and usually a higher peak on the right-hand side, because small particles on the hot vapor side (high $\phi$) tend to grow faster. Along the streamwise direction (for $x^*$ from 125 to 581) , the PSD evolves from bi-modal to tri-modal, and then to bi-modal finally. The tri-modal distribution is believed to relate to the big roller structure of the mixing layer (see Fig. 3 in Ref. ^[13]) in the region transition to turbulence. In the fully turbulent region, the distribution is bi-modal. In the measurement of aerosol PSD in a turbulent jet^[1], both bi-modal and tri-modal distributions were observed. Lesniewski and Friedlander^[1] considered the formation of tri-modal distribution is due to coagulation. However, the number density is too low (they were aware of this point) in their experiment for coagulation to be significant. Therefore, it is anticipated that the tri-modal distribution observed in the experiment^[1] may be also due to the large flow structure.

4 Conclusions

In a former work, the evolution of DBP particles in a turbulent mixing layer was investigated by means of DNS for the flow field, combined with the QMOM for aerosol dynamics. In this work, further MC simulation of the evolution of aerosol particle along the Lagrangian trajectories obtained in the former simulation is performed.

Along a fluid parcel moving through the turbulent flow, temperature and vapor concentration exhibit complex fluctuations, which trigger complicate aerosol processes. The particle size distribution exhibits complex modality due to the synergistic effect of nucleation and condensation (coagulation is not significant here). However, the mean PSD is found to be bi-modal in most of the mixing layer, except that a tri-modal distribution is found in the turbulent transition region due to the large roller structure. In the measurement of aerosol PSD in a turbulent jet^[1], both bi-modal and tri-modal distributions are observed. Lesniewski and Friedlander^[1] considered the formation of tri-modal distribution is due to coagulation. However, the number density is too low (they were aware of this point) in their experiment for coagulation to be significant. Therefore, it is anticipated that the tri-modal distribution observed in the experiment^[1] may be also due to the large flow structure.

Acknowledgements The DNS simulation was done with the valuable help of Anotonio Attili and Fabrizio Bisetti at King Abduallah University of Science and Technology.