**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 Ω*<sup>x</sup>* subdomain, the dimensionless weak form of the energy equation is given by the relation [43]:

$$\left(k\_{\rm rp} - 1\right) \int\_{\partial\Omega\_{\rm x}} \Phi \nabla T\_{\hbar} d(\partial\Omega\_{\rm x}) + \int\_{\partial\Omega\_{\rm x}} \nabla T\_{\hbar} d(\partial\Omega\_{\rm x}) = 0,\tag{6}$$

where Φ is a spatial step function defined as unit in the particle phase and zero elsewhere, and *kr p* = *kp k f* 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 *ke f f* = R *S k ∂T ∂n dS*, where the surface *S* 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 ~5 mil. nodes and ~1.5 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, *ke f f* , of a suspension of solid spheres in liquid [51]:

$$k\_{eff} = \frac{k\_p + 2k\_f + 2\left(k\_p - k\_f\right)f\_p}{k\_p + 2k\_f - 2\left(k\_p - k\_f\right)f\_p}k\_{f\prime} \tag{7}$$

where *k <sup>f</sup>* is the conductivity of the base fluid, *f<sup>p</sup>* is the volume fraction of the particles, and *k<sup>p</sup>* 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 *k<sup>p</sup>* with *k <sup>f</sup>* and *f<sup>p</sup>* with 1 − *f<sup>p</sup>* , the effective conductivity calculation changes drastically [11]:

$$k\_{eff} = \frac{k\_f + 2k\_p + 2\left(k\_f - k\_p\right)\left(1 - f\_p\right)}{k\_f + 2k\_p - 2\left(k\_f - k\_p\right)\left(1 - f\_p\right)} k\_{p\prime} \tag{8}$$

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 (*ke*) and an effective volume fraction (*fe*). For the calculation of *k<sup>e</sup>* , many relations from the effective medium theory have been applied [11,52–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 (*fe*), can be expressed as:

$$f\_{\varepsilon} = \frac{4\pi}{3} \sum\_{i}^{N\_{\varepsilon}} R\_{\mathbb{S}^{I}} \,^3 \prime \,\tag{9}$$

where *N<sup>c</sup>* is the number of aggregates in the solution and *Rgi* is the (dimensionless) radius of gyration of each aggregate. The volume fraction of the solid phase inside the aggregates is *f<sup>i</sup>* = *fp*/ *f<sup>e</sup>* [17]. In this work, *k<sup>e</sup>* is calculated from the upper limit (Equation (8)) by swapping *f<sup>p</sup>* with *f<sup>i</sup>* , while the Maxwell relation is used (Equation (7)) in a second step, by replacing *k<sup>p</sup>* with *k<sup>e</sup>* and *f<sup>p</sup>* with *f<sup>e</sup>* .
