Green–Kubo assessments of thermal transport in nanocolloids based on interfacial effects

Abstract Thermal transport in a water–Cu nanocolloid system was investigated using equilibrium molecular dynamics. A systematic analysis of the Green–Kubo calculations is presented to clarify the effect of simulation parameters. Several sources of error were identified and quantified for the thermal conductivity estimations, and the effect of the base fluid potential was investigated. Simulations were carried out with a single copper particle for different diameters and water potentials, and thermal enhancements exceeding both theoretical and experimental results were observed in parallel with some other studies in the literature. The anomalous Green–Kubo thermal enhancement results could be explained by the interfacial dynamics and the neglect of calibrating the interaction potential to satisfy the physically-observed energy flow at the interface.


I. INTRODUCTION
Colloidal suspensions have recently been the focus of numerous scientific studies due to their potential engineering applications [1]. Studies of particulate systems can involve rheology, interface science and colloidal chemistry [2], and could lead to advances in foams, gels, paints, coatings and wetting-dispersing agents [3]. One recently-developed colloidal system is known as a nanofluid, and is a suspension of nanoparticles designed for enhanced heat transfer [4]. The physics of thermal transport in these systems requires further study, and simulations are frequently used for this purpose because of experimental limitations at the relevant length scale.

Molecular Dynamics (MD) is a powerful approach to investigate nanoscale dynamics
where interatomic forces govern the system behaviour. There are two different schemes that can be used to measure the heat transport properties in MD simulations: Non-equilibrium Molecular Dynamics (NEMD) [5] and Equilibrium Molecular Dynamics (EMD) [6]. A heat flux is imposed on the system in NEMD and the thermal conductivity is calculated based on the resulting temperature gradients. However, finite size effects are severe as a result of very high temperature gradients, and careful consideration is required to obtain reliable temperature profiles [7]. Equilibrium Molecular Dynamics instead calculates thermal conductivity from the time decay of heat flux fluctuations based on the fluctuation-dissipation theorem. While this requires more computational power, the method does not suffer from the drawbacks of NEMD and is widely used in the literature.
The thermal conductivities of Lennard-Jones liquids [8], water [9][10][11], methane hydrate [12], Ar-Cu nanofluids [13][14][15], and a water-platinum nanofluid [16] have been estimated with EMD, though most of these studies do not include a detailed error analysis. This is significant because the fluctuation-dissipation theorem is based on the chaotic movements of particles or atoms, which introduces statistical noise into the calculation of the autocorrelation function (ACF). Porter and Yip [17] and Chen et al. [18] proposed to truncate the integration time of the heat current autocorrelation function (HACF) to minimize the statistical noise. Other approaches include curve fitting [19], block averaging [20], random walk modelling [21], and spectral analysis of time series [22]. These studies provide fundamental insights about the statistical nature of the EMD method, and should be carefully considered to estimate the errors associated with Green-Kubo (GK) calculations.
In most nanofluid studies, Maxwell's relation [23] is used as a reference for the effective thermal conductivity of a simulated system. The Maxwell model neglects Brownian motion, nanolayering and agglomeration [24] as possible heat transfer enhancement mechanisms for nanofluids, but provides a lower limit by only considering the thermal conductivity of the mediums. Anomalous thermal enhancement values with respect to the Maxwell limit have been reported in some nanofluid studies that use the GK method, where the contribution of various mechanisms to this enhancement has been investigated by comparing simulations and theoretical calculations [15,16,[25][26][27]. However, inconsistent results and recent developments [18,19,22] concerning GK calculations in the literature indicate that there is no generally accepted GK algorithm, and that the effect of the GK parameters on the thermal transport has not yet been fully understood.
The anomalous thermal conductivity of nanofluids as found by GK calculations has been hypothesized to be the effect of Brownian motion [15] or nanolayering [26]. Recent studies have found an insignificant contribution of Brownian motion of the nanoparticles to the heat transfer enhancement of nanocolloids [28,29]. As for the nanolayering effect, both experimental and theoretical investigations have tried to quantify the thermal conductivity of the adsorbed liquid layer [30][31][32]. However, investigating nanolayering with MD simulations can be difficult since the interface potential is expected to have a considerable effect on the thermal transport at the solid-liquid surface. Most EMD studies use the Lorentz-Berthelot rules at the interface for lack of an alternative without adequately studying the effect of the surface parameters, even though the Lorentz-Berthelot rule has no physical justification.
It has been recently shown that Lorentz-Berthelot can lead to an overestimation of the interfacial thermal resistance for a hBN-water system [33], and an optimization process carried out by fitting potential parameters is required to accurately represent the interface.
This study presents a comprehensive evalution of the Green-Kubo approach by considering the thermal transport in a nanocolloid system, where several inconsistencies and anomalous thermal conductivity results have been reported in the literature. Using a standard force-field as the interatomic potential and a common modeling scheme for interface interactions, possible sources of the observed anomalities are revealed. A statistical assessment of several sources of error is presented and the effect of different water potentials on the thermal transport calculations is tested. The effect of the interface parameters on the thermal conductivity calculations is investigated. Finally, thermal resistances at the solid-3 liquid interfaces are quantified to further investigate interfacial effects. We believe that the presented results clarify the contributions of surface phenomena and the associated interfacial thermal resistance to the anomalous thermal enhancement found with Green-Kubo calculations.

A. Green-Kubo Relations
The Green-Kubo method relies on the fluctuation-dissipation theorem and linear-response theory [6], and describes the system behavior using time autocorrelation functions (ACF).
Integrating an ACF in time allows transport properties in the equilibrium state to be extracted [34]. The thermal conductivity is given in the GK formalism as: where k B is Boltzmann's constant, T is the absolute temperature, V is the system volume, J is the heat current vector, and the integrand is the heat autocorrelation function (HACF).
The HCAF is an ensemble average, which relates the thermal conductivity to the heat current in an equilibrated system. Equation (1) involves integrals in the continuous case, but since time is discrete in MD simulations, the integral is replaced by a summation in practice. The heat current vector is defined as: wherev i is the velocity of atom i,r ij andf ij are distance and force vectors between atoms i and j, α is the atomic species, and h α is the enthalpy of that species. The total energy e i is the sum of the kinetic and potential energies which can be expressed as: where Φ(r ij ) is the interatomic potential energy.
The enthalpy exclusion in the last term of Equation (3) for multi-phase systems was introduced by Babaei et al. [35] who studied a methane-Cu colloidal system using EMD, and was not common in earlier GK studies [15,16,27]. This term is subtracted from the heat flux because it represents the energy carried by the mediums but not transported between them. The species enthalpy is defined as: where N α is the number of atoms of species α.

B. Problem and Simulation Details
Pure water and water-Cu models with a single copper nanoparticle (diameter 1.3 or 1.8 nm) in the middle of cubic simulation domains with a side length of 2.8 nm and 4 nm were generated as shown in Fig. 1. Volume fraction calculations are not trivial at the nanoscale, considering that the clearence at the solid-liquid interface is not well-defined. We calculated the volume fraction to be 5% using mass fractions and the bulk densities of the phases, and 4.5% using the sum of volumes of the Voronoi polyhedra. This was kept fixed for all nanofluid systems, and is reported as 5% throughout the text. TIP3P (flexible) [36], TIP4P/2005 (rigid) [37] and SPC/E (rigid) [38] potentials were used to define the water interactions. Pure water systems of the same dimensions were used to estimate the thermal conductivity of the base fluid. The SHAKE algorithm [39] was used to enforce the bond and angle constraints in the rigid water models, and long-range Coulombic interactions were solved using particle-particle-particle-mesh (PPPM) summations [40]. Atomic interactions within a copper nanoparticle were defined using a Lennard-Jones 6-12 potential with Cu−Cu = 9.4353 kcal/mole and σ Cu−Cu = 0.2338 nm [41]. The Lorentz-Berthelot rule was used for interactions between oxygen atoms and copper atoms, and the timestep was set to 1 fs.
Periodic boundary conditions were applied in all directions.

B. Thermal Conductivity Calculations
The thermal conductivities of pure water and water-Cu nanocolloid systems are calculated following the procedure described in Section II A. The results for different water potentials, particle diameters and number of particles are presented in Table I The anomalous thermal enhancement observed with rigid water models exceeds both theoretical and experimental results, and has been attributed to an artificial particle correlation effect arising from the use of a single particle with periodic boundary conditions [25]. The suggested solution is to use multiple nanoparticles to break the symmetry of the system. This was tested for the rigid SPC/E water model using a three-nanoparticle system with the same volume fraction and particle diameter. Contrary to the hypothesis that artificial particle correlations are responsible for the anomalous enhancement, a further increase in the thermal conductivity was observed for the rigid SPC/E water model (Simulation 8).
Moreover, an insignificant difference was found between the one-and three-nanoparticles cases with the flexible TIP3P model (Simulations 4 and 7), suggesting that there may be some other effects responsible for the observed anomalous increase. The difference between water models persists and is even exacerbated for larger system sizes at a fixed volume fraction, with the larger single particle system having a slightly higher thermal conductivity for the TIP3P model (Simulations 4 and 9) and more than double the thermal conductivity for 9 the SPC/E model (Simulations 5 and 10). The abnormal thermal conductivity values associated with SPC/E models tend to increase with increasing particle surface area within the system as additional particles are introduced or particle size is increased (Simulations 8 and   10). Given the reasonable thermal conductivities of the pure water models, these findings suggest that the interfacial interactions that obey the Lorentz-Berthelot rules and that are different for each system are the cause of observed anomalous thermal conductivities.
This argument is supported by the results presented in Table II, where the estimated thermal conductivities for the three different water models are compared with the Maxwell limit [23] and experimental results of Xuan and Li [47] at the same volume fraction. The significant difference between our simulations and experiment regarding the particle size could be interpreted based on the findings of the experimental studies that investigate the effect of particle size to the nanofluid thermal conductivity. It has been showed that the thermal conductivity is increased with smaller size of nanoparticles due to the effect of enhanced Brownian motion [48,49]. Considering the larger particle diameter in the experiment, our effective thermal conductivity result with TIP3P model is quite reasonable.
However, switching the water potential to SPC/E or TIP4P results in anomalous thermal conductivity ratios exceeding both experiment and theory, implying that the parametrization of the system dynamics is responsible for the overestimation rather than any physical mechanism. When the results in Table I and Table II are considered together, our conclusion is that the surface interactions of solid-liquid interfaces in nanocolloids need to be accurately represented to avoid any unphysical interfacial thermal transport.

C. Interfacial Effects
The effect of surface interaction parameters on the effective thermal conductivity of solidliquid systems has been investigated before in several NEMD studies [33,50], but the effect of these parameters on GK calculations has not been studied to our knowledge. Our earlier NEMD study [33] showed up to a 12% change in thermal conductivity using an interface coupling tuned to satisfy the experimentally observed contact angle as compared to Lorentz-Berthelot. Surface interactions could have an even larger effect for EMD calculations since GK considers the energy flow associated with the motion of each atom at each time step by means of Equation (2), whereas NEMD uses an average temperature calculated by means TABLE I: Green-Kubo estimations of the thermal conductivity of water-Cu nanosuspensions with a single nanoparticle.k is the mean and σk is the standard error of the mean for different water models (Model), volume fractions (φ), number of particles (N p ), and particle diameter (D p ). Simulations 1-3 are pure water. Standard error of the means are calculated based on the long-time errors as mentioned above. σ r in the last column is the ratio of σk tok.   [47]. N p , D p and φ are the number of particles, particle diameter and volume fraction, respectively. k ef f is the thermal conductivity ratio of the nanosuspension, which is the ratio of nanofluid thermal conductivity to base fluid thermal conductivity. A Lennard-Jones 6-12 function is employed for the interfacial interactions: where is an energy scale and σ is a characteristic length. Ideally, and σ should be derived identify whether this could be the source of the anomalous thermal conductivity rather than some other simulation input. Existing investigations consider the effects of varying [50][51][52] and the same procedure is followed here. The effect on the thermal conductivity of varying at the water-Cu interface is reported in Table III, indicates that the interfacial energy parameter would easily be able to account for the overestimation of the thermal conductivity in Table I. Specifically, changing from the 1.21 kcal/mole given by the Lorentz-Berthelot rule for the SPC/E potential to the 0.98 kcal/mole for the TIP3P potential (Simulations 1 and 3) significantly decreased the thermal conductivity. Although these values of do not have any particular pyhsical basis, they do establish that the thermal conductivity is very sensitive to surface , and that using the values predicted by the Lorentz-Berthelot rule is likely a significant source of error. This reinforces the necessity of calibrating the interface potential based on experimental findings before using the Green-Kubo method to quantify thermal transport in nanocolloids.
Another observation is that the anomalous thermal conductivity values observed with both rigid models in Tables I and II were not observed with the flexible TIP3P model.
While this could appear as a possible reason for the abnormally high thermal conductivites, this conclusion is not supported by the reasonable thermal conductivities of the base fluids in Table I. Moreover, Table III provides direct evidence that an inaccurate interfacial potential could be responsible for the anomalous thermal conductivities.
The significance of the interface to the observed dependence of thermal transport on 12 FIG. 4: Two water blocks have been created in between three copper blocks to apply NEMD.
the water model could also be expressed in terms of interfacial thermal resistance. Based on prior work [53], this is measured using the NEMD method. Three different water-Cu systems were generated with the geometry shown in Fig. 4, and the Muller-Plathe algorithm [5] was applied to create a heat flux in the system in the z-direction. The temperature was calculated using a methodology presented elsewhere [53]. The interfacial thermal resistance was calculated as R = ∆T /q where R is the thermal resistance, ∆T is the temperature difference at the water-Cu interface, and q is the heat flux. Each system was equilibrated for 100 ps in the NPT ensemble, followed by a 1 ns production run in the NVE ensemble.
The thermal resistance results were obtained by averaging the results of 100 independent runs, and are presented in Table IV. A higher thermal resistance was found for the TIP3P model compared to the rigid models, and is consistent with our thermal conductivity results.
The discrepancy between the relative difference in the thermal resistances (around 60 %) and the thermal conductivities (more than 200 %) could be due to the fact that the thermal resistances were calculated using the NEMD Muller-Plathe algorithm, whereas the thermal conductivities were calculated using the EMD Green-Kubo algorithm. A second contributing factor could be the difference in geometry: whereas the nanocolloid contains highly curved surfaces, the thermal resistances are calculated for flat surfaces.

IV. CONCLUSIONS
The thermal conductivity of a water-Cu nanocolloid system was studied using the Green-Kubo method, with the intention of identifying the source of the anomalously high nanofluid sources of statistical errors, denoted as short-time errors and long-time errors. The magnitude of the error was estimated for both, and the long-time error associated with differences in velocity seeding was found to be larger than the short-time error associated with the fluctuations of a single autocorrelation function.
The anomalous thermal enhancement in the literature was reproduced using rigid SPC/E and TIP4P/2005 water models, but was not observed for the flexible TIP3P water model. Although the rigidity of the water model has appeared to be a possible reason for the observed anomalies, further simulations revealed that the discrepancy is related to the difference in the interfacial potential, pointing to the unsuitability of the Lorentz-Berthelot rules when calculating thermal conductivity with the Green-Kubo formulation. We conclude that the interfacial potential parameters should be carefully optimized to correctly simulate heat transport for solid-liquid systems when using the Green-Kubo method.