Modeling Nanoparticle Agglomeration using Local Interactions

Abstract

Nanoparticle agglomeration plays an important role in processes such as spray drying and particle flame synthesis. These processes have in common that nanoparticles collide at low concentrations and get irreversibly linked at the point of contact due to plastic deformation. In this article, we investigate several models of irreversible connections, which require only local interactions between the colliding nanoparticles and thus allow for scalable simulations. The models investigated here connect the particles upon collision by non-bonded strongly attractive interactions, bonded interactions or by binding agents placed at the point of contact. Models using spherically symmetric interactions form compact agglomerates and are therefore unsuitable to study agglomeration. In contrast, models that are either based on both central and angular potentials (type one) or on binding agents (type two) efficiently prevent restructuring of the agglomerates, and are therefore useful for modeling contacts formed by plastic deformation. Moreover, both types of models allow to control the rigidity and by that the degree of restructuring. The first type of model is computationally more efficient at low fractional dimensions of the aggregates, while the second gives easy access to local shear forces, which is important when breaking of agglomerates is to be considered. As example applications, we reproduce the well-known diffusion-limited agglomeration (DLA) and report results on soot aggregation.

Copyright 2014 American Association for Aerosol Research

1. INTRODUCTION

Particle agglomeration is a process in which small particles with a diameter of 1–100 nm collide and stick to form larger agglomerates. Agglomeration is the main growth mechanism in spray drying and particle flame synthesis, but also plays an important role in water purification, mineral beneficiation and biological separation processes (Kodas and Hampden-Smith 1999; Gregory 2005). Experimental studies provide information on the size and shape of selected agglomerates, but not on the history of the agglomeration process or the intermolecular forces within the agglomerates that are of importance for the dynamics of the formation of the final product (Forrest and Witten 1979; Witten and Sander 1981; Mountain and Mulholland 1988; Lall and Friedlander 2006). Computer simulations can help to gain more insight, since they allow for tracking of individual particles and by that provide a better understanding of the agglomerate’s growth and morphology.

The relevance of simulations crucially depends on the underlying model, since the agglomerating particles can be of very different nature, from crystallites to small protein aggregates. The particles feel short-range attraction due to induced dipole forces that results in an elastic deformation and thus strengthens the contact, or they create bridges from one particle to another (Kinloch 1987; Gay 2002; Li et al. 2011). Usually, this attraction is so strong that the contacts can be considered permanent. In addition, the particles are dispersed in a fluid that mediates flow and Brownian motion (Zaichik et al. 2008; Das and Garrick 2002; Babler et al. 2005), so that particles adhere only when the attractive force overcomes the thermal fluctuations and the hydrodynamic drag (Sander et al. 2009; Vanni and Baldi 2002). The properties of the formed agglomerates therefore depend both on the modeling of the attractive particle interactions and the fluid.

Peng et al. (2010) used a discrete element method (DEM) to simulate the nanoparticle agglomeration due to random Brownian diffusion. Even though the simulation results were found to be in good agreement with experimental data, the coordination numbers were too high. Binder et al. (2008) showed that hydrodynamic interactions lead to compact agglomerates that are exposed to high drag forces. Recently, Isella and Drossinos (2010) investigated numerically the agglomeration of spherically symmetric nanoparticles that undergo Brownian motion and stick irreversibly when coming into contact. The aggregates were found to be compact, tubular and elongated.

These numerical studies highlight a common issue, namely that the agglomerates obtained from simulations are often too compact. The reason for this is sliding and/or the rolling of particles within an agglomerate. This rearrangement, which leads to rather compact agglomerates as time proceeds, may or may not be physical. Particles that are tied together by weak physical van der Waals forces (soft agglomerates) may restructure when heated or exposed to high tension (Schmidt-Ott 1988). However, the particle bonds within hard agglomerates are sufficiently strong and restructuring after collision does not appear, so that compact clusters are not found. This may also be true for soft agglomerates such as agglomerates of real silica or titania that tend to break rather than to restructure (Dominik and Tielens 1997).

Thus, in order to model strongly binding particles in computer simulation, contact bonds that prevent rearrangement are necessary. Kusaka et al. (2011), Kempf et al. (1999), and Iglberger et al. (2008) coupled the nanoparticle motion models with a rigid body motion of the formed aggregates. Although rigid body dynamics solves the undesired restructuring, it is a nonlocal interaction that requires to exchange data between all processors that a particular agglomerate occupies, and thus limits the scalability to large systems.

In our study, we investigate several models of particle binding that rely only on local interactions, and by that naturally allow for scalable simulations. The contact binding is achieved by simple van der Waals attraction, formation of harmonic and/or angular bonds, or by binding agents at the contact points (cp), so that the particles loose their rotational symmetry (Table 1). The Brownian motion of the particles is modeled by Langevin dynamics (LD) (Limbach et al. 2006), and for simplicity, the particles only feel a weak van der Waals attraction, apart from the contact binding interactions. We show that central potentials always lead to compact agglomerates, independently of whether they are introduced as an always present nonbonded interaction with finite well depth or as harmonic bonds. In contrast, models that either introduce angular bonds or binding agents and thus break rotational symmetry can very effectively prevent restructuring at low computational cost. Moreover, the models allow to control the rigidity of the agglomerates and thus control restructuring accurately.

TABLE 1. Summary of the four different agglomeration formation models examined in this report

Model number	Model name	Model description	Acronym
i	Nonbonded potential model	Collision model only with nonbonded potential	NB
ii	Single bond model	Single bond potential between collided particles model	SB
iii	All bonds model	Single bond and angular bond potentials model	AB
iv	Angle bonded virtual sites model	Collision model with two virtual particles connected through angle potential	AnBV

This article is organized as follows: Section 2 introduces the potentials between the particles and the binding models. Section 3 reports simulation results on the fractal and morphological properties of the agglomerates for each model. Section 3 then investigates the noncentral interaction models in more detail, and compares present results to known results for the classical diffusion-limited aggregation (DLA) and to results of other numerical and experimental studies on the formation and growth of soot aggregates.

2. NUMERICAL METHODS

Simulations are performed using the Molecular Dynamics (MD) software package ESPResSo (Arnold et al. 2013). It provides all the required interactions plus a number of additional features such as electrostatic or hydrodynamic interactions, that might be of interest for future studies of agglomeration.

In this section, we first describe the interactions between the particles, then explain our models of binding of the particles that have collided and finally give details on our test simulations such as the equations of motion employed.

2.1. Relevant Interactions Between Particles

The potential between two primary particles is modeled by integrating the intermolecular potential over the volume of these particles. The traditional (12-6) Lennard-Jones potential is used to model the intermolecular potential. It combines a long-range attractive force (the 1/r⁶-term) with a short-range repulsive force (the 1/r¹²-term) between the particles, [1] where r is the particles’ distance, σ_LJ is the distance at which the potential equals zero and ε is the depth of the attractive potential, defining the maximum attractive energy between two particles.

The potential between two particles (U(r)) is obtained from the integration of Equation (1[1] ) over two spherical particles of diameter σ (Lazaridis and Drossinos 1998; Narsimhan and Ruckenstein 2008). It is given by [2a] [2b] [2c]

In Equation (2[2a] ), A is the Hamaker constant and expressed by [3] where n is the molecular number density in the solid (Isella and Drossinos 2010).

In the present study, the interaction between two clusters is obtained via the particle potentials between the respective particles that belong to the different clusters and are closer in distance than a cut-off distance r_cut. In addition, adjacent particles can be bound via a 2-body spring potential and 3-body angular potentials. The 2-body spring bond (harmonic bond) potential between two particles in contact is given by [4] where gives the distance between the particles, r₀ is the equilibrium distance, and k is the spring constant. The 3-body angular bond is used to fix the angle between the position vectors from the center particle to the two other particles. The energy configuration of the angular bond potential is [5] where θ is the angle in radians between vectors and , θ₀ is the equilibrium angle, and k is the angle bond constant.

2.2. Equations of Motion

LD is used to model the particle motion including Brownian motion. We use a Velocity-Verlet-integrator plus a Langevin thermostat to integrate the Langevin equation (Limbach et al. 2006). At each time step all particles are subjected to a random force and a frictional force such that these two forces satisfy the fluctuation-dissipation theorem and balance each other. In this formalism, the Langevin equation for the i-th particle (Kubo 1966) is given by [6] where F_Ci are conservative forces that arise from inter-particle potentials that the ith particle feels (F_Ci= –▽ U_i), γ is the friction constant between the primary particles and surrounding fluid, u_i is the velocity of the ith particle and W_i(t) is Einstein’s white noise term, which is a Gaussian random source (Lemons 1908; Fox 1978). The noise models the random kicks of the fluid molecules to the particles with zero-mean and satisfies [7] The exact modeling of the friction coefficient may also depend on the Knudsen number Kn, which is the ratio of the fluid molecules mean free path length l_fluid to the particle radius, Kn = 2l_fluid/σ (Knudsen and Weber 1911). For large particles (Kn ≪ 1), the flow is in the continuum regime and the well known Stokes law applies, viz. [8] with μ_f denoting the dynamic viscosity of the surrounding fluid. For smaller particles (Kn > 1), however, a slip velocity between the particle and the surrounding fluid will exist. The friction coefficient will be lower than the Stokes law would predict and the so-called Cunningham correction, C_C, needs to be introduced as is indicated by Equations (9a[9a] ) and (9b) (Lee and Chen 1984), [9a] [9b]

where the constants A₁, A₂, and A₃ are 1.257, 0.4, and 0.55, respectively (Lee and Chen 1984).

TABLE 2. Simulation parameters

	Cross-shape	Randomly placed particles	DLA	Soot particle simulations
Case names	Case-1	Case-2	Case-3	Case-4
Number of particles	17	500	215	1000
Normalized box length	36σ	40σ	20σ	100σ
NB model	190 σ = 100 nm	—	—	—
SB model	19 σ = 100 nm = 1000	—	—	—
AB model	19 σ = 100 nm 1000 1000	—	—	—
AnBV model	19 σ = 100 nm 1000 1000	19 σ = 50 nm 1 ⩽ ⩽ 1000 10 ⩽  ⩽ 1000	19 σ = 100 nm 500 500	19 σ = 20 nm 500 500

Some of the models that are described in the next section (Section 2.3) require rotationally noninvariant particles. Their orientations can be expressed by quaternions (Martys and Mountain 1999) that are defined in terms of Euler angles for an individual particle as [10] and satisfy the equality [11] The time derivative of the quaternions can then be expressed via the principal angular velocity w_p as [12] This equation for the quaternions is integrated using the same velocity–Verlet integrator as for the solution of the momentum equation. The necessary thermalization of the rotational degrees of freedom is based on LD and is equivalent to the procedure used for the translational degrees of freedom. More details can be found in (Martys and Mountain 1999).

2.3. Particle Agglomeration Models

Two particles are assumed to link permanently and form a cluster when their center-to-center distance, r, becomes smaller than the particle diameter σ. Clusters with at least one common particle constitute an agglomerate. ESPResSo does not explicitly identify agglomerates. We therefore use an external postprocessing tool that performs cluster labeling based on bond information from ESPResSo. Particles are processed serially. If a particle is not connected to a (preceding) particle that had been identified as part of a cluster, this particle forms a new cluster. Otherwise the particle inherits the cluster label of the particle it is connected to. If the particle is connected to more than one particle with different cluster labels, these clusters are merged by assigning a single cluster label to all particles within the clusters.

In this article, four different models that simulate the simultaneous formation and growth of agglomerates, are compared. For the first two models, only the adhesion force between the contacting particles is considered. However, this force acts along the line connecting the centers, and it therefore does not exert torque on the spherical particles. In a further step, we therefore introduce a range of forces and torques, in order to avoid particles from sliding and/or rolling motions. These models are summarized in Table 1 and described in the 4 following subsections.

2.3.1. Model Using Nonbonded Potentials (NB)

In this model, when particles come close, a potential (Equation (2[2a] )) with a highly attractive well (such as as given in Table 2) would make them very unlikely to separate. However, the interactions are still rotationally invariant.

2.3.2. Single Bond Model (SB)

For the “single bond” model, a distance-dependent bond potential (i.e., harmonic bond potential) is used to ensure a binding energy between two contacting particles according to the JKR theory (Chokshi et al. 1913). The binding energy due to the surface deformation is given by [13] where a is the radius of the circular contact region and expressed by [14] Here, Ψ is the surface energy of the particles, is Poisson’s ratio, E_Y is Young’s modulus and D₀ = 0.165 nm (Pantina and Furst 2006). In Equation (14[14] ), the Hamaker constant A is estimated from experimentally measured fragmentation curves for diesel soot (Chokshi et al. 1913; Rothenbacher et al. 2008). Equations (13[13] )–(14[14] ) can thus be used to relate the modeled bond potential in Equation (4[4] ) with the actual physical process by k_h = 2E^JKR_adh/a². The equilibrium distance is set to r₀ = σ which ensures direct particle contact.

2.3.3. All Bonds Model (AB)

The “all bonds” model uses 3-body angular bond potentials to represent tangential forces in addition to the 2-body spring force. In this model, the 2-body potential is used to keep particles together, and the 3-body potentials are responsible for avoiding the sliding/rolling of connected particles. Similar to Becker and Briesen’s study (2008) the angle bond constant k_a (Equation (5[5] )) is estimated using a single-bond bending rigidity, viz. [15] where κ₀ is the bending rigidity per bond and given by (Pantina and Furst 2006) [16]

2.3.4. Angle Bonded Virtual Particles Model (AnBV)

Our goal is that the particles should always be connected at the particular point where they initially touched. This can be achieved by using so-called virtual particles, that is, the particles’ evolution is not subject to the integrated equations of motion, but they are placed at the surface of a given real particle and its relative position to this particle is fixed. This requires that particles carry a co-moving reference frame that is integrated by rotational LD. In our model, the virtual particles are dynamically created during the simulation and move with the Brownian particles they are attached to.

For the “angle bonded virtual particles” model, two virtual particles are created at the point of contact and fixed in the reference frame of the two colliding particles. They are connected to both particles that have collided by angle-dependent potentials with angle θ₀ = π (Equation (5[5] )) as illustrated in Figure 1.

FIG. 1. AnBV model, two virtual particles are placed at the collision point and the angle potential (with θ₀ = π) is defined by two collided particles and virtual particles. Note that we use an exploded view for better illustration.

A 2-body spring bond potential (Equation (4[4] )) with bond length σ between the centers of the particles that collided and a second spring with zero bond length between the virtual sites prevent significant motion around the point of contact, and the angular potential through particles and virtual particles prevents the particles from sliding around each other. This model now is part of the opensource distribution of ESPResSo as a general dynamic bonding feature. Documentation on the creation of virtual particles and dynamic bonding can be found in the ESPResSo user’s guide.

2.4. Simulation Parameters

In this study, the simulations are carried out with a system consisting of spherical particles placed in a cubic box with periodic boundary conditions for all sides. The surrounding fluid is assumed to be air.

We use “reduced units,” in which three reference quantities define the unit system. These are the units of energy (ε*), length (σ*) and mass (m*). All other quantities are expressed in units that can be derived from ε*, σ*, and m*. In the following, the reference values of energy, length and mass are set to the thermal energy of the system, the diameter of the particle, and the mass of the particle, respectively. The normalized key quantities of the LD simulation are then calculated by dividing the real values by the reference units (e.g., , ).

The time step () of the LD simulations plays an important role in solving the equation of particle motion. For the reduction of the computational time it should be as large as possible and preserve accuracy. To obtain reliable simulations, the simulation time step should be much smaller than the relaxation time of the particles, which is the characteristic time for a particle to adjust its velocity. The relaxation time of a Brownian particle is given by [17]

The normalized time step of all test case simulations is set to = 0.01 which is smaller than the relaxation time of the Brownian particles, and air at T = 600K is typically assumed as the surrounding fluid.

The cut-off distance of r_cut = 5σ has been introduced and beyond this distance the attraction between particles is negligible compared to the thermal energy. The nondimensionalized value of the Hamaker constant () corresponds to a dimensional quantity of A = 2.38 × 10^{− 19}J which is a typical value derived for soot particles (Rothenbacher et al. 2008). The corresponding values for the harmonic and the angular bond stiffness constants have been computed from Equations (14[14] )–(16[16] ), assuming the Young’s modulus and Poisson’s ratio of a soot particle to be E_Y = 2 GPa and ν = 0.3, respectively (Bhowmick and Biswas 2011). The binding energy due to surface deformation of the two contacting particles with a diameter of σ = 50 nm is calculated to approximate E_JKR = 1.3 × 10^{− 17} J. This energy is used to estimate the harmonic and angular bond stiffness constants k_h and k_a. This leads to a harmonic bond constant for the SB model of k_h = 0.7 J/m² assuming r_ij = σ − a and r₀ = σ (compare Equations (4[4] )–(14[14] )). The corresponding values for the AB and AnBV models can be estimated to be k_h = 0.35 J/m² and k_a = 0.002 J/m², respectively. This assumes an even distribution of the binding energy to the angular and harmonic bonds. The bond stiffness constants are then normalized by the reference quantities (t^{× 2}/m*). Note that the derived bond constants are rather large and could be reduced by one order of magnitude. This does not affect the simulation results due to the irreversible nature of the bonds but removes unnecessary numerical stiffness.

The other simulation parameters for the different cases are listed in Table 2. For all simulations, the canonical ensemble is assumed, where the number of particles N, the volume V and the temperature T are fixed.

3. RESULTS AND DISCUSSION

In this section, we present the results of the different test cases summarized in Table 2. The first test cases, Case-1 and Case-2, are designed to analyze the capabilities of the different agglomeration models to represent a realistic formation and growth process of agglomerates in dynamic simulations. Case-3 and Case-4 are designed to investigate the cp model in more detail.

3.1. Simulations of the Pre-Defined Agglomerate—Case-1

In this case, the initial positions of particles are given by the shape of a cross. The distance between neighboring particles is closer than the predefined collision criteria in our simulations. Hence, the agglomerates are formed before the iterations start, and the initial shape should be preserved.

Figure 2 shows the configuration of the cross-shaped agglomerate at the end of the simulations using the different models. With the exception of the AB and the AnBV models, the agglomerates do not preserve their initial configuration and change their connectivity. For the NB model, small values of the potential (Equation (2[2a] )) cannot avoid breakage. Larger Hamaker constants yield rather compact agglomerates (D_f = 3) that do not preserve the initial shape.

FIG. 2. (a) The initial configuration of the cross-shaped agglomerate, and its final configurations obtained from (b) NB model with = 190, (c) SB model, (d) AB model, and (e) AnBV model simulations.

We use the longest distance (L_max) between two particles as characteristic quantity to describe shape changes. Figure 3 shows its variation with time for all the models used at different parameters. Independent of the Hamakar constant (A) and bond strength, the NB and SB models cannot preserve the initial configuration of the clusters, which leads to a very noticeable change in L_max. The NB model with a low energy minimum leads to the break-up of the agglomerate and L_max is determined by the size of the simulation box. Larger potential values, however, result in compact aggregates, so that L_max decreases to the distance of around L_max = 1.6σ. In contrast, the change in L_max is small for the AB and AnBV models with angular bond constants , which indicates that the shape is preserved. Even for , L_max decreases only slightly, which is in fact only due to increased flexibility of the side arms, while the overall cross shape is still preserved. However, if the agglomerate would be more compact, this flexibility would allow for more contacts and thus more bonds to form. Bond constants can therefore be used to tune the compactification probability.

FIG. 3. Change in longest distance L_max of initially cross-shaped agglomerates for different simulation models and parameters and averaged over 20 simulations. If the cross shape is preserved, L_max is constant. Thus, the NB and SB models are unable to preserve the shape independently of parameters. In the AB and AnBV model simulations harmonic bond constant is set as , and they require angular bond constant > 10 to preserve the shape.

Figure 4 illustrates the close agreement of the AB and AnBV models. Both models use identical attraction and harmonic bonds to form the bonds. However, for each connection, only one harmonic bond is formed in the AB model, while the AnBV model uses two bonds. Thus, the AnBV model requires just half the bending constant. It is noted here that the dynamic behavior of the two models can differ slightly (as has been observed for the case of randomly placed particles which is introduced in the next section), since the AnBV model gradually introduces additional degrees of freedom, namely the rotational ones, which are thermalized. As a consequence the effective temperature of the AnBV model is slightly higher, but the models can still be matched, if the system temperature is reduced accordingly. While the computational results agree well, the models and their implementations are technically distinctly different. The additional virtual particles and the associated bonds lead to an increase in computational time by a factor of two to three for the computation of the resultant forces on the (real) particles and the respective solution of the Langevin equations. However, the angular bond potential describes interactions between three particles in the AB model, while the interactions in the AnBV model are restricted to two physical particles and their respective virtual particles only. Thus, the contact forces between the particles can be directly accessed, and this will be important when large agglomerates (such as soot agglomerats) are subject to shear and breakage is to be considered.

FIG. 4. Probability density of the deviation (∣π − θ∣) from the equilibrium angle between three beads bonded by the AB and the AnBV models. Note that the AB and AnBV models show the same bending characteristics for and , respectively. This reflects that the AnBV model is constructed using two parallel bonds and therefore requiring only half of the bending rigidity per bond.

FIG. 5. Dynamics of agglomerates obtained from the AnBV model simulations with varying harmonic () and angular () bond constants. (a) Change in the number of agglomerates, and (b) average fractal dimensions of the systems.

3.2. Simulations of Randomly Placed Particles—Case-2

Further studies have been performed to investigate the AnBV model’s characteristics in more detail and to identify the differences between the AB and AnBV models. The fractal dimension, D_f, is now used as a measure of the agglomerates’ structure. It is measured from the slope of the double logarithmic plot of the radius of gyration, R_g, versus the number of particles, N (Xiong and Friedlander 2001). This definition is consistent with the common interpretation of a dimension if the object is compact, that is, the fractal dimension of a sphere is D_f = 3, of a disk is D_f = 2 and of a straight chain is D_f = 1. Structures that are not compact can be characterized by fractal dimensions.

The simulations of Case-2 are used for a quantitative and qualitative assessment of the structure of agglomerates that are formed from primary particles subject to Brownian motion. The agglomerates are formed by the collisions of single Brownian particles and/or collisions of clusters during the simulation. Since the clusters are formed by bonded particles, they move naturally following the same unshielded Langevin equation of motion as the nanoparticles (Isella and Drossinos 2011). A total of 500 particles have been placed randomly in a box with periodic boundary conditions. These values represent reference values for our sensitivity studies with respect to the bond constants in Sections 3.1 and 3.2.

The growth and fractal dimension of formed agglomerates (N > 15) are shown in Figures 5a and b, respectively, for different and values. At the beginning of the simulations the particles collide and create clusters as indicated by the rapid increase of the number of agglomerates. This rapid increase is followed by a much longer period of reducing the number of agglomerates due to particle-cluster and cluster–cluster collisions that form larger and larger clusters. The dynamics of cluster formation are only little affected by the bond constants chosen for the respective computations.

The agglomerates’ average structures, however, are much more sensitive to the choice of bond constants as illustrated in the same figure and quantified by the different values of D_f. The reduced bond constants provide more flexibility, and the agglomerates can bend, fold, elongate and/or shrink during the simulation. The mean fractal dimension varies between D_f = 1.73 and D_f = 2.22 for the simulations with smaller bond constants (i.e., and ). This agrees well with the numerical studies by Chen et al. (1988) and Jullien and Meakin (1989) for cluster–cluster aggregates (CCA) grown by DLA including impact restructuring. In their studies, the aggregates might fold, bend and twist. The fractal dimension was found to equal D_f = 2.09 if the agglomerates were allowed to bend, D_f = 2.17 if both bending and folding were allowed and D_f = 2.19 if the complete restructuring stages were included for CCA grown by DLA.

In contrast, the mean fractal dimension varies between D_f = 1.65 and D_f = 1.71 for the AnBV model simulation with higher bond constants (i.e., = 1000 and = 1000). This indicates that no significant restructuring occurs and results match well with the measured fractal dimensions of the CCA obtained by various numerical and experimental studies (Bourrat et al. 1987; Samson et al. 1987; Zhang et al. 1988; Megaridis and Dobbins 1990). It is thus apparent that the correct choice of bond constants allows for the representation of different macroscopic interactions between the particles at the cp. The effect of the different bond constants can also be quantified in terms of changes in fractal dimension and number of cp with time as listed for some representative agglomerates in Table 3.

TABLE 3. Time sequence for changes in D_f and in the number of cp of representative agglomerates, A_i, during the AnBV model simulations of Case-2. The harmonic and angular bond constants are given as: A₁ (=1 =10), A₂ (=10 k_a=10), A₃ (=1000 =10), A₄ (=10 =1000), A₅ (=1000 =1000). N gives the number of primary particles in each cluster

	A1(N = 136)		A2(N = 130)		A3(N = 140)		A4(N = 124)		A5(N = 136)
	Df	cp	Df	cp	Df	cp	Df	cp	Df	cp
25	1.75	533	2.0	524	2.0	433	1.75	245	1.65	280
27	1.77	535	2.17	527	2.0	434	1.44	245	1.65	280
29	1.5	537	2.15	528	2.0	436	1.41	245	1.65	280
31	2.19	556	2.15	528	2.0	437	1.54	246	1.65	280
33	1.94	563	2.16	528	2.0	440	1.48	246	1.65	280
35	2.16	568	2.16	529	2.0	442	1.43	246	1.65	280

The computations demonstrate that the resulting morphology depends largely on the harmonic bond potential . Small values, that is, allow bending of the bond and thus lead to much more compact cluster structures as time advances (agglomerate A₁) while large values preserve the gobal structure (agglomerate A₃). In addition, high values of support preservation of the agglomerates structure by preventing the particles to fold and thus to create new cp (agglomerates A₄ and A₅ with and a nearly constant number of cp). Note that both and thus allow to control the flexibility and by that the compactness of the agglomerates. However, reducing facilitates a rolling motion which may or may not be physical and will be dependent on the properties of the investigated system/agglomerates. Reducing , however, allows the primary particles to separate temporarily, which is unphysical under van der Waals forces. Therefore, reducing should be preferred and should be kept at its physically determined value of .

FIG. 6. Comparison of the size of the aggregates characterized by radius of gyration, R_g during the growth process. Stars (red) and solid line represent the conventional model, squares (blue) and dotted line represent the AnBV model. The slopes of the linear fits correspond to the average fractal dimensions of the agglomerates.

3.3. Diffusion Limited Aggregation—Case-3

The third test case (Case-3) is the simulation of the well-known DLA that can be used to further validate the developed AnBV model. We study DLA of the irreversible growth of a single cluster grown from a seed particle fixed in three-dimensional space. The growth rule is simple, in which the first particle is fixed in the center of the system, and the following particle is then released from a random position far away and is allowed to move due to Brownian motion. If it collides with the first particle, it is connected according to the rules in our AnBV model algorithm and becomes part of the agglomerate. Then, further particles are launched one-by-one and each of them is connected to the cluster after hitting any of the particles belonging to the cluster.

The growth mechanisms are same for the both cases and explained above. However, in the “conventional DLA” model, when the released particles hit any particle of the cluster, they stop as in the classical DLA model, that is, the formed clusters are immobile. In the “AnBV DLA” model, the aggregate itself is moving, and therefore increasing the probability of collision. The size of the agglomerates are characterized by the radius of gyration, R_g that is given by [18a] with [18b] being the agglomerate’s center of mass. Figure 6 compares log (R_g) plotted versus log (N) for the two models. The fractal dimensions of the DLAs can be obtained from linear fits of the data as D_f = 2.386 and D_f = 2.392 for the AnBV and the conventional DLA models, respectively. This confirms our expectation that the motion of the cluster has no significant influence on the DLA.

3.4. The Growth of Soot Agglomerates—Case-4

For a final validation we compare results from our LD simulations with soot aggregates generated in ethylene diffusion flame and measured by Megaridis and Dobbins (1990). They collected samples using a thermophoretic sampling technique by transmission electron microscopy and determined an average fractal dimension of D_{f, exp} = 1.69. Our LD shall mimic the experimental conditions and we used primary particle diameters of σ = 40 nm and the molecular mean free path length is set to l_fluid = 600 nm. The computations show average fractal dimensions of D_f = 1.68 and analysis is based on all aggregates with more than 15 particles. The maximum agglomerate composes n = 157 primary particles. The numerical value agrees very well with the value given in Megaridis and Dobbins (1990) and a qualitative comparison of the shape of the agglomerates (Figure 7) demonstrates the realistic structures of the simulated soot aggregates and supports the need for a novel implementation such as the AnBV model for the simulation of aggregation where particles do not move relative to each other after collision.

FIG. 7. Clusters at the end of the LD simulations leading to agglomerates with an average fractal dimension of D_f = 1.68.

4. CONCLUSIONS

In this study, we studied the capability of four different models to simulate the nanoparticle agglomeration process. All four models require local bonding only. Analyzing the shape and structure of both pre-defined and randomly formed agglomerates in terms of fractal dimension, it became evident that only models with noncentral interactions, the so-called AB and AnBV models, are capable of modeling nanoparticle agglomeration without undergoing any restructuring process. Using only nonbonded central interactions (NB model) is not sufficient to hold particles together and to prevent sliding around the cp of individual particles. Also, forming only harmonic bonds on contact, the SB model, does not fix the particles at their cp and allows restructuring. Hence, only bulky and compact agglomerates are formed, as one would expect at high process temperatures. We have demonstrated that both the AB and AnBV model can prevent restructuring with tunable degree of compactification. Moreover, the AB and AnBV models can be parameter matched to exhibit identical bending, and thus restructuring properties. However, technically the models are very different. The AB model is computationally faster, but only the AnBV model gives direct access to shear forces, which is required if breakage has to be included. Finally, the AnBV model’s ability to reproduce diffusion limited particle–cluster and cluster–cluster agglomerates has been shown. The simulations show remarkable agreement with experimentally determined fractal dimensions of soot aggregates, and we may conclude that the AB and AnBV models provide a suitable implementation of nanoparticle aggregation based entirely on local interactions.

FUNDING

The authors acknowledge the financial support by the Deutsche Forschungsgemeinschaft (DFG) as part of the Collaborative Research Center (SFB) 716 held by the University of Stuttgart.