1. Introduction
Numerous contemporary applications are based on the incorporation of nanoparticles into convectional fluids [
1,
2]. The resulting nanofluids may have drastically increased thermal properties and reduced sedimentation. Nanofluids and nanoparticles are increasingly used in a variety of fields. Heating and cooling systems using nanofluids show significant improvement in energy consumption, while new medical techniques are developing using nanoparticles in bioliquids. [
3,
4]. A large number of recent publications study the potential use of nanofluids in multidisciplinary fields [
5,
6]. Although much effort has been placed on the correlation of the heat transfer properties of nanofluids with the underlying phenomena, it appears that there is no widely accepted explanation of their behaviour or a reliable way to predict their heat conduction properties [
7,
8].
Many models have been developed to predict the effective conductivity of nanoparticles. The effect of Brownian motion, interfacial resistance, the existence of nanolayers, and the aggregation mechanism have been discussed in detail in the literature [
9,
10,
11,
12,
13]. A large increase in effective thermal conductivity has been detected experimentally when nanoparticles are organized in small aggregates [
14,
15]. The increased contact of the particles within the aggregate was found to facilitate heat transfer compared to fully dispersed particles. On the contrary, larger mass aggregates have a negative effect on the stability of the nanofluid and, therefore, on heat transport properties [
16]. Lotfizadeh et al. [
11], Prasher et al. [
17], Evans et al. [
18], and Liao et al. [
19], among others, developed models to predict the thermal conductivity of nanofluids based on the morphology of the aggregates. These works showed that the configuration of the nanoparticles and the morphological parameters of the aggregates can alter the effective conductivity of the nanofluids noticeably. A typical assumption made in several studies is that the nanoparticles are of the same size and the neighbouring monomers are mainly at single point contact.
In real conditions, samples always come with a certain distribution in particle size [
20,
21,
22], which is held responsible for altering the thermal properties of nanofluids [
21]. It has been noticed that polydispersity may occur during dispersion of the nanoparticles in the base fluid [
23,
24]. Chon et al. [
21] used commercial, uniform nanoparticles for the preparation of nanofluids. They measured the size distribution after dispersion and found a significant deviation. In fact, it is technically challenging to synthesize and disperse a large quantity of highly monodispersed nanoparticles [
25]. Zhiting et al. [
26] studied the effect of polydispersity, among other parameters, on the heat transfer coefficient of nanocomposites with molecular dynamics simulations. They concluded that polydispersity negatively affects the effective conductivity. However, the nature of the used method does not allow the simulation of large systems, such as that of aggregated nanoparticles.
Strong electrostatic forces between particles, collision of particles during the formation of the aggregates, and high-temperature environments are some of the factors that contribute to a certain degree of overlapping between particles [
27,
28]. Two methods have been widely used for the preparation of nanofluids. The one-step nanoparticle production and blending method produces nanofluids with increased stability and offers elevated thermal conductivity, but has a relatively high production cost and is not yet suitable for large-scale production [
29]. The two-step method is one that proceeds in two sequential steps, namely, the step of separate production of nanoparticles followed by suspension in a base fluid [
29]. Commercial nanoparticles are usually found in powder form. The production of nanoparticles commences with the creation of a nanoparticle suspension from their precursors. This suspension is dried using various methods and, eventually, a powdered form of nanoparticles is obtained [
23,
30]. Thermal decomposition of organic precursors is a well-established process for the fabrication of solid nanoparticles [
23]. It is temperature dependent and has been reported to result in sintered and/or polydispersed particles [
23,
24,
30]. Other methods include nanoparticle production in the gas phase [
31,
32]. During the creation of these nanoparticles, formation processes, including surface reactions, condensation, coagulation and sintering, are some of the key mechanisms that take place [
32,
33,
34,
35]. The kinetics of each process determines the final structure morphology, which can vary among spherical particles, agglomerates, or compact aggregates [
34]. The resulting system usually includes partially coalesced particles with sintering necks [
32,
35].
Attempts to model the aforementioned process usually start with the formation of an aggregate using a stochastic method. Eggersdorfer et al. [
36] modelled the sintering process in aggregates, formed by Diffusion Limited Aggregation (DLA), Diffusion Limited Cluster–Cluster Aggregation (DLCCA), and Ballistic Aggregation (BA). The driving force for sintering was the minimization of free energy. They noted that, during sintering, primary particles approach each other. Sander et al. [
31] presented an analytical description of the underlying phenomena during the production of nanoparticles, such as coagulation and sintering. The primary particles were modelled having a spherical shape and a polydispersed size, with each particle described by its sintering level and radius. The final structure has been compared with experimental data and transmission electron microscopy (TEM) images. An overlapping algorithm, developed by Brasil et al. [
27], has studied the effect of sintering on the morphological properties of the aggregates, such as the fractal dimension and the radius of gyration. Schmid et al. [
33,
37] have developed a model for aggregates subjected to coagulation and sintering. The sintered aggregates were presented as the result of successive overlapping of spherical, primary particles.
The previous discussion underlines two major open issues in the study of the thermal conductivity of aggregated nanoparticles. Even though the effect of aggregation has been extensively studied, there is a dearth of research dealing with the effect of polydispersity of the nanoparticles within the aggregate. Moreover, sintering and the concomitant partial coalescence are most likely to occur in nanoparticle systems, yet the study of their effects on heat conduction remains a challenging, open field. These configurations complicate the determination of heat transfer properties and, therefore, reliable simulations of transport phenomena are required.
The present work examines the effects of polydispersity and sintering of particles on the effective thermal conductivity of nanofluids that contain particle aggregates. To this end, the method developed by the authors [
38] for reconstructing particle aggregates is extended to include particle overlapping due to sintering as well as non-uniform particle size as a realistic outcome of nanofluid preparation. Among the merits of this method is the algorithmically rapid reconstruction of agglomerated systems with predetermined properties, namely, the fractal dimension and the average number of particles in the aggregates. As a case study, the particle size here follows the normal distribution and the standard deviation is expressed as a fraction of the mean size. Moreover, a technique has been developed to simulate sintered aggregates. The sintering process is expected to change the particle position and size while, naturally, preserving the mass of the working sample. An overlap parameter and the morphology of the primary aggregate determine the final morphology. The effective thermal conductivity is calculated through the temperature distribution obtained from the solution of the heat transfer equation. The Meshless Local Petrov–Galerkin (MLPG) method [
39,
40,
41] is used here as it was shown to provide stable and fast solutions to particulate systems even with point contact. The Discretisation-Corrected Particle Strength Exchange (DC PSE) method [
42,
43] is used to approach the field function and its derivatives, while the meshless nature of the method allows local increase of the domain discretisation at the interface between the base fluid and solid particles.
The effect of the overlap parameter and the polydispersity level of particle aggregates on the thermal conductivity is studied by changing the number of particles in the aggregate, the fractal dimension of the aggregates, and the volume fraction of the particles. The effective conductivity of the polydispersed nanoaggregates, as predicted by the present method, is compared to the effective conductivity of the corresponding systems, as these result from the Diffusion Limited Aggregation (DLA) method. Moreover, the effective conductivity of aggregates consisting of polydispersed particles is compared with that of aggregates of monodispersed particles, keeping all other morphological parameters constant. Notable deviations between monodispersed and polydispersed cases are observed and discussed. In addition, the effect of sintering is examined by varying the overlap parameter. The results are compared with predictions of analytic expressions from the literature.
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,
, using the relation between the number of particles in the aggregate,
, and basic cluster-size characteristics [
45,
46]:
N.
where
is the mean radius of the primary particles,
is the structure factor, and
is the radius of gyration of the aggregate [
12,
32,
47]:
where
is the position vector of the centre of mass of particle
,
is its mass, and
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
. All spatial parameters and variables are normalized with this quantity. Thus, the particle radius is related with the volume fraction (
), the number of particles in the aggregate (
), and the number of aggregates (
), as follows:
For monodispersed particles, the definition of
is straightforward. For polydispersed particles, in the present extended approach, this value is set as the mean radius of the particles,
. 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
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
. The maximum deviation of the particle radius, in this study, is
. 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.
The structure factor,
, 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,
[
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]:
where
is the final distance of the centres of the neighbour particles, and
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,
, while the sizes of the particles remain constant. The penetration coefficient relates the final distance of the neighbouring particles to the initial radius (
). 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 (
). 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:
The overlapping coefficient, δ, is the sum of the penetration () and the growth () coefficients: . If the aggregates are totally sintered (i.e., every aggregate merged into a single particle), whereas indicates that particles are in point contact. The permissible values for the penetration () and growth () 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
and
. Evidently, each combination results in a different volume fraction (
). In order to determine the appropriate combination, for each aggregate the value of the penetration coefficient varies from zero to
. For each
value, the final radius (Equation (5)) is calculated. Finally, the volume fraction of the resulting system is compared with the initial volume fraction (
), (
) and its dependence on
, as shown in
Figure 2. The
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 (
), for different values of the overlapping coefficient and the morphological characteristics of the initial aggregate. It is shown that a unique combination of
and
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.
3. Effective Conductivity Calculation
The reconstructed aggregates that are obtained following the algorithm of the previous section are used as input to heat transport modelling. A constant temperature difference is imposed along the vertical axis of the nanofluid, whereas the rest of the boundaries are considered periodic. For the numerical solution of the heat transport equation, the Meshless Local Petrov–Galerkin (MLPG) method is used [
41]. It has been shown to offer concrete advantages to more conventional methods in particulate systems with several contact areas, as is the case here. Differential equations are integrated into local subdomains, a fact that facilitates the increase of discretisation in regions of the geometries where steep gradients are developed. In the nanofluid case, the effective conductivity changes drastically at the interface of nanoparticles and the base fluid. The shape of the subdomains alters the performance of the method, with cubic sectors having been proved to increase the stability of the method [
40]. The DC PSE approach is chosen as the trial function, while a step function is used as test function for the integration [
42,
43]. A set of cubic grids digitizes the domain and the integrals are calculated with the Gauss quadrature method [
43]. In each
subdomain, the dimensionless weak form of the energy equation is given by the relation [
43]:
where
is a spatial step function defined as unit in the particle phase and zero elsewhere, and
is the ratio of the conductivity of the particles to that of the base fluid. The solution of the heat transfer equation determines the temperature throughout the computational domain. Then, the calculation of the dimensionless effective conductivity is straightforward from
, where the surface
is vertical to the heat flow. A detailed description of the approach, the respective variables and integrals, the mesh construction, and the conductivity calculation can be found in a previous work by the authors [
43].
This method is capable of calculating the effective conductivity of large particle systems. The aggregates are considered stationary and the heat conduction equation is solved within the computational domain. A typical simulation contains about 500–1000 particles organized into aggregates. Modelling of aggregates in heat transfer processes is performed with in-house meshless CFD methods implemented in Matlab kernels, as described in [
38]. The computational time for the reconstruction of the aggregates is also provided in [
38], along with comparison with other aggregation models. The effective conductivity calculations used herein have been shown in [
43] to reduce the computational cost, compared with other numerical models and commercial software. A typical run for 1000 particles requires
mil. nodes and
ks on an Intel(R) Xeon(R) Silver 4116 CPU at 2.10 GHz using 12 cores.
Moreover, the present method can be extended to include calculations of the effective conductivity of different nanoparticle shapes. If the equation of the external surface of the particles is simple, the application is straightforward; otherwise the aggregation algorithm should be modified rather drastically, especially for non-convex surfaces.
The corresponding predictions of analytical models for the effective thermal conductivity are presented next. A well-known model for the conductivity of dispersed particles is Maxwell’s effective medium theory. Maxwell developed an expression for the effective conductivity,
, of a suspension of solid spheres in liquid [
51]:
where
is the conductivity of the base fluid,
is the volume fraction of the particles, and
is the conductivity of the particles. According to this equation, the process is apparently not sensitive to the size or the arrangement of the dispersed phase. However, this relation is not symmetric. By replacing
with
and
with
, the effective conductivity calculation changes drastically [
11]:
Equation (7) accurately predicts well-dispersed particles at low volume fraction, while Equation (8) has been used to describe the conductivity of aggregate structures. In this perspective, the particles are considered as a solid network, with the base fluid enclosed in some regions. In any case, Equations (7) and (8) estimate the lower and upper bounds of the conductivity of an inhomogeneous medium, respectively [
52,
53].
The majority of the attempts to develop a model for the conductivity of colloidal clusters include a two-step approximation. The clusters are considered as spheres, with an effective conductivity (
) and an effective volume fraction (
). For the calculation of
, many relations from the effective medium theory have been applied [
11,
52,
53,
54]. The size of these effective particles is usually set equal to the radius of gyration of the aggregates [
11,
55]. In this case, the effective volume fraction (
), can be expressed as:
where
is the number of aggregates in the solution and
is the (dimensionless) radius of gyration of each aggregate. The volume fraction of the solid phase inside the aggregates is
[
17]. In this work,
is calculated from the upper limit (Equation (8)) by swapping
with
, while the Maxwell relation is used (Equation (7)) in a second step, by replacing
with
and
with
.
5. Conclusions
The effect of particle size polydispersity and the sintering level on the thermal conductivity of aggregated nanoparticles was studied in the present paper. It was shown that both parameters examined here have the potential to change the heat performance of nanofluids drastically.
A method for reconstructing aggregates with the desired polydispersity degree was developed, satisfying simultaneously the requirements for certain morphological characteristics of the aggregate, namely, the fractal dimension and the number of particles in the aggregate. Particle sintering in aggregates was simulated for monodispersed cases and encoded as an overlapping mechanism in two steps: a penetration step and a growth step. In order to ensure mass conservation, the progression of each step was controlled through the minimisation of the error in the volume fraction of the sintered aggregate compared with the volume fraction of the initial aggregate. A meshless method with local refinement was used for the solution of the heat transfer equation and was found to be stable for the complex systems that were studied here. This is of key importance in the present problem as it allows using relatively large working domains that contain overlapping particles or particles at point contact with others and being able to extract statistically meaningful conclusions.
The effective thermal conductivity was calculated for aggregates that resulted from the present method of aggregation of polydispersed particles, then compared with the thermal conductivity of aggregates that were constructed with the Diffusion Limited Aggregation (DLA) method. The dependence of the effective thermal conductivity on the fractal dimension was found to be in good agreement with that in DLA method aggregates, which indicates that the proposed method produces aggregates that are thermally equivalent to those resulting from methods that describe the physical process of particle aggregation. Consequently, one can employ the present method for the investigation of the behaviour of nanofluids in heat transport problems, taking advantage of the increased simplicity of the aggregation algorithm and its rapid convergence to the final configuration.
The variation of the effective thermal conductivity was investigated over a wide range of fractal dimension values, number of particles per aggregate, and standard deviation of the particle size. Compared to fully dispersed particles, aggregation was shown to increase the thermal conductivity in all cases studied here. Small radius deviation does not substantially change the thermal conductance compared to monodispersed cases; however, a further increase of polydispersity leads to a clear reduction of the effective conductivity. On the contrary, strong polydispersity leads to an increase in the projected area, which implies an increase of heat transfer. A possible explanation for our result could be the existence of small particles within the aggregate that hinder heat transfer. This result is qualitatively confirmed by experimental measurements [
56,
57] according to which nanofluids consisting of particles with low polydispersity levels have higher heat performance compared to particles with high polydispersity.
The two-step Maxwell model predicts a monotonic decrease of the effective conductivity with increasing fractal dimension; however, large deviations from the numerical results were found for most of the cases examined here.
The effect of sintering of the aggregates was investigated and quantified as a function of the overlapping coefficient. Sintered aggregates have a lower effective size than the original aggregates, so a reduction in the effective conductivity should be expected. At the same time, sintering increases the heat conduction by forming larger heat pathways. This interplay yields a maximum in the thermal conductivity as a function of the degree of coalescence. The precise value of the overlapping coefficient that provides the highest conductivity increase depends on the morphological properties and the volume fraction of the initial aggregates. The present study indicates that the conditions of the production and dispersion of nanoparticles have a major impact on the thermal properties of the nanofluids. This is a possible explanation for the large deviations that have been observed between experimental works. Nanofluids with monodispersed particles, which are organized into aggregates with small overlapping, offer the highest heat transfer coefficient over the range of parameter values that were examined here. The results and conclusions of this work are also relevant to nanocomposite materials that contain polydispersed particle inclusions, which are organized in aggregates either at simple contact or in sintered form.