Interface-inspired formulation and molecular-level perspectives on heat conduction and energy storage of nanofluids

Aiming for the introduction of stability requirements in nanofluids processing, an interface-based three-step method is proposed in this work. It is theory-based design framework for nanofluids that aims for a minimum tension at the solid-liquid interface by adjusting the polar and dispersive components of the base fluid to meet those of disperse nanomaterial. The method was successfully tested in the preparation of aqueous nanofluids containing single-walled carbon nanotubes that resulted to be stable and to provide good thermal properties, i.e. thermal conductivity increases by 79.5% and isobaric specific heat by 8.6% for a 0.087 vol.% load of nanotubes at 70 °C. Besides, a system for these nanofluids was modelled. It was found to be thermodynamically consistent and computationally efficient, providing consistent response to changes in the state variable temperature in a classical Molecular Dynamics environment. From an analysis of the spatial components of the heat flux autocorrelation function, using the equilibrium approach, it was possible to elucidate that heat conduction through the host fluid is enhanced by phonon propagation along nanotubes longitudinal axes. From an analysis of the structural features described by radial distribution functions, it was concluded that additional heat storage arises from the hydrophobic effect.

where γ S p , γ S d , γ L p and γ L d are the polar and dispersive components of the solid and liquid surface tensions, respectively. Thus, in order to minimize γ SL , the magnitudes of γ S p and γ S d and also the magnitudes of γ L p and γ L d should be as close as possible. Therefore, the basic idea is to modify and to adjust with additives the components of the base liquid to meet those of the disperse solid because in such case both phases are expected to be compatible and so a stable colloid can be prepared. Given the case of carbon nanotubes (γ CNT  [25][26][27] . We prepared a series of aqueous standards with different concentrations (1.0·10 −3 , 2.0·10 −3 , 3.0·10 −3 , 4.0·10 −3 and 5.0·10 −3 vol.%) of Triton X-100 (Panreac © , surfactant for automatical analysis) and determined their polar and dispersive components. Although Triton X-100 surfactant has been previously used for the preparation of similar nanofluids [28][29][30] , its use has not been rationalised from this perspective. Young's equation provides a relation of them with the contact angle θ RL of a drop of these standards onto a surface (R) of reference: γ = γ − γ θ = γ + γ − γ + γ θ cos ( ) ( )cos A neat polished polytetrafluoroethylene (PTFE) stratum was chosen as surface of reference. Its polar and dispersive components were previously characterized by measuring the contact angles of a series of pure liquids whose polar and dispersive components are well known. At 25 °C and 1.0 atm in air atmosphere with 55% relative humidity, we found γ R p = 0. 8  so that the dispersive components of each standard can be directly acquired from its absolute surface tension γ L , which is determined using the drop-weight method, and its contact angle θ RL , using a goniometer and an optical subsystem 31 . The same stratum and ambient conditions were used. Table 1 summarises the results for each standard. From these values, we predicted the dependence of the solid-liquid interface tension of the SWCNT/TX/H 2 O system for increasing volumetric fractions of Triton X-100. This is the plot in Fig. 1. From the fitting function, we determined a 1.7·10 −3 vol.% of Triton X-100 in water that satisfies the minimum tension, which defines the host fluid of use. The last step, dispersion of the nanomaterial into the host fluid, comprises the immersion, disaggregation and mixing processes. Approximately 1 mg of as-received powder-like SWCNT was dispersed within 100 ml of host fluid by means of stirring and pulsed sonication (Sonics&Materials © , VCX-500, 13 mm probe tip) set at 50% amplitude with alternating on/off bursts of 2 s length (47 W per pulse). Aliquots were extracted after four, six and eight hours of sonication and centrifuged at 7500 rpm for one hour in order to remove aggregates of large size. Supernatants were collected for characterisation and labelled as NF4, NF6 and NF8 (a nomenclature that references sonication time in hours). They were the nanofluids analysed in this work.
Characterisation of nanofluids. Goodness criteria for nanofluid stability are constant load and size of heat carriers in suspension within the host fluid over a long-time period. The load was assessed in terms of www.nature.com/scientificreports www.nature.com/scientificreports/ spectral extinction using a UV-Vis extinction spectrometer (OceanOptics © , DH-2000-BAL light source and USB2000 + general purpose spectrometer), whereas the size was estimated by means of the dynamic light scattering technique (Malvern Instruments © , Zetasizer Nano ZS), both measured in triplicate every day during the first week and once a week during a full month. Also, the degree of dispersion was qualitatively studied with scanning electron microscopy (FEI © , Nova NanoSEM-450) images of SWCNT solid samples directly extracted from the host fluid by evaporation at 150 °C.
Since forced convection is the mechanism that rules macroscopic heat flow in a heat exchanger, quantifying the convection coefficient (h) is a minimum requirement to assess whether thermal performance is enhanced or not. The dimensionless Figure of Merit (FoM) h nf /h bf may be used as an indicator to quantify the relative enhancement of thermal performance which is solely due to changes in those properties of the nanofluid (nf) with respect to its base fluid (bf), which is the HTF in use. According to the Dittus-Boelter equation 32 (valid for fully developed turbulent flow in circular tubes under heating conditions), h nf /h bf can be defined as where ρ is the density,η is the dynamic viscosity, κ is the thermal conductivity and C p is the isobaric specific heat. Henceforth, these four properties were determined for each nanofluid. Density was measured using a pycnometer (Pobel © , 10 ml) at 25 °C and 1 atm. Dynamic viscosity was measured using a sine-wave viscometer (A&D Company © , SV-10) at 25 °C and 1 atm. Density and viscosity measurements were performed in triplicate. Thermal conductivity was measured with a transient hot-bridge (Linseis © , THB-100) operating at 25, 50 and 70 °C and 1 atm, with increasing input powers (10, 20 and 30 mW, respectively) for a better signal-to-noise ratio. Up to ten replicas were recorded for each sample and temperature. Lastly, isobaric specific heat was measured by means of temperature modulated differential scanning calorimetry (Netzsch © , DSC 214 Polyma), setting a temperature program that involves (i) a 5 °C/min ramp up to 40 °C, (ii) isothermal equilibration for 10 min, (iii) a 1 °C/min ramp down to 10 °C and (iv) a 1 °C/min ramp from 10 °C to 70 °C under temperature-modulated conditions with oscillations of ±1 °C amplitude and 1 min period.

Computational Framework
Modellisation of nanofluids. LAMMPS   www.nature.com/scientificreports www.nature.com/scientificreports/ in its all-anti conformation) were considered, with all their sites being explicit and a maximum occupancy of one atom per site. The nanofluid initial configuration was created by removing a cubic domain from an arrangement of equally spaced water molecules and then inserting a nanotube and two surfactant units. This ensures no overlap of atomic positions in the initial configuration and so that the simulation runs without instabilities. The OPLS-AA force field 35,36 was used to compute bonding and non-bonding atomic interactions within the model in liquid phase simulations. It is described by the following functional of classical potentials where K 1 is the harmonic force constant of each two-body interactions of covalently bonded atoms stretching from their equilibrium length 1 0 ; K θ is the harmonic force constant of each three-body interactions of covalently bonded atoms bending from their equilibrium angle θ 0 ; V m are the torsion barriers for the dihedral angle φ of each four-body interactions of covalently bonded atoms rotating around a bond; r ij is the distance between each pair non-bonded atoms, i and j, that interact as described by the energetics of the standard 12/6 Lennard-Jones potential, which is a function of the well depth ε ij and the collision diameter σ ij , and the Coulomb potential, which is a function of the atomic charges q i and q j . The parameters used for the OPLS-AA force field were adopted from literature [37][38][39][40][41] , as they were found to be compatible and transferable to our model. Lorentz-Berthelot combining rules were used to define the Lennard-Jones potential parameters for pair interactions between odd atoms, except for water-nanotube interactions, whose parameters were derived by Kaukonen et al. 38 from van der Waals density functional (vdW-DF) calculations with no empirical parameters.
Simulation of nanofluids. Equilibrium MD simulations were performed in a cubic environment with periodic boundary conditions imposed in all Cartesian directions. Velocities were rescaled with the Nosé-Hoover thermostat and barostat [42][43][44] (with damps of 10 and 100 fs, respectively) and updated using the Verlet integration algorithm 45 with a timestep of 1.0 fs. For the sake of computational time, Lennard-Jones pairwise interactions were limited to a cut-off distance of 10.0 Å, Ewald summation 46 was applied to compute long-range electrostatics beyond that cut-off and TIP3P water geometry was constrained using the SHAKE algorithm, so that dynamics run freely but bonds and angles are reset to their equilibrium length and angular values after each timestep 47 . Systems were initially forced to change their volume until density was equal to 1000 kg·m −3 (by doing so, velocities were rescaled for the atoms to match the box change and keep their natural motion), and then thermalized using the NVT ensemble at 25, 50 or 70 °C during 100 ps. Keeping up those conditions for 500 ps, transport coefficients were computed using the Green-Kubo formalisms 48,49 , considering a 10 ps correlation length and a 100 ps −1 sampling frequency. Autocorrelation functions were integrated using the trapezoidal rule numerical integration scheme. Lastly, the systems were pressurized using the isobaric-isothermal (NPT) ensemble at 1 atm during 100 ps. Now radial distribution functions (RDF) were computed, using the same correlation length and sampling frequency stated above. Simultaneously, density and isobaric specific heat values were computed every 10 ps up to 500 ps. Also, we dumped *.xyz files every 10 ps for the creation of 3D graphics in VMD 50 (NIH Center for Macromolecular Modelling & Bioinformatics, 1.9.3 version).

Results and Discussion
stability assessment. Changes in load and size of heat carriers were monitored as stability indicators. The load of SWCNT in nanofluids was related to their spectral extinction, since the intensity of light that is extinguished due to absorption and scattering phenomena as radiation passes through the liquid sample and interacts with them. Figure 2A shows the evolution of extinction values registered at 566 nm during a full month period after nanofluids were prepared. Wavelength choice is not arbitrary but meaningfully justified by S 22 transitions www.nature.com/scientificreports www.nature.com/scientificreports/ between van Hove singularities in [6,5]-SWCNT 51,52 , which is the most abundant chirality according to the information in supplier data sheet.
Attending to the initial values of extinction, the load of SWCNT introduced in the host fluid is greater for NF8, which is the sample treated under sonication for a larger time. However, there are no significant changes in the way they evolve. All nanofluids undergo a minimal aggregation process that is perceptible in a very early extinction decay that may be interpreted as a spontaneous change towards the equilibrium. Surface tension at the solid-liquid interface is low but still positive (see Fig. 1), so it can be expected from the system to reduce the interfacial surface. Indeed, we know from dynamic light scattering determinations, shown in Fig. 2B, that the average size (in terms of the hydrodynamic diameter) of heat carriers peaks a maximum value over the first two weeks in all cases. Size distribution width slightly increases over time, from about ±50 nm to ±100 nm. No phase segregation was observed by the end of the tests. So, the discussion above allows to conclude that NF4, NF6 and NF8 samples are stable once the equilibrium situation is reached, so their rheological and thermal properties are expected to be sustainable in time. What is more, the fulfilment of both stability criteria validates the theory-based design framework for nanofluids proposed in this three-step method for the preparation of stable carbon-based nanofluids in a polar host liquid. Figure 3A contains SEM images of as-received SWCNT commercial powder, revealing that carbon nanotubes are originally entangled and compact. Alternatively, an array of randomly oriented carbon nanotubes can be seen in Fig. 3B, which corresponds to a powder-like sample extracted from the nanofluid. The entanglement is significantly reduced after sonication, even though they are forming some bundles after the extraction process. More sharpness and higher magnification in Fig. 3C allow the identification of small groups of carbon nanotubes in solid state samples. By comparing all SEM images, it can be clearly appreciated that the sonication treatment included in the preparation protocol provides a good degree of dispersion and that the properly controlled addition of stabilizers favours the availability of individual SWCNT or bundles with a small number of them. thermal performance assessment. Density values for each nanofluid are available in Table 2. Given a standard volumetric flow in a heat exchanger system, that HTF with a higher density will provide a greater mass transfer coefficient and so a better heat transfer rate. This is numerically considered in the Dittus-Boelter equation introducing density values to the power of a positive exponent. However, we found density to increase up to 0.1% with respect to the base fluid, so it is no expected to be decisive for heat transfer enhancement. We estimated the volume fractions of solid (for discussions attending to the content of nanomaterial in each sample) from density by assuming that nf bf nt bf where ρ nf , ρ bf and ρ nt are the densities of the nanofluid, the base fluid and the nanotube.  www.nature.com/scientificreports www.nature.com/scientificreports/ Dynamic viscosity values are also represented in Table 2. It can be appreciated that they increase for increasing volume fractions of nanomaterial, up to 13.9%, which is still affordable for conventional pumping systems. We checked that the diluted amount of Triton X-100 does not cause a measurable difference in viscosity, so it is solely due to the presence of nanotubes. In Fig. 4 we compared our results with those predicted by those numerical models presented in Table 3. Einstein's model 53 for the viscosity of suspensions (Eq. 8) is based on the resolution of the hydrodynamic equations for suspensions of non-interacting rigid particles. Alternatively, Batchelor's model 54 (Eq. 9) takes into account the influence of two-particle interactions by introducing the Huggins coefficient k H , whose value was calculated by Berry et Russel 55 as 2/5 for interactions between rigid rods under steady shear flow conditions. In both cases the intrinsic viscosity {η} parameter is included, which is sensitive to the particle shape. Its theoretical expression, which can be found in literature 56 , is a function of the particle diameter-to-length ratio d np /L np . For the calculation of such a ratio, average diameter and length of SWCNT (available in the supplier data sheet) were considered. Despite this, none of these models was found to fit our experimental data, and the main limitation here may be the fact of nanotubes being considered as rigid bodies (as both models are constructed on this assumption) with a diameter-to-length ratio of 7.8·10 −4 whereas they can be stacked or partially folded in suspension, which should be properly represented by a larger aspect ratio. Still, the main conclusion that can be drawn from the comparison of these models in Fig. 4 is that the rheological behaviour of these nanofluids is 'better' described (since it is reasonably linear) by the Einstein's model and so that SWCNT, for low volume fractions in a host fluid with optimum surface tension, are more likely to behave as non-interacting (but not necessarily rigid) particles. This means minimum aggregation and maximum availability of heat carriers in suspension. Figure 5A shows thermal conductivity values of the base fluid, the host fluid and SWCNT/TX/H 2 O nanofluids, as a function of temperature. As it increases, that nanofluid with the highest volume fraction of SWCNT seems to be more sensitive to changes in temperature, peaking a 79.5% enhancement in thermal conductivity at 70 °C. These enhancements are better than expected by numerical predictions using models from Table 3, at any temperature. This is better appreciated in Fig. 5B-D. Hamilton-Crosser's model 57 (Eq. 11) is supposed to calculate the effective thermal conductivity of a biphasic mixture as a function of their individual thermal conductivities by considering the sphericity Ψ and the volume fraction ϕ of the disperse particles but, for this case, it estimates an enhancement below 2.0%. The applicability of this model is hardly limited by its poor approach to non-spheroid particle shapes. Aiming for a wider applicability of thermal conductivity predicting formulae, Yamada-Ota's model 58 (Eq. 12) was originally proposed for parallelepiped particle suspensions by introducing the length-to-diameter ratio L np /d np . For the calculation of such a ratio, once again, average length and diameter of SWCNT (available in the supplier data sheet) were considered. The prediction is closer indeed, but this model also fails to predict thermal conductivities of the nanofluids presented in this study. A discussion of transport properties of nanofluids on the basis of numerical models is clearly limited because they are systems of high

Reference Equation
Yamada-Ota 58 Table 3. Models to calculate the effective viscosity (Eqs 8 and 9) and thermal conductivity (Eqs 10 and 11), as a function of the volume fraction and shape of the disperse solid.
www.nature.com/scientificreports www.nature.com/scientificreports/ complexity that cannot be properly understood as common solid-liquid mixtures. For that we think the use of MD simulations helps for a better understanding of these systems and their properties. Now considering the low volume fractions involved in this work, these are, to our knowledge, the highest enhancements ever reported in literature for nanofluids containing carbon nanotubes. For instance, Choi et al. 59 reported a 159% enhancement of thermal conductivity in a synthetic poly-α-olefin oil with a 1.0 vol.% load of MWCNT at room temperature. Ding et al. 60  The values of isobaric specific heat are discussed. This property plays an important role in heat transport phenomena, as it determines the thermal energy storage capacity of the HTF in use. The specific heat of nanofluids is expected to decrease for increasing volume fractions, as the specific heat of solids is generally lower than that of liquids, but some experimental observations are opposed to this assumption 63 . Figure 6 illustrates that the isobaric specific heat of NF4 sample is lower than those of the base fluid and the host fluid at low temperatures, but  www.nature.com/scientificreports www.nature.com/scientificreports/ it significantly increases with temperature. NF6 and NF8 samples have higher specific heats over the entire range. A maximum 8.6% enhancement was found for the 0.087 vol.% nanofluids at 70 °C.
Contributions to heat transfer from these parameters must be assessed attending to a single criterion, what gives sense to the FoM defined in Eq. 5. Figure 7 gathers the values calculated for each temperature and volume fraction according to this formula. Any FoM value higher than 1.0 implies that nanofluids enhances the heat transfer process with respect to the base fluid. Best performance was found for the 0.087 vol.% SWCNT/TX/H 2 O nanofluid, which is more efficient than water by 31.3% under heating conditions at 70 °C, but still the criterion is satisfied in all cases, so any of these nanofluids could be a choice of preference over its base fluid to be used as HTF for thermal management. In overall terms, the conception of the three-step method leads us to the production of a whole series of SWCNT/TX/H 2 O nanofluids that were found to be stable and to provide an exceptional heat transfer efficiency for both conduction and convection processes according to the findings of an exhaustive systematic characterisation. So, the three-step method proposed here can be a promising advance in the field of nanofluids.
Understanding heat conduction. Before using the simulation approach for understanding and rationalising a system, it is important to validate the model with well-known experimental results. Figure 8A-D show computed values of the physical properties for both the base fluid and the nanofluid models, and how the presence of a carbon nanotube induces changes in the properties of the base fluid, which are clearly exalted. These values are slightly overestimated if compared with empirical evidences, and this is due to the limitation of the model size (the occupancy of carbon nanotube in the simulation box is equivalent to 1.4 vol.%, which is higher than experimental volume fractions by two orders of magnitude). Note that it is built as a finite box under periodic  www.nature.com/scientificreports www.nature.com/scientificreports/ boundary conditions that pretends to emulate an infinite bulk domain. The major concern here was the number of atoms within that finite box and for that a trade-off between system representativeness (i.e. the nanotube must be a minor component with respect to whole number of water molecules), calculation convergence and processing time was required for the simulation to be feasible. However, it is noteworthy that thermal conductivity (see Figs 5A and 8C) and isobaric specific heat (see Figs 6 and 8D) values follow experimental-like trends. This proves that the model offers good response to state variable changes and good representativeness of nanofluids, and also demonstrates the thermodynamic consistency of the chosen force field parameters for the simulation of its dynamics 64 . All these facts lead to a successful and unambiguous validation of the model for further analysis.
Defining the mechanisms that rule heat transport phenomena in nanofluids is a key factor towards the comprehension of their thermal properties. This is a controversial issue among researchers. Keblinski et al. 65 , on the basis of MD simulations using a hypothetical nanofluid, suggested that the major contribution to thermal conductivity is the ballistic heat transport through in solid nanoparticles promoted by phonons and also that solid-liquid interfaces may influence the overall process in nanofluids. Evans et al. 66 , following a similar scheme, stated that thermal conductivity enhancements are not associated to hydrodynamic effects induced by Brownian motion. Alternatively, Sarkar et al. 67 , using copper-argon system model, concluded that is not the Brownian motion of nanoparticles but the enhanced movement of the surrounding base fluid atoms what promotes localized heat transport. These proposals, although they are far from being conciliated in a heat transfer theory for nanofluids, share a common point: nanoparticles modify the natural dynamics of the base fluid.
Our thermal conductivity results encouraged the need for understanding the underlying heat conduction mechanisms in nanofluids. In this work, we proceeded with the equilibrium approach for the analysis of heat conduction through nanofluids using a model whose dynamics were observed in a MD environment. For such approach, the Green-Kubo formalism 48,49 provides the following time-correlation function to calculate thermal conductivity where V is the system volume, T is the system temperature, k B is the Boltzmann constant and 〈J(t) · J(0)〉 is the heat flux autocorrelation function (HFAF). This function describes how long it takes for the fluctuations (about zero at equilibrium) of the heat flux vector J(t) to dissipate, which implicitly contains information about the heat conduction process through the nanofluid 68 . Figure 9A contains the HFAF of the base fluid and the nanofluid at the same temperature. Decay pattern to zero was found to be monotone for the base fluid but oscillating for the nanofluid, in agreement with Keblinski et al. and Sarkar et al. 65,67 Correlation lasted longer due to these thermal fluctuations, thus leading to a higher thermal conductivity of the nanofluid upon integration. For increasing temperatures, as it is shown in Fig. 9B, the decay ratio is lower and that explains why thermal conductivity of these nanofluids is sensitive to changes of this state variable. Figure 9C shows the spatial components of the HFAF (this is J x (t) · J x (0), J y (t) · J y (0) and J z (t) · J z (0)) for the nanofluid system at the highest temperature. Thermal fluctuations within the x-and z-directions, which are orthogonal to the nanotube longitudinal axis, depict a random sequence but those within the y-direction, which contains the nanotube longitudinal axis, depict some periodic features that resemble much of 'heat beats' , as it is shown in Fig. 9C. We verified that the frequency of those features in the HFAF fits with that expected for the superposition of resonant n = 1 and n = 2 vibration modes of phonons in a well whose boundary width is equal to the nanotube length 69 , together with a higher frequency pattern that matches the oscillations in the HFAF of the base fluid, but with a much larger amplitude. This could be due to phonon-enhanced water flow, as we found water molecules occupy the nanotube inner cavity, as is discussed below. For better appreciation, we found convenient to plot the superposition of the temporal parts of the classical waves with the matching frequencies.
Given the statements above, we concluded that presence of nanotubes disturbs the natural dynamics of the base fluid by increasing the phonon mean free path along the nanotube longitudinal axis (local anisotropy), thus promoting the correlation of thermal fluctuations and leading to a higher thermal conductivity. If extended to a bulk system, we suggest that coherent and consecutive phonon propagation events between randomly oriented nanotubes (average isotropy), together with phonon-enhanced water flow, are responsible of the macroscopic thermal conductivity of the nanofluid. It is important to note that it is the Green-Kubo equilibrium methodology what provides access to this information (which is otherwise missed by non-equilibrium methodologies that impose a temperature gradient in a certain direction) and so we assert that it meets the requirement for the analysis of heat conduction at a non-continuum scale in nanofluids.
Understanding heat storage. The enhancement of isobaric specific heat in nanofluids has received much less attention than thermal conductivity from theoretical research. Shin et al. 70 , in a consideration on the mechanisms responsible of this anomalous enhancement, proposed three independent modes for heat storage in nanofluids: (i) higher specific heat of nanoparticles compared to the bulk solid, due to the high surface energy per unit mass; (ii) additional storage arising from solid-liquid interactions at interfaces; (iii) layering of liquid molecules around nanoparticles, configuring a solid-like layer with higher specific heat than the bulk liquid. Jo et al. 71 , on the basis of MD simulations, were able to estimate the thickness and the density of a compressed liquid layer on the surface of graphite nanoparticles in a binary carbonate molten-salt. Navas et al. 13 , by analysing radial and spatial distribution functions from MD simulations, observed the arrangement of base fluid (a eutectic mixture of biphenyl and diphenyl oxide) molecules in a rigid layer surrounding copper and nickel nanoparticles and also concluded that the nature of interactions determines how tightly the layer is bound to the surface and how that boost heat transfer through the interface.
Since our model was found to reproduce that anomalous heat storage enhancement, we proceeded with an analysis of the equilibrium molecular architecture from MD simulations. Computed RDFs for the C nt ··· O bf pair, which are shown in Fig. 10A, revealed an anomalous population of water molecules within a virtual shell defined between 3.0-4.0 Å from the nanotube surface. This accounts for both water inside and outside the nanotube (see Fig. 11A for better appreciation). C nt ··· O bf RDFs were again computed for a system that was forced not to have water molecules inside the nanotube and plotted in Fig. 10B. Number density of water molecules outside the nanotube is not consistent with sharp layering due to specific interactions but with some short-range order induced by the presence of the nanotube, that is distinguishable of bulk disorder and not altered by temperature  www.nature.com/scientificreports www.nature.com/scientificreports/ (although loss of order would be expected as kinetic energy increases, it seems not to be significant for the studied range of temperatures). This observation allows a molecular interpretation of the isobaric specific heat increase in terms of the 'hydrophobic effect' 72,73 : since nanotubes are non-polar, its inclusion in water forces the formation of a cage of molecules around them (i.e. some sort of clathrate), with those molecules being in a more 'solid-like' (yet dynamic) situation due to constrains arising from a higher number atom-atom interactions per molecule. Compared to liquid water, more heat is required to break up that structure, hence isobaric specific heat is increased.
On the other hand, RDF revealed no significant number density for any interaction of Triton X-100 atoms within a 10.0 Å cut-off radius from nanotube carbon atoms and, given its spatial configuration (see Fig. 11B,C), it seems not to be fully adsorbed onto the surface. The polyether chain is indeed preferentially located in the aqueous phase; as nanotube carbon atoms exhibit no natural charge separation, the extension of its adsorption is limited to van der Waals interactions to the nanotube surface, which are weaker than electrostatics interactions making water more competitive for solvation. The non-polar head is more akin to the surface but still attached to the polyether chain, which makes it more steric demanding and unable to approach towards the surface. The cumulative result of the above-stated forces suggests that Triton X-100 is partially adsorbed, yet located in the vicinity of the solid-liquid interface and altering its energetics. This fact may serve as a piece of evidence to conclude that the stabilisation provided to SWCNT/TX/H 2 O nanofluids is not arising from electrosteric events but from its influence on those energetics, changing cohesive forces and thus change surface tension.

Conclusions
Next stage on the roadmap of nanofluids is the improvement of colloidal stability for sustainable thermal properties. An interface-based three-step method, arising from the rationalisation of the minimum tension at the interface between the disperse solid and the host fluid, was proposed in this paper. It was tested for the preparation of aqueous nanofluids containing single-walled carbon nanotubes, using the non-ionic surfactant Triton X-100 for the purpose of adjusting the polar and dispersive components of the base liquid until convergence with those of the disperse solid, thus making both phases compatible. Characterisation procedures revealed good dispersion and stability, in terms of constant particle load and size distribution, during at least one month. No significant changes in rheological properties were found but an extraordinary 79.5% enhancement in thermal conductivity and 8.6% in isobaric specific heat were observed for that nanofluid with a 0.087 vol.% load of nanotubes in suspension at 70 °C. Also, the heat convection process, as assessed by the Dittus-Boelter Figure of Merit, was enhanced by 31.3% for same nanofluid at the same temperature. Overall stability and thermal performance endorse the proposed three-step method and the method itself fulfils the target of sustainable thermal properties.
These results promoted the interest for the understanding of structural and dynamic features explaining thermal properties of nanofluids. For that, a homologous nanofluid system was modelled and validated for its simulation in a classical Molecular Dynamics environment. The spatial components of the heat flux autocorrelation function of the system revealed that heat conduction is significantly enhanced along the nanotube longitudinal axis due to a larger phonon mean free path. We suggested that consecutive phonon propagation events between randomly oriented nanotubes are responsible of the macroscopic thermal conductivity observed for these type of nanofluids. It was possible to elucidate from an analysis of the radial distribution functions that liquid layering phenomena occurs and may be related to isobaric specific heat enhancement due to the hydrophobic effect. Also, that Triton X-100 surfactant is partially adsorbed on the solid-liquid interface and so the stabilisation mechanism is not arising from electrosteric hindrance but from surface tension minimisation. Here modelling and simulation have supported experimental characterisation by providing information that cannot be obtained empirically.