**2. Modelling Polydispersed and Sintered Aggregates**

Modelling of fractal aggregates is the first stage in the process of correlating numerically their physical properties with their morphological parameters. In reality, aggregation is a very complex process, being sensitive to parameters such as the temperature, the physical properties of the particles and the solvent, the polydispersity extent, the primary particle shapes, etc. [11,12]. A well-established approach to generate fractal structures numerically is the usage of various stochastic methods, such as Diffusion Limited Aggregation (DLA), Diffusion Limited Cluster–Cluster Aggregation (DLCCA), and Reaction Limited Aggregation (RLA) [43,44]. The validity of these methods has been verified through comparison of the resulting structures with experimental data [45,46]. Different physical processes are engaged in different aggregation models. A very common description of such clustering relies on the fractal dimension, *d<sup>f</sup>* , using the relation between the number of particles in the aggregate, *N*, and basic cluster-size characteristics [45,46]: *N*.

*N* = *k<sup>g</sup> Rg*/*r<sup>p</sup> d<sup>f</sup>* , (1)

where *r<sup>p</sup>* is the mean radius of the primary particles, *k<sup>g</sup>* is the structure factor, and *R<sup>g</sup>* is the radius of gyration of the aggregate [12,32,47]:

$$R\_{\mathcal{S}} = \sqrt{\frac{\sum\_{i}^{N} m\_{i} (r\_{i} - r\_{c})^{2}}{\sum\_{i} m\_{i}}},\tag{2}$$

where *r<sup>i</sup>* is the position vector of the centre of mass of particle *i*, *m<sup>i</sup>* is its mass, and *r<sup>c</sup>* is the position vector of the centre of mass of the aggregate.

In order to study the effect of polydispersity on thermal conductivity, the work developed in [38] has been extended to include polydispersed particles. This technique offered fast convergence of the algorithm for the representation of agglomerated systems with predetermined properties. For the sake of completeness, the major steps of the algorithm are mentioned below. The primary input of the algorithm consists of the volume fraction of the particles, the fractal dimension or a range of values around it, and the average number of particles per aggregate or a range of values around it. A random deposition of a particle initiates the process. A new particle stochastically appears on the surface of the particle, and the process is repeated, with the restriction of no overlapping between any pair of particles. During the process, certain restrictions are imposed as described in [38], aiming at the convergence of the fractal dimension to the target value or range of values. It has been shown that the desired fractal dimension can be achieved with only a few particles. The process ends when the aggregate acquires a predefined number of particles. Then, another particle appears at a random place in the computational domain, and the aforementioned process is repeated for the formation of a second cluster. The whole algorithm is repeated as many times as needed to satisfy the desired number of aggregates in the working domain. Further analytical descriptions of the method are presented in [38].

In the present study, all simulations take place in a cubic box of length *l*. All spatial parameters and variables are normalized with this quantity. Thus, the particle radius is related with the volume fraction (*fp*), the number of particles in the aggregate (*N*), and the number of aggregates (*Nc*), as follows:

*r<sup>p</sup>* = <sup>3</sup> s 3 *f<sup>p</sup>* 4*π* ∑ *Nc <sup>i</sup> N<sup>i</sup>* , (3)

For monodispersed particles, the definition of *r<sup>p</sup>* is straightforward. For polydispersed particles, in the present extended approach, this value is set as the mean radius of the particles, *r*0. The exact radius of each particle is randomly sampled from a prescribed normal distribution. The deviation of the particle size distribution is expressed as a fraction of the particle radius. The probability density function (*p*) for the particle radius is shown in Figure 1 for two different standard deviations. In order to avoid negative values, a threshold was imposed at zero radius. For symmetry reasons, another threshold is set at the value *rp*,*max* = 2*r*0. The maximum deviation of the particle radius, in this study, is *σ* = 0.5*r*0. For higher polydispersity levels a lognormal distribution can be used instead, in order to avoid a large number of negative values and adhere to more realistic particle size distributions.

The resulting system may have a volume fraction different from the desired one. This is an important issue in polydispersed particle systems since large sizes may be sampled as the tail additions to a cluster. Depending on whether the desired volume fraction is smaller or greater than expected, particles will be added or removed. To remove a particle, an aggregate is randomly selected and the last-added particle is removed. To add a particle, an aggregate is randomly selected and a new particle with a size that is sampled from the prescribed distribution is added at a random location on the surface of a randomly selected particle of the aggregate, with restrictions in order to satisfy the fractal dimension and the non-overlapping condition, as described in [38] for uniformly sized particles. The process is repeated until the predetermined volume fraction is achieved. With this methodology, critical quantities such as the volume fraction, the mean radius of the particles, and the fractal dimension remain within their prescribed bounds while a small deviation is maintained in the number of particles per aggregate.

**Figure 1.** The probability density function, *p*, of the particle radius distribution for different standard deviation values (red *σ* = 0.2*r*0, blue *σ* = 0.5*r*<sup>0</sup> ) and a representative visualization of the resulting aggregates.

The structure factor, *kg*, is strongly dependent on the polydispersity levels and the degree of overlapping of the particles [48,49]. Eggersdorfer et al. [49] and Tomchuk et al. [48] studied the effect of polydispersity on the fractal dimension and the structure factor, considering a wide range in the number of particles in each aggregate, the prescribed deviation of the particle size, and the aggregation model that is used. They noted a reduction in the structure factor for increased polydispersity for aggregates formed by the DLCCA algorithm and the Ballistic Aggregation model. Independently of the agglomeration mechanism, in the limiting case of infinitely polydispersed particles, the structure factor tends to unity [49]. For monodispersed particles with fractal dimension ranging between 1.7 and 2.5, the structure factor can be considered constant, *k<sup>g</sup>* = 1.5 [38,50]. In the present case, for polydispersed particles, the structure factor changes linearly with the deviation of the particle radius, taking values between 1.2 and 1.5 [49].

During sintering, particles are expected to increase their radius and come closer to each other [37]. Typical simulations of coagulation and sintering of nanoaggregates include an overlapping step, where neighbouring particles penetrate each other, and a growth step, during which particle size increases to maintain mass and volume. This process captures the redistribution of mass in the free surface of the aggregate and offers a realistic representation of the final morphology [33]. Assuming aggregates consisting of monodispersed particles and following this methodology, an overlapping coefficient is defined as [27]:

$$
\delta = 1 - \frac{d}{2R\_p} \, , \tag{4}
$$

where *d* is the final distance of the centres of the neighbour particles, and *R<sup>p</sup>* is the final radius of the particles. The initial radius of the particles can be calculated from Equation (3).

At initial stages of the sintering algorithm, aggregates are forced to collapse to their centre of mass by the penetration coefficient, *δp*, while the sizes of the particles remain constant. The penetration coefficient relates the final distance of the neighbouring particles to the initial radius (*δ<sup>p</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>d</sup>* 2*rp* ). Obviously, this process causes a mass loss. In a second step, the particle sizes are increased, in order to reproduce the volume fraction at its prescribed value. The growth coefficient relates the initial and the final particle radius to the final distance of the neighbouring particles (*δ<sup>R</sup>* = *Rp*−*r<sup>p</sup> Rprp d*/2). Combining these definitions, the final radius of the particles can be related to the initial radius, the overlapping coefficient, and the penetration coefficient, as follows:

$$R\_p = \frac{r\_p}{1 - \delta} (1 - \delta\_p) \tag{5}$$

The overlapping coefficient, *δ*, is the sum of the penetration (*δp*) and the growth (*δR*) coefficients: *δ* = *δ<sup>p</sup>* + *δR*. If *δ* = 1, the aggregates are totally sintered (i.e., every aggregate merged into a single particle), whereas *δ* = 0 indicates that particles are in point contact. The permissible values for the penetration (*δp*) and growth (*δR*) coefficients range from zero to *δ*.

The generation of sintered aggregates initiate with the determination of the volume fraction, the number of particles in the aggregate, the fractal dimension, and the overlapping coefficient. After the formation of each aggregate, a series of trial simulation scenarios are used to evaluate the values of *δ<sup>p</sup>* and *Rp*. Evidently, each combination results in a different volume fraction (*fp*,*<sup>δ</sup>* ). In order to determine the appropriate combination, for each aggregate the value of the penetration coefficient varies from zero to *δ*. For each *δ<sup>p</sup>* value, the final radius (Equation (5)) is calculated. Finally, the volume fraction of the resulting system is compared with the initial volume fraction (*fp*), (*β* = 100· *fp*,*<sup>δ</sup>* − *f<sup>p</sup>* / *fp*,*<sup>δ</sup>* ) and its dependence on *δp*, as shown in Figure 2. The *δ<sup>p</sup>* value with the smallest acceptable error is selected. Following the aforementioned technique, mass conservation is secured for each aggregate of the system, for the entire range of the overlapping coefficient. In Figure 2, the percentage error in volume fraction (*β*) is represented as a function of the penetration coefficient (*δp*), for different values of the overlapping coefficient and the morphological characteristics of the initial aggregate. It is shown that a unique combination of *δ<sup>p</sup>* and *R<sup>p</sup>* results in a system with the same volume fraction. This methodology, in addition to achieving the desired overlapping coefficient and volume fraction, has the advantage of being straightforward and fast during calculations. Needless to say, in real conditions the final structure during sintering may be different from that of the overlapping spheres, due to the appearance of neck effects and redistribution of mass that will eventually differentiate the structure from the one that is simulated here.

**Figure 2.** The percentage error in volume fraction (*β*) varying with the penetration coefficient for different values of the fractal dimension, the number of particles in the aggregate, and the overlapping coefficient. A representative visualization of each aggregate structure is also shown.
