A Monodisperse Population Balance Model for Nanoparticle Agglomeration in the Transition Regime

Nanoparticle agglomeration in the transition regime (e.g. at high pressures or low temperatures) is commonly simulated by population balance models for volume-equivalent spheres or agglomerates with a constant fractal-like structure. However, neglecting the fractal-like morphology of agglomerates or their evolving structure during coagulation results in an underestimation or overestimation of the mean mobility diameter, dm, by up to 93 or 49%, repectively. Here, a monodisperse population balance model (MPBM) is interfaced with robust relations derived by mesoscale discrete element modeling (DEM) that account for the realistic agglomerate structure and size distribution during coagulation in the transition regime. For example, the DEM-derived collision frequency, β, for polydisperse agglomerates is 82 ± 35% larger than that of monodisperse ones and in excellent agreement with measurements of flame-made TiO2 nanoparticles. Therefore, the number density, NAg, mean, dm, and volume-equivalent diameter, dv, estimated here by coupling the MPBM with this β and power laws for the evolving agglomerate morphology are on par with those obtained by DEM during the coagulation of monodisperse and polydisperse primary particles at pressures between 1 and 5 bar. Most importantly, the MPBM-derived NAg, dm, and dv are in excellent agreement with the data for soot coagulation during low temperature sampling. As a result, the computationally affordable MPBM derived here accounting for the realistic nanoparticle agglomerate structure can be readily interfaced with computational fluid dynamics in order to accurately simulate nanoparticle agglomeration at high pressures or low temperatures that are present in engines or during sampling and atmospheric aging.


Introduction
Nanoparticles made by gas-phase manufacturing or emitted by incomplete combustion of fossil fuels are abundant in our everyday lives [1]. Incipient nanoparticles grow by coagulation, sintering and surface growth. Inception [2], surface growth [3] and sintering [4] are active only in a narrow window of time, when the temperature is very high. Thus, coagulation is the dominant process in controlling nanoparticle morphology and number concentration, forming fractal-like agglomerates [5]. Such agglomerates are often present at very high concentrations during the gas-phase synthesis of carbon black, ceramic (TiO 2 and SiO 2 ) and metallic (Ni, Fe and Cu) powders, as well as soot emissions from engines, fires and volcanic plumes [1]. At high pressures or low temperatures that are present in engines [6] or during sampling [7] and atmospheric aging of aerosols [8], these agglomerates are formed by coagulation in the transition regime.
Agglomeration dynamics could be simulated by population balance models if an accurate collision frequency, β, and realistic particle morphologies are employed [5]. For example, sectional population balance models (SPBMs) are used to simulate agglomeration in the transition regime, with β obtained from the harmonic mean of those in the free molecular and continuum regimes [9] based on the agglomerate gyration, d g , mobility, d m , and primary particle, d p , diameters. The agglomerate structure, quantified by the relation between d g and d m [10], changes depending on the number [11], diameter and polydispersity of their constituent primary particles [12]. Accounting for the evolving agglomerate structure with SPBM is not trivial, as multiple equations per section need to be solved [13,14]. Therefore, models for the combustion synthesis of nanomaterials assume that d m = d g = d p n 0.56 p , based on a constant fractal dimension, D f = 1.8 [15,16] often with a constant value for d p [17,18]. Similarly, climate models simulate the coagulation dynamics of aerosols assuming that the spheres have the same volume-equivalent diameter, d v , with nanoparticle agglomerates [19].
With mesoscale discrete element modeling (DEM), the evolution of particle size distribution [3], morphology [10] and collision frequency [20] can be obtained from first principles. For example, it was recently shown that DEM-derived agglomerates reach a quasi-self-preserving size distribution (SPSD) with a mobility-based geometric standard deviation, σ g,m = 1.48 ± 0.03 [21], that is narrower than the SPSD attained in the free molecular regime [12] and in excellent agreement with data for SiO 2 [22] and soot [23,24] agglomerates. The DEM-derived collision frequency of polydisperse agglomerates at their quasi-SPSD was on average 82% higher than that of monodisperse ones, regardless of the chemical bonding and polydispersity of their constituent primary particles [21]. This DEM-derived collision frequency enhancement is on par with those measured from flamemade TiO 2 nanoparticle agglomerates in the transition regime [21]. Furthermore, power laws relating the mobility diameter normalized by the primary particle diamater, d m /d p , or gyration diameter, d m /d g , with the number of primary particles per agglomerate, n p , are derived by DEM simulations in the free molecular and transition regimes [3]. Such power laws can be used in population balance models to estimate accurately the evolving structure of nanoparticles during agglomeration in the transition regime.
Sectional models have been interfaced with computational fluid dynamics (CFD) in order to explain soot formation in diffusion flames [13,14] or predict agglomeration of soot nanoparticles from diesel engine exhausts [7,25]. Finite element methods have also been coupled with SPBMs to simulate coagulation of spheres [26], fractal-like agglomerates [27] and linear stacks (rouleaux) [28]. However, both SPBM and mesoscale simulations are computationally expensive and interfacing them with CFD is not trivial. Monodisperse population balance models (MPBMs) are computationally affordable and easy to use [16,29]. However, they apply best when particles have attained their SPSD and asymptotic fractallike structure by coagulation [30]. This is typically the case when high concentrations of nanoparticles coagulate and rapidly attain their quasi-SPSD [21]. Kruis et al. [16] used a three equation MPBM assuming monodisperse agglomerate and primary particle size distributions to track their coagulation and sintering and compared it to a two dimensional SPBM [31]. The influence of the fractal-like agglomerate structure on its β was accounted for by estimating the agglomerate collision diameter, d c , using D f , d p and n p . The MPBM predictions agreed well with those of the SPBM for the equivalent primary particle diameter. However, the agglomerate concentration was overpredicted, because the enhancement of β due to the polydispersity of the agglomerate size distribution was not considered. Goudeli et al. [11] interfaced a MPBM with an evolving DEM-derived D f to elucidate the agglomerate dynamics during coagulation and sintering in the free molecular regime. Neglecting the evolution of D f hardly affected d p , but overpredicted the agglomerate collision diameter up to 30% during the transition from hard-to soft-agglomeration (e.g., when the characteristic collision and coalescence times were comparable) [11]. In this regard, the MPBM for coagulation and surface growth [32] was coupled with DEM-derived power laws and coagulation rates in order to accurately describe the dynamics of soot nanoparticles in the free molecular regime. Extending this MPBM for coagulation in the transition regime is essential for the design of cleaner combustion engines, robust sampling devices for flame-made nanomaterials and accurate climate models for the estimation of the aerosol environmental impact.
The objective of this work is to develop an accurate but computationally affordable MPBM for nanoparticle agglomeration in the transition regime accounting for the evolving agglomerate structure and size distribution. Thus, the DEM-derived quasi-selfpreserving size distributions of agglomerates [21] and data-proven power laws governing their structure during their formation in the transition regime [10] are interfaced with a simple and accurate MPBM of their coagulation dynamics. In particular, the coagulation of nanoparticles is elucidated by a MPBM at low and high temperatures and pressures that are relevant during particle formation in engines [6] or sampling [7] and atmospheric aging of aerosols [8]. The collision frequency, number concentration, mobility and volumeequivalent diameters are validated with DEM simulations and data of soot nanoparticles at identical conditions.

Discrete Element Modeling and Agglomerate Structure
Ballistic and Brownian coagulation dynamics of nanoparticles in the absence of rotation, convection, deposition, van der Waals, electric or hydrodynamic forces are simulated by DEM of agglomeration in the free molecular and transition regimes [20]. This is a serial, event-driven algorithm implemented in C++, as detailed in [20]. In brief, spheres with mean diameter, d p , of 20 nm and volume fraction of 1 ppm are randomly distributed in a cubic cell at 1 bar and 1830 K. The bulk primary particle density is set to 1800 kg/m 3 , which is commonly used to simulate the dynamics of mature soot [3] and is close to that of SiO 2 [11]. Particles are in equilibrium with the surrounding gas and change direction after each collision or after traveling their persistence length [33]. The time between collisions is calculated with an event-driven method [34]. Particles stick to each other after collisions and the trajectory of the newly formed agglomerate is calculated by the momentum conservation principle. The surrounding gas pressure, P, is increased from 1 to 5 bar to simulate conditions close to those in combustion engines [6]. The geometric standard deviation, σ g,p , of the initial particle size distribution is also varied from 1 to 1.5 in order to be consistent with the measured soot [35] and metal oxide [22] primary particle size distributions.
The mobility diameter, d m , of DEM-derived agglomerates in the free molecular and transition regimes is related to their number of primary particles per agglomerate, n p , and d p by [10]: (1) and to their diameter of gyration, d g , by [10]: These power laws are obtained for soot particles formed in flames with maximum soot volume fractions spanning two orders of magnitude [10]. Furthermore, Equation (1) has been validated with data of soot, SiO 2 , ZrO 2 , Au and Ag aerosols [36].
The evolution of the detailed DEM-derived d m distribution from the free molecular to the transition regime has been validated with data of organic and inorganic nanoparticles from premixed [3], diffusion [21] and spray flames [21]. Figure 1 shows the DEM-derived geometric standard deviation, σ g,m , of the d m distribution as a function of the normalized mean, d m /d p , of agglomerates, with primary particles having σ g,p = 1 (a, b: broken lines), 1.2 (a: dotted line) and 1.5 (a: solid line) coagulating at 1830 K, and P = 1 (a, b: broken lines), 3 (b: dotted line) and 5 bar (b: solid line). Small σ g,p = 1 and 1.2 are representative for organic nanoparticles (e.g., soot [3,35]), while σ g,p = 1.5 is common for inorganic nanoparticles (e.g., SiO 2 [22] and ZrO 2 [37]) that attain their self-preserving size distribution by coalescence. Regardless of σ g,p , the σ g,m of the agglomerate size distribution increases up to 1.7 ± 0.05 in the free molecular regime at d m /d p = 3, consistent with DEM simulations of agglomeration and surface growth [10]. At d m /d p > 3, agglomerates coagulate in the transition regime and their σ g,m decreases, asymptotically attaining its quasi-selfpreserving σ g,m = 1.48 ± 0.03 [21] (thin solid line and shaded area). Increasing P from 1 to 5 bar accelerates the attainment of the σ g,m = 1.48 ± 0.03, as agglomerates coagulate mostly in the transition regime at these conditions. The large σ g,m attained by coagulation in the free molecular and transition regimes enhances the agglomerate coagulation frequency, as elaborated in the next section. Figure 1. Evolution of geometric standard deviation of agglomerate mobility diameter, σ g,m , as a function of the normalized mobility diameter, d m /d p , during coagulation (a) at P = 1 bar for agglomerates consisting of polydisperse primary particles with 1 ≤ σ g,p ≤ 1.5, and (b) for agglomerates of monodisperse primary particles coagulating at different pressures, 1 ≤ P ≤ 5 bar.

Agglomerate Dynamics With a Monodisperse Population Balance Model
Here, a monodisperse population balance model [16] is implemented in MATLAB (R2020a, MathWorks, Inc., Natick, MA, USA) and used to describe the coagulation dynamics of agglomerates in the transition regime. As the total mass and surface area of the agglomerates are conserved, one equation is enough to track the agglomerate number density and structure. At isothermal conditions, the rate of change in the total agglomerate number concentration, N Ag , is [5]: where β m is the collision frequency of monodisperse agglomerates given by the harmonic mean [9]: where β f m is the collision frequency in the free molecular regime [9]: where β co is the collision frequency in the continuum regime [9]: where k B is the Boltzmann constant, T is the temperature, d c is the agglomerate collision diameter that is equal to d p for spheres and d g for agglomerates (Equations (1) and (2)) to account for their fractal-like structure on β [38], m Ag is the agglomerate mass: where n p is given by: where N Ag,0 is the initial number density of the spherical primary particles. The collision frequency of polydisperse agglomerates, β p , is related to β m by [21]: The average enhancement factor of β p and its standard deviation have been derived based on 10 DEM simulations of soot and SiO 2 nanoparticle agglomeration at temperatures ranging from 1400 to 1830 K and pressures ranging from 1 to 10 bar [21]. Figure 2 shows the evolution of β as a function of d m measured (symbols) or estimated (lines and shade) for monodisperse (β m , Equation (4); broken line) or polydisperse (β p , Equation (9); solid line and shade) agglomerates of flame-made uncharged or weakly charged TiO 2 nanoparticles with d p = 20 nm at T = 295 K [39]. Agglomerates are initially formed by coagulation in the free molecular regime (Kn > 2.6). Therefore, their β rapidly increases with increasing d m . Agglomerates with d m larger than 100 nm coagulate in the transition regime (Kn < 2.6) and their β decreases towards the asymptotic β co in the continuum regime (Equation (6)). The measured β is underestimated up to a factor of 2 by β m that neglects the polydispersity of the agglomerate size distribution. The Fuchs interpolation for β m of monodisperse agglomerates in the transition regime results in a similar underestimation of the measured β [39]. In contrast, β p derived for polydisperse agglomerates [21] is 82 ± 35% larger than β m and in excellent agreement with the data. This validates the use of Equation (9) in the MPBM of nanoparticle coagulation in the transition regime.  Figure 3 shows the dynamics of the agglomerate structure quantified by the ratio of the mean mobility diameter over that of gyration, d m /d g , as a function of the mean number of primary particles per agglomerate, n p , obtained by DEM (dotted line and shade) and MPBM (solid line). Particles rapidly evolve from spheres with d m = 1.29 d g [40] into agglomerates with d m = d g having n p ∼ 10, and asymptotically reach d m = 0.7 d g [41] as n p increases. This rapid reduction of d m /d g is induced by the enhancement of the agglomerate inertia that determines d g [42]. Therefore, the agglomerate inertia becomes larger than its drag force at n p > 10 for monodisperse primary particles, resulting in d m /d g < 1. The DEM-derived evolution of d m /d g has been validated with data from wood combustion [10,43]. Most importantly, the agglomerate d m /d g estimated here by the MPBM using Equation (2) is in excellent agreement with that obtained by DEM from first principles. This confirms that the MPBM interfaced with DEM-derived power laws accounts for the agglomerate morphology dynamics in detail. Assuming a fixed value for d m /d g in the MPBM, such as 1.29 [18], 1 [15] or 0.7 [41], could significantly reduce its accuracy, as d m /d g = 1.29 is only valid for spheres and d m /d g = 0.7 is only reached for very large agglomerates [12]. In this regard, using d m /d g = 1 in the MPBM results in a 49% overestimation of the agglomerate d m , as discussed in Section 3.4. Figure 3. Evolution of the ratio of the mean mobility diameter over that of gyration, d m /d g , as a function of n p during agglomeration derived by DEM (dotted line and green shade) and the MPBM interfaced with Equation (2) (solid line). The d m /d g is initially equal to 1.29 for single spheres [40], but rapidly decreases to 1 for agglomerates with n p ∼ 10 and σ g,p = 1 and asymptotically reaches 0.7 (dashed line) for large agglomerates with n p ≥ 700 [41]. The DEM-derived relation [10] (solid line, Equation (2)) quantifies the evolution of d m /d g within 6% of the DEM simulations. Figure 4 shows the evolution of (a) collision frequency, β, (b) number density, N Ag , (c) mean mobility, d m , and (d) volume-equivalent, d v , diameters as a function of time, t, for agglomerates consisting of primary particles with mean d p = 20 nm and geometric standard deviation, σ g,p = 1 (broken lines), 1.2 (dotted lines) and 1.5 (solid lines) at T = 1830 K and P = 1 bar derived by DEM (thick lines) and the MPBM (thin lines and shades). The shades quantify the statistical variation of β p (Equation (9)). Increasing σ g,p enhances β, resulting in a faster reduction of N Ag with time as shown in Figure 4a. However, increasing σ g,p also delays the attainment of the tranisition regime, as agglomerates consisting of primary particles with larger σ g,p have smaller d m and d v . This is due to the small primary particles in those agglomerates that have marginal impact on increasing d m and d v . The MPBM-derived N Ag , β, d m and d v are in excellent agreement with those obtained by DEM for all σ g,p investigated here. Therefore, the present MPBM could be used to simulate the coagulation dynamics of organic nanoparticles with rather monodisperse primary particle size distributions [3,35], and also those of metals [44] and metal oxides [22,37] that attain the self-preserving size distribution by coalescence in practical applications. The MPBM-derived agglomeration dynamics are also in excellent agreement with those obtained by DEM for nanoparticles with d p = 10, 20 and 40 nm (Supplementary Materials: Figure S1). The present MPBM has been also validated with DEM simulations and measurements of non-spherical aggregated soot nanoparticles coagulating in the free molecular regime [32]. Figure 5 shows the evolution of (a) β, (b) N Ag , (c) d m and (d) d v as a function of t at P = 1 (broken lines), 3 (dotted lines) and 5 bars (solid lines) during the coagulation of monodisperse primary particles with d p = 20 nm derived by DEM (thick lines) and the MPBM (thin lines and shades). Agglomerates at 1 bar are intially in the free molecular regime. Thus, their β increases as their n p and d g increase rapidly by coagulation. As the agglomerates grow, they enter the transition regime where their β gradually decreases towards an asymptotic value that is rather independent of agglomerate size in the continuum regime [45]. Increasing pressure from 1 to 5 bar decreases β by up to a factor of 2.6, as agglomerate diffusivity decreases at higher pressures [46]. This results in larger N Ag at 3 and 5 bars compared to those obtained at 1 bar. Agglomerate d m and d v derived at 1 bar are at least 50% larger than those obtained at 3 and 5 bar. This is due to their large collistion frequency at 1 bar (Figure 5a) that takes place initially in the free molecular regime. Agglomerate d v is on average 35% smaller than its d m at all pressures, quantifying the ramified non-spherical agglomerate structure [42]. The MPBM-derived N Ag , β, d m and d v are in excellent agreement with those obtained by DEM for all P relevant for soot formation in engines investigated here [6]. The MPBM-derived agglomeration dynamics are also in excellent agreement with those obtained by DEM for initial N Ag,0 and (incubation) residence times spanning 6 and 8 orders of magnitude, respectively ( Figure S2).  Figure 6 illustrates the agglomerate effective density, ρ e f f , as a function of the normalized mobility diameter, d m /d p , measured for soot particles sampled from premixed ethylene flames with an equivalence ratio (EQR) of 2 (triangles) or 2.4 (squares), and estimated here accounting for the evolving fractal-like agglomerate structure (Equations (1) and (2) [10], solid line) or assuming a constant agglomerate structure with d m = d g = d p n 0.56 p (broken line) [16]. Ramified agglomerates are formed at EQR = 2 and 2.4 having average d p = 9 [35] and 19.6 nm [47], respectively. The agglomerate d m /d p increases during coagulation, reducing ρ e f f and d m /d g from 1.29 to 0.7 (as shown in Figure 3). Neglecting the evolving agglomerate structure and assuming d m = d g = d p n 0.56 p underpredicts the measured ρ e f f by up to 50%. In contrast, the agglomerate ρ e f f obtained here using DEM-derived power laws (Equations (1) and (2)) that account for the realistic agglomerate morphology is in excellent agreement with data of soot nanoparticles from different combustion conditions. The impact of the agglomerate morphology on the coagulation dynamics of nanoparticles in the transition regime is investigated next. In particular, Figure 7 shows the agglomerate N Ag (a, b), mean d m (c, d) and d v (e, f) of soot nanoparticles measured as a function of t during coagulation at T = 295 K (symbols). Soot nanoparticles are sampled from the premixed flames with EQR = 2 (a, c, e) or 2.4 (b, d, f) shown in Figure 7. The measured N Ag , d m and d v (symbols [7]) are compared to those estimated by the MPBM (lines, shades) assuming volume-equivalent spheres (dotted lines) or agglomerates with constant (broken lines) or evolving morphology (solid lines).

Validation of Low-Temperature Coagulation Dynamics with Experiments
Approximating agglomerates with volume-equivalent spheres underpredicts their d m and d v by up to 93 and 18%, respectively, and overpredicts their N Ag by up to 55%. The N Ag and d v obtained by the MPBM for agglomerates with a constant or evolving structure are in excellent agreement with the data of soot obtained from both EQRs. However, neglecting the evolving agglomerate structure during coagulation in the transition regime overpredicts d m by up to 49%. Thus, interfacing the MPBM with DEM-derived power laws (Equations (1) and (2)) is essential to account for the realistic agglomerate structure and accurately estimate the d m dynamics so that they are in excellent agreement with data.

Conclusions
A simple monodisperse population balance model (MPBM) is derived here for the agglomeration of nanoparticles in the transition regime. The MPBM uses relations derived from detailed discrete element modeling (DEM) simulations in order to obtain the evolving structure of agglomerates quantified by their mobility and gyration diameters. As a result, the DEM-derived collision frequency, β, that accounts for the agglomerate size polydispersity is 82 ± 35% larger than that of monodisperse agglomerates and in excellent agreement with measurements of flame-made TiO2 nanoparticles. Therefore, the N Ag , d m and d v derived by the MPBM, accounting for the evolving fractal-like structure of agglomerates (Equations (1) and (2)) and the impact of their polydispersity on β (Equation (9)), are on par with those obtained by detailed DEM simulations for both monodisperse and polydisperse primary particles coagulating at pressures P = 1-5 bar. Most importantly, the soot agglomerate N Ag , d m and d v derived during coagulation at low temperatures are in excellent agreement with data from various premixed flame conditions. In contrast, neglecting the fractal-like morphology of agglomerates or their evolving structure during coagulation results in an underestimation or overestimation of the mean d m by up to 93% or 49%, repectively. Thus, the MPBM derived here accounting for the realistic nanoparticle agglomerate structure can be readily interfaced with CFD in order to accurately simulate the agglomeration dynamics of nanoparticles at high pressures or low temperatures that are present in engines or during sampling and atmospheric aging.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/ma14143882/s1, Figure S1: Evolution of (a) collision frequency, β, (b) number density, N Ag , (c) mean mobility, d m , and (d) volume-equivalent, d v , diameters as a function of time, t, for agglomerates consisting of monodisperse primary particles with d p = 10 (solid lines), 20 (broken lines) and 40 nm (dotted lines) derived by DEM (thick lines) and MPBM (thin lines and shaded areas) at T = 1830 K and P = 1 bar. The MPBM-derived agglomerate dynamics are in excellent agreement with those obtained by DEM for the wide range of d p studied here. Figure S2: Evolution of (a) β, (b) N Ag , (c) d m , and (d) d v as a function of t for monodisperse primary particles with initial N Ag,0 = 2 · 10 15 (solid lines), 2 · 10 17 (broken lines) or 2 · 10 19 m −3 (dotted lines) and dp = 20 nm derived by DEM (thick lines) and MPBM (thin lines and shaded areas) at T = 1830 K and P = 1 bar. The MPBM-derived agglomeration dynamics are in excellent agreement with those obtained by DEM for N Ag,0 and (incubation) residence times spanning 6 and 8 orders of magnitude, respectively.