Specific heat capacity enhancement studied in silica doped potassium nitrate via molecular dynamics simulation

Abstract

Molten salts serve an important purpose for short term heat energy storage and as heat transfer fluids in solar power plants. Different experimental groups have shown that certain mixtures containing salts doped with small amounts of nanoparticles exhibit much greater specific heat capacities compared to the same base salts without nanoparticles. This effect is technically interesting and economically important. Thus far, however, it is not understood. Our aim is the theoretical investigation of the specific heat capacity in the aforementioned nanofluids on the molecular level using simulations. Here we present results for liquid potassium nitrate doped with silica nanoparticles. We discuss the observed increase of the specific heat in terms of the particle induced hydrodynamic reinforcement and liquid structure. The theoretical background of this discussion is a ω-space resolved phonon theory of liquids in conjunction with differential spectral densities, computed for the different systems with and without nanoparticles.

**Subject ** Atomistic models, Nanoparticles

Introduction

Nanofluids consist of solvents containing suspended nanoparticles. They latter may cause strongly altered thermophysical properties in the base liquids. These changes are sometimes quite remarkable in the sense that the effects are much greater than expected on the basis of models, which describe the properties of interests in terms of a weighted sum of the corresponding properties of the individual components^1^. Typically the nanoparticle content varies between 0.1% to 5% by weight. Among these thermophysical properties are the thermal conductivity, viscosity, or the specific heat capacity^2^. This makes nanofluids particularly interesting in heat transfer applications - especially in the context of solar thermal technologies^3,4^.

In this study we are interested in the specific heat capacity of certain types of heat transfer fluids. This thermophysical quantity, measured in units of J/gK, should be as high as possible^5^. In a previous paper^6^ one of us has given a comprehensive overview of the attendant literature. If small amounts, typically around 1% wt., of different nanoparticle, e.g., alumina^7–9^, silica and copper oxide^8^ are added to water the result is negative, which means that the studies show no increase of the specific heat capacity. In the case of molten salts as base fluids the results are quite different and significant increases of the isobaric specific heat capacity, cP, between a several percent to more than 20% are obtained - again for nanoparticle concentrations of around or even less than 1% wt.

Most studies focus on eutectic mixtures because of their practical significance, e.g., Li2CO3/K2CO3 doped with carbon nonotubes^10,11^, graphite nanoparticles^12^, SiO2^13–15^, or alumina^16^. These authors have reported specific heat capacities about 120% above the base in eutectic (62:38) mixtures of Li2CO3 with KCO3 containing 1.5% wt. silica nanoparticles in the range from 2 to 20 nm^17^. Another frequently studied base salt mixture is the NaNO3/KNO3 (60:40) eutectic. Again, small concentrations of alumina^18,19^, silicon dioxide^19–21^, copper and titanium^19,22,23^ as well as mixtures of silicon dioxide and alumina^19^ yield pronounced increases of cP. Some of these papers report that the size of the nanoparticles has a significant effect, e.g.,^17,20,21^ (Li2CO3/K2CO3 eutectic) or^18^ (NaNO3/KNO3 (60:40) eutectic). Another observation in some of the studies is that the heat capacity exhibits a maximum. Chieruzzi^19^ and Lasfargues et al.^22,23^ observe this maximum in the range between 0.1 to 1% wt. nanoparticle. A similar maximum is observed in ref.^24^, where the authors use the same base salt doped with Al2O3 nanoparticles. Yet another confirmation of the existence of a maximum heat capacity is found in ref.^25^. The authors study the eutectic mixture of NaNO3 with KNO3 containing silica nanoparticles in the range of concentrations from 0.5 to 2% wt. Their maximum enhancement of the heat capacity is 25%. It is important to mention that the observation of an increased heat capacity is not limited to nanofluids based on binary salt mixtures. In ref.^26^ the authors study the single component base fluid KNO3 containing silica and alumina nanoparticles of variable size, and observe an increase of cP in the solid as well as in the liquid phase. Addition of 1% wt. nanoparticles increases cP in the solid phase by between 5 to 10%. The corresponding increase in the liquid phase is roughly 6%. In ref.^27^ the authors investigate a molten ternary salt mixture. Focusing on the optimal concentration of the alumina nanoparticles, they obtain a specific heat capacity increase of about 20% at only 0.063% wt. nanoparticle concentration. An even larger enhancement of cP was found in ref.^28^, where the author studies various ionic liquids containing nanoparticles. In ionic liquids the cations are large organic species and anions are organic or inorganic species. It is worth noting that the marked increases of cP reported in this reference appear to be monotonic in the nanoparticle volume fraction over the entire range of volume fractions ranging up to 2.5% wt. Aside from water, molten salts or ionic liquids, other systems were studied as well. Some show no increase or even exhibit a decrease of cP. Examples can be found in refs^29–31^. Here the authors do find reduced specific heat capacities of nanofluids containing zinc oxide, silicon dioxide, and alumina nanoparticles, respectively, dispersed in water/ethylene glycol mixtures compared to the base fluid. Another example is ref.^32^, reporting studies on silica-ethylene glycol, silica-glycerol, and silica-glycerol/ethylene glycol (60:40) mixtures (by mass), where the nanoparticle concentration ranges from 1.0 to 4.0% wt. Still other nanofluids show a considerable increase of cP, for instance Zhou et al.^33^. They study the specific heat capacity of ethylene glycol-based CuO nanofluids and find a maximum increase of roughly 6%. Nelson et al.^34^ report a cP of polyalphaolefin, which is increased by 50% due to the presence of graphite nanoparticles, when their mass fraction is around 0.6%.

We can summarize these results with the observation that the specific heat capacity of salt mixtures is independent of the type of nanoparticle. However, the mass concentration of particles should not be greater than 1% wt. Nanofluid preparation may vary, e.g., the temperatures during drying of the samples. But despite of this, cP-enhancements are between 10 to 30% on average. What is lacking thus far, however, is a molecular theory of the specific heat capacity of the above heat transfer fluids and its dependence on nanoparticle content.

Here we present the results of molecular dynamics (MD) simulations of KNO3 containing silica particles. We choose this system, because it exhibits an increase of cP as has been shown experimentally^26,35^. A one-component base liquid is preferable, because it is simpler to parameterize as compared to a two or even multi-component base liquid. This eliminates one source of possible error. In addition, the system has already been examined using MD before^35^. We present results for the specific isochoric heat capacity, cV, the specific isobaric heat capacity, cP, as well as the isobaric expansion, αP, in the concentration range from 0 to 40% wt. Our suspended nanoparticle possesses a diameter of about 1.8 nm. The simulations produce a clearly discernible local maximum or ‘hump’ of cV as well as of cP at around 2% wt. The increase of cP at the local maximum is around 3%, which is about half the increase observed in ref.^26^. In addition we observe an approximately linear increase of cV, which results in a 10% increase at the high end of the concentration range compared to a linear interpolation between the respective pure components, i.e., neat KNO3 and isolated silica particles. This suggests a second effect, aside from the one causing the ‘hump’ at low concentrations, possibly due to additional vibrational modes inside the highly structured base liquid separating the nanoparticles. Because αP decreases with increasing nanoparticle concentration, we find that cP, aside from the ‘hump’ at low concentrations, does not exhibit such a pronounced linear increase as exhibited by cV with rising particle concentration. We discuss the observed increase of the specific heats in terms of the particle induced liquid structure and the differential spectral densities, which we compute for the different systems with and without nanoparticles. The theoretical background of this discussion is the phonon theory of liquids put forward in ref.^36^.

The paper is organized as follows. The next section contains the technical details of the simulation methodology and systems studied here. Subsequently we present the simulation results including their discussion. A concise overview of theoretical concepts serving as basis for this discussion is compiled in an appendix. The final section is the conclusion.

Results

Specific heat capacity of the pure components

The upper panel in Fig. 1 shows the isochoric heat capacities of isolated silica particles. The three different data sets were obtained as follows. The classical result, cV,cl, follows via

where m is the mass of the particle, kB is Boltzmann’s constant, and 〈δE^2^〉 is the classical mean square energy fluctuation in the canonical ensemble realized using a Nosé-Hoover thermostat^37^. As a solution for the undisturbed quantum mechanical harmonic oscillator, cV,qm is computed via

where u=ℏω/kBT. The spectral distribution g(ω), obtained also from the molecular dynamics simulation, is explained in detail below. Finally, the quantum corrected classical specific heat is given by cV,cl + ΔcV, where

Figure 1 (a) Specific isochoric heat capacity, cV, of SiO2-particles at T = 700 K vs. particle radius, r. The arrow indicates the experimental bulk value discussed in the text, i.e., 1.12 J/gK. (b) cV vs. T for four particles sizes, i.e., r = 5, 9, 13, and 17 Å. The size of the symbols correspond to the size of the particles.

Notice that the first term in brackets is the contribution to the isochoric heat capacity of a harmonic oscillator with the frequency ω as in Eq. (2), whereas the last term in the brackets is the classical limit of this expression. Another approach utilizes cV=(kB/m)∑i(βℏωi)2eβℏωi(eβℏωi−1)2 (β = (kBT)^−1^). The summation is over all normal mode frequencies calculated via a standard normal mode analysis of the particles. The resulting data are not included, as the findings are close to cV,qm computed via Eq. (2).

Note that the reduction of cV with increasing particle radius, r is in accord with a previous theoretical work^38^ as well as with the experimental data. A simple argument for the observed decrease of cV with increasing r is the following. The frequency of a mode, ω, behaves according to ω^2^ = k/mk, where k is an average ‘spring’ constant (or stiffness) and mk is a typical mass associated with this mode. The latter is limited by the linear dimension r of the nanoparticles. A large particle can support longer wavelength involving more mass, thus mk ~ r. The average k is smaller for a small particle, because of its comparatively larger surface (where the coupling is weaker), thus k ~ r^2^ (at least approximately) and therefore ω^2^ ~ r. Thus, we expect more ‘soft’ and therefore more highly excited modes for smaller particles - mainly due to the weaker coupling at near the surface (with the exception of the silanol groups of course).

The bulk value of cV for silica at T = 700 K is close to 1.12 J/gK, which we find via a quadratic extrapolation based on the experimental values provided in ref.^39^ at 400, 500 and 600 K. The value is in good accord with our results for cV,qm and cV,cl + ΔcV (see also ref.^40^). The classical result is higher, as expected. Note that the Debye temperature for silicon dioxide is slightly less than 500 K in both the amorphous and crystalline states^41,42^, which is still not too far below our temperature, i.e., at 700 K we have not yet reached the Dulong-Petit limit. The bottom panel of Fig. 1 shows cV vs. T for four particles sizes, i.e., r = 5, 9, 13, and 17Å. The temperature range is the range studied in this work. Notice that the upper set of classical heat capacities are close to horizontal, as they should. Their quantum counterparts, however, increase with increasing temperature, due to the increasing excitation of higher frequency modes.

The following figure, Fig. 2, depicts the isochoric heat capacity of neat KNO3 vs. temperature in the liquid phase. Let us look to the experiments first and briefly estimate cP for neat KNO3 via its Dulong-Petit (or high temperature) limit, i.e., cV = RNf/(2mmol). R is the gas constant, Nf is the number of degrees of freedom, which here is close to 6nA, where nA is the number of atoms per mole substance and mmol is the molar mass. In the case of KNO3 the number of atoms is nA = 5 and the molar mass mmol = 101 g mol^−1^, which yields 3RnA/mmol = 1.23 J/gK. The corresponding experimental value is 1.14 J/gK. We calculate this number via cV=cP−TαP2/(ρliqκT) from cPexp above the melting point, i.e., cPexp≈1.39J/gK in the temperature range from 620 K to 730 K^43^ (in the (older) literature cPexp varies roughly between 1.35 and 1.41 J/gK (for instance^44–46^). Chieruzzi et al.^26^, who also study the effect of different types of nanoparticles, obtain the lowest value, i.e., cP = 1.18 J/gK.). Using the values for the isobaric thermal expansion αP and the liquid density ρliq from ref.^47^, i.e., 3.89 < 10^4^αP < 4.07 K^−1^ and 1.778 < ρliq < 1.857 g/cm^3^ for 620 < T < 730 K, and the isothermal compressibility from ref.^48^, i.e., 2.1 < 10^10^κT < 2.7 Pa^−1^ for 620 < T < 730 K, we find 0.24<TαP2/(ρliqκT)<0.25 J/gK for 620 < T < 730 K. We note that the density of KNO3 at 700 K in the simulation is found to be 1.73 ± 0.01^3^ g/cm^3^, which is less than 2% below the corresponding experimental value. In Fig. 2 we find an overall reasonable agreement with the experiment if we include the quantum corrections ΔcV computed according to Eq. (3). Notice that cV is almost constant, which again agrees with the experimental evidence. In the following discussion we shall frequently refer to different concepts of heat capacity in liquids, all of which are summarized in the appendix.

Figure 2 Isochoric heat capacity of neat KNO3 vs. temperature in the liquid phase. The horizontal line indicates cV based on experimental data, which shows no clear evidence for a significant temperature dependence in this temperature range.

Firstly, it is interesting to note that generally a decrease of cV with increasing temperature is observed (Here we observe a slight decrease in the uncorrected data only.). In Eq. (11) for cV, in the appendix, the temperature enters mainly through the two factors (1 − r^3^/3), where r is the ratio of the Frenkel frequency, ωF, to the Debye frequency, and (1 + αPT). The latter generally causes an increase with increasing temperature. The former, however, decreases with increasing temperature. This is because r∝ωF and the Frenkel frequency, which really is due to a gap in momentum space rather than being a frequency gap as a recent simulation study has shown^49^, is expected to increase with increasing temperature. This in turn can give rise to the usually observed decrease of cV with increasing temperature T in liquids. It is this explanation of the decrease of cV, that lends particular support to the phonon theory of liquids^50^ (cf. the extensive attendant discussion in ref.^36^). In the present case, as already mentioned, cV appears to increase slightly with increasing T (cf. Fig. 2). But this increase is much less than what is expected based on the factor (1 + αPT).

Specific heat capacity of the nanofluid

The next figure, Fig. 3, presents a key result, i.e., the isochoric heat capacity, cV, of KNO3 doped with SiO2 vs. wnp, the weight fraction nanoparticles, at T = 700 K. The upper panel shows an enlarged view of the low end of the wnp-axis. The radius of the silica particle used here is 0.9 nm. There are two sets of data points in each panel. The upper data set is the direct result of the MD simulations, which yields the classical cV via the equilibrium fluctuations of the energy in the canonical ensemble. The straight line is a linear interpolation formula connecting cV,liq with cV at around wnp = 0.4. Note that the data in the range 0.05≤wnp<0.1 fall systematically above this line. Thus, our simulations do show a small but clearly discernible increase of cV beyond the linear interpolation connecting the limits of low and high wnp. This is highlighted by the shading bounded by the linear interpolation from below and the solid line from above. The latter is a fit through the simulated cV-values based on the interacting mesolayer model, i.e., Eq. (14) in the appendix, developed in ref.^6^. The numerical values of the adjustable parameters, i.e., Δ/R und κmax, are included in the figure (here: κmin = 1). The quantum corrections mainly cause a constant vertical scaling of the data points. Notice again that the resulting cV of neat KNO3, including the correction, is close to its experimental value discussed above, i.e., 1.14 J/gK. The maximum of the apparent local increase of cV in the doped system, beyond the linear interpolation, is close to 3%. However, the wnp-range over which we observe this increase is larger than expected from previous experiments, which usually show no discernible effect when wnp is close to or greater than 0.05 (cf. ref.^6^).

Figure 3 (a) and (b) Isochoric heat capacity, cV, of KNO3 doped with SiO2 vs. wnp, the weight fraction nanoparticles. (c) Isobaric heat capacity, cP, vs. wnp. The lines are discussed in the text.

It is important to note that the limiting cV-value for wnp → 1, including quantum corrections, according to the linear interpolation formula is significantly larger, i.e., 1.39 ± 0.03 J/gK, than the value expected according to Fig. 1 showing the results for isolated nanoparticles. Because the limit wnp → 1 is not realistic for spherical particles, it is better to make this comparison near the high end of the wnp-range considered here, i.e., wnp = 0.4. But even there we find that cV is 10% higher than what we expect on the basis of a linear interpolation between the respective pure components, i.e., neat KNO3 and isolated silica particles. This suggests a second effect, aside from the one causing the ‘hump’ in the range below wnp ≈ 0.1. As will be discussed below in more detail, we attribute the increase of the specific heat capacity at high wnp to additional vibrational modes made possible by the pronounced density modulations induced in the liquid salt between nanoparticles. However, a simple example may serve to illustrate this point. Consider a molecule adsorbing on a surface from the gas phase. This molecule gains additional heat capacity due to the ‘conversion’ of translational motion into oscillatory motion within the surface potential well.

The bottom panel in Fig. 3 shows our estimate of cP vs. wnp computed via the data for cV in the middle panel in conjunction with the results for the isothermal expansion coefficient, αP and TαP2/(κTρ) both shown in Fig. 4. Here αP is computed via two distinct routes. The first one is a heating-cooling-schedule, which yields the volume at fixed wnp at three different temperatures including and bracketing 700 K Subsequently, αP is computed on the basis of a polynomial fit through these points (the average is over 4 distinct systems). The second approach utilizes the fluctuation formula αP = kBβ^2^〈δVδH〉/V. The difference is largest for wnp = 0.01, but otherwise the agreement between the two methods is quite good. In order to obtain TαP2/(κTρ) we compute the isothermal compressibility via κT = β〈δV^2^〉/V. Subsequently we use a simple linear fit to these data to calculate cP. We note that the quantum corrected cP in Fig. 3 is in good accord with the experimental result in the limit wnp = 0. As in the case of cV, we do observe the ‘hump’-like feature, exceeding the linear interpolation by 3% at its maximum. Figure 5 shows the radial pair correlation function, g2(r), where here r is the distance from the center of the particle, for the ions at different wnp. We do observe substantial ordering imposed by the nanoparticle on the ions of the liquid. Of course, the periodicity of the system will enhance this effect, but so will the presence of neighboring particles in a true experimental situation^51^. The real difference to most experiments is the small size of the particles in our simulation, which is more than an order of magnitude less than common experimental particle sizes. In other words, in a real experiment, when the particles are well dispersed, we can expect that their separation at the same wnp is greater than in the simulation. Thus, on average, the surface induced ordering may be less than what is observed here (however, experimental particle sizes are average values and the widths of the attendant distributions are large, i.e., a sizable fraction of the particles may indeed be much smaller).

Figure 4 (a) Isobaric thermal expansion, αP, vs. wnp, the weight fraction nanoparticles. (b) Difference between specific isobaric and isochoric heat capacity, cP − cV, vs. wnp. The straight line is a fit to the data used to compute cP in the previous figure from the simulated cV-values.

Figure 5 Radial distribution function g2(r) of silica particles with radius 9 Å obtained for the indicated weight fractions (r in this figure should not be confused with the particle radius). Note that the origin is the center of mass of the particle. The shading indicates the extend of the nanoparticle.

A useful quantity, which allows to study the influence of the nanoparticles on vibrational modes, is the spectral density g(ω), which already appeared in Eqs (2) and (3). Following ref.^52^ we compute g(ω) via

where

Note that ℱ〈x˙j(to)x˙j(to+t)〉 is the Fourier transform of the autocorrelation function of the atomic velocity components. In particular, for ω = 0, this will be the diffusion coefficient. Figure 6 shows the spectral density vs. frequency, ν = ω/(2π), in neat KNO3 at 700 K. Notice the rather broad peak at small frequencies, corresponding to longer wavelength collective modes, in contrast to the rather narrow peaks at higher frequencies due to molecular modes. The effect of the presence of the nanoparticles on the spectral distribution g(ω) is best seen if we study the difference spectra Δg(ω)=gdoped(ω)−gbasefluid(ω), which are shown enlarged in Fig. 7. The figure is composed of four panels, showing Δg(ω) vs. ν in enlarged sections on the frequency axis of the previous figure.

Figure 6 Spectral density, g(ω), vs. frequency, ν = ω/(2π), according to ref.^52^ obtained for neat KNO3 at 700 K. Also shown are difference spectra for doped KNO3. The same spectra are shown enlarged in the next figure. Note that W^Δ^(ν), which is the term in brackets in the argument of the integral on the right side of Eq. (3), is the statistical weight of the classical contribution to cV at that frequency.

Figure 7 Difference spectral density, Δg(ω), vs. frequency, ν = ω/(2π), following ref.^52^ obtained for doped KNO3. Here the spectrum for neat KNO3 obtained under otherwise identical conditions is subtracted from the the velocity spectra of the doped systems.

We first look at the upper left panel in Fig. 7. Notice that the difference spectral density in the limit of small frequency is negative for the smallest concentration of nanoparticles, i.e., wnp = 0.01. For the larger concentrations, i.e., wnp = 0.05 and 0.1 it is positive. This indicates an excess of low frequency, i.e., classically excited, modes in the doped fluid at sufficiently high concentrations. The remaining panels show the difference spectra at the positions of the narrow peaks in Fig. 6. Here we find that the difference is overwhelmingly positive for wnp = 0.01. Overall this is much less the case for the higher particle concentrations.

As pointed out before, in the simulation we observe two apparent types of cV-enhancement. One is the ‘hump’ at low wnp. The other one is an overall near linear increase with wnp. Both types of of cV-enhancement are likely due to additional vibrational modes in excess to those present in the base liquid or the nanoparticles, individually. Additional modes can be generated by the ‘stiffening’ or hydrodynamic ‘reinforced’ of the liquid due to the nanoparticles, i.e., by a decrease of the Frenkel frequency, ωF. Additional modes may also occur due to the surface induced structuring of the liquid between particles, as mentioned above. Unfortunately, if we compare the cV-enhancement for the three particle concentrations, we notice that they differ very little. Note that we compare cV,cl + ΔcV at wnp = 0.01, 0.05 and 0.1 in Fig. 3 (middle panel) relative to an almost horizontal line intersecting with the vertical axis at the cV of the pure salt. Figure 7 tells us that in the case of wnp = 0.01 the enhancement most likely is due to the frequencies around 1.8/ps and 3.25/ps. In the case of wnp = 0.05 and 0.1 a significant contribution is due to very low frequencies. However, notice also the following point. The experimental velocity of sound, cs, in KNO3 at about 650 K is approximately 1700 ms^−1 ^^53^. In addition, the simulation box dimension, L, for wnp = 0.05 it is about 3.8 nm, which also roughly is the distance between particles. If we take this number to be the wavelength of phonon modes, then the attendant frequency is 0.45/ps. This shows that phonon modes in the system where wnp = 0.05 begin to suffer frustration due to the finite system size (or, in a similar vein, due to the presence of nanoparticles). In essence this means that modes contributing to the upper left panel in Fig. 7 strongly scatter from the particles. The other panels show much more localized modes, associated with the molecular vibrations of the base fluid affected by the nanoparticles. If we carry this over to the experiments, then it is this effect of the nanoparticles on the molecular vibrations of the base fluid, which appear to be the cause for the observed enhancement of the specific heat.

Discussion

Nanofluids, consisting of a base liquid doped with small amounts of nanoparticles (often less than 1% wt.), frequently exhibit specific heat capacities, which cannot be explained by mere additivity of the heat capacities of the constituents. In the experiments, where this has been studied, apparently the effect is the same in the solid as well as in the liquid phase. If for instance the specific heat capacity is found to increase in the liquid phase, then there is a corresponding increase in the solid phase as well. The same holds true when the specific heat capacity is diminished by adding nanoparticles. Looking at the specific heat capacity as function of nanoparticle concentration, a local maximum is observed between 0.5 to 1.5%. However, quite generally the measurements are plagued by considerable scatter - sometimes comparable to the effect itself. Surprisingly this can be true also for the pure base fluid. When nanoparticles are present, the scatter may in part be attributed to the tendency of the nanoparticles to aggregate and an attendant lack of equilibration. In the field of filled elastomers, where the underlying physics has much in common with the physics of nanofluids, this aggregation or rather flocculation is well known and the subject of extensive research^54^. One particular information, which one can carried over from rubber research, is that there is no theory thus far describing the specific interaction between the matrix material, polymers there and base liquids here, with the different types of nanoparticles. A second piece of information is the long-range hydrodynamic interaction between the particles and the surrounding matrix material.

Recently the specific heat capacity of nanofluids has been studied via molecular simulation. The authors do indeed find an enhanced heat capacity, but the reason for this is not uncovered. We also do find that SiO2 nanoparticles enhance the heat capacity of liquid KNO3 including a maximum at low nanoparticle concentrations. By studying the effect of nanoparticles of the liquid’s spectral density distribution, we conclude that the presence of the particles causes additional low frequency vibrational modes. However, these modes require high particle densities and thus cannot be responsible for the ‘hump’-like feature in the heat capacity at low particle densities. And it is this ‘hump’-like feature, which we associate with the experimentally observed specific heat capacity maximum at low particle concentrations. In addition to the aforementioned low frequency modes we find that the presence of nanoparticles enhances existing molecular modes in the base liquid. It appears that this is causing the observed ‘hump’ of the specific heat capacity. Nevertheless, it still remains unclear how exactly the particles affect the molecular modes, e.g., by distortion of the intermolecular potentials.

The differential spectral densities do not permit to probe the gap in k-space, which is one of the governing factors in the phonon theory of liquids. However, the heat capacity increase in doped nanofluids is present also in the solid phase. This means that the effect the nanoparticles may have on the k-space gap should not be the true cause of the observed heat capacity enhancement. So how about anharmonicity? Unfortunately, in this work, we do not find a significant effect of the nanoparticles on αP aside from a near to linear decrease with increasing wnp. A distortion of the intermolecular potentials should however affect αP. On the other hand, the effects observed in the present system are small and somewhat difficult to quantify. Future work therefore should focus on systems, which experimentally exhibit much larger enhancement of the specific heat capacity.

Methods

Simulation details

Here we use the molecular dynamics simulation technique carried out using the widely used computation package LAMMPS^55^. The force field for KNO3 and its parameters were taken from ref.^56^ (see also ref.^57^). All intermolecular interactions are described by a Buckingham/Coulomb potential, i.e.,

with a cutoff radius of 9 Å. In the case of the Coulomb interaction long-range corrections are calculated using the PPPM method^58^ with an accuracy of 10^−6^. All bonded interactions, i.e., bond, angle and an improper dihedral potential keeping the nitrate planar, are described via harmonic potentials. Figure 8 shows an example of a typical system configuration. We model a single SiO2-particle immersed in liquid KNO3 inside a cubic volume of linear dimension L, applying periodic boundaries. This approach is analogous to two previous simulation studies, i.e., ref.^59^, focussing on Cu-nanoparticles in water, and ref.^35^, focussing on SiO2 in sodium nitrate, potassium nitrate and lithium nitrate. Nevertheless, we have considered a series of selected systems of size 2L, while keeping the concentration of nanoparticles constant, in order to check for finite size effects. Note that L and the diameter of the nanoparticle, D, are related via

Figure 8 System configuration showing the actual simulation cell surrounded by its replicas corresponding to a nanoparticle weight fraction wnp = 0.05.

The density of the nanoparticles, ρnp, is roughly 2.3 g/cm^3^. The density of the base salt, ρliq, is about 1.73 g/cm^3^. Finally, wnp is the mass fraction nanoparticles in the nanofluid. Focussing on a single nanoparticle rather than on a collection of nanoparticles immersed in a liquid, has practical reasons. Because the particle mass fractions of interest are small, modeling a collection of nanoparticles would either be prohibitively time consuming or require the particles to be unrealistically tiny. In addition, nanoparticles are known to aggregate, which probably contributes significantly to the scatter observed on the same system when it is studied by different experimental groups. This aggregation is a slow process, even on the time scale of real experiments, and inevitably causes a great deal of uncertainty in the simulations due to the limited statistics.

Our spherical SiO2-particles, with radii r, are cut from β-cristobalite. Subsequently the truncated valences are saturated with OH-groups. Even though industrial silica consists of amorphous particles, there is justification for this approach as has been discussed in previous work^60,61^. In the following we use the partial charges qSi = 1.91, qO = −0.9352, and qH = 0.4238. Otherwise the silica force field is modeled according to ref.^62^, using the reparametrization procedure in ref.^63^. Interactions between the silica and the salt are modeled using the Lorentz-Berthelot mixing rules^37^. Additional details can be found in ref.^64^ as well as in the contexts of specific results.

Selected theoretical concepts of heat capacity in liquids

Figure 9 shows cP vs. T for an eutectic salt mixture (NaNO3/KNO3 (60:40)) doped with SiO2 nanoparticles at 1% wt. (based on Figs 2 and 3 in Dudda and Shin^20^). Even though here the specific heat capacity enhancement is rather large, the figure is an overall typical example of the temperature dependence of cP in the pure and in the doped system. The immediate solid-to-liquid transition region is excluded. Addition of SiO2 nanoparticles increases cP in both phases and also increases the slope of the lines. Notice that qualitatively the solid state is not much different from the liquid state, suggesting that the underlying mechanism is at least similar in both phases.

Figure 9 cP vs. T for an eutectic salt mixture (NaNO3/KNO3 (60:40)) doped with SiO2 nanoparticles at 1% wt. from ref.^20^. Dashed pure salt mixture; solid salt mixture doped with nanoparticles of 60 nm diameter.

Overall the behavior of cP in the neat salt mixture is very much the same as in most common liquids, like water for instance, and even polymer melts below and above the glass transition (cf. Figure 2 in ref.^65^). The pronounced increases of the heat capacity, comparing its value on both sides of the transition (for instance 2 J/gK for water), may be understood qualitatively in terms of particles oscillating in ‘cages’. This idea dates back to Eyring et al. in the 1930 and was developed in the context of so-called cell models of liquids (cf. ref.^66^; chapter 4). On the solid side the molecules are more tightly constrained and corresponding oscillation frequencies are higher than on the liquid side. Therefore the Debye temperature is quite high on the solid side and (a fraction of) the vibrational modes are not fully excited below the transition temperature, resulting in a lower heat capacity of the solid phase.

Inspection reveals that the specific isochoric heat capacity of certain common liquids just above the melting transition is close to the Dulong-Petit (or high temperature) limit, i.e., cV = RNf/(2mmol). R of course is the gas constant, Nf is the number of degrees of freedom, which here is close to 6nA, where nA is the number of atoms per mole substance, and mmol is the molar mass. This appears to be true for most of the base fluids in the present context.

Here we are interested in the behavior of the heat capacity far from any phase transition. At a second order phase transition, for instance, the specific heat capacity diverges due to the divergence of the fluctuation correlation length ξ. Note that quite generally^67^ we expect

(in the continuum limit), where δe(r) is the local energy density fluctuation at a distance r from the origin (in an isotropic system). At a first order transition the specific heat capacity remains finite but shows a strong temperature dependence as well. Far from a phase transition the present nanofluids typically show a weak linear T-dependence over the entire range of temperatures covered by the experiments. Focusing momentarily on the solid state, this observation can be explained by considering the anharmonicity of the molecular interactions. If we express the free energy F of the solid at high temperatures as F=Fo+3(RnA/mmol)Tln(ℏω¯/kBT)^68^, where the second term is the phonon contribution, then, via E = F − T∂TF|V, we find E=Eo+3(RnA/mmol)T(1−T∂Tlnω¯|V). Note that ω¯ is defined via the geometric average ω¯3nA=∏i=13nAωi, where ωi are the normal mode frequencies of the solid. It can be shown^69^ that

and thus

where αP is considered independent of temperature. To see this we study the following two-step process in the P-T-plane. In a first step a system undergoes a thermal expansion at constant pressure, leading to a volume increase δV1 = αPVδT. In a second step the system’s volume is reduced by increasing the pressure at constant temperature, i.e., δV2 = −(V/B)δP|T. Here B is the system’s bulk modulus. We require the net volume change to vanish, i.e., 0 = δV1 + δV2 or αPBδT|P = δP|T. Finally we use V/B = −δV/δP|T once again. Defining a δB via δBV = BδV, we have δB = −δP|T and thus −αPBδT|P = δB or −αPB = δB/δT at constant volume. What this describes is the softening of the modulus by δB, which allows the volume to remain constant when the temperature increases by δT. The final ingredient is the assumption B∝ω¯2. On a very basic level this may be understood in terms of a harmonic oscillator for with k ∝ ω^2^, where k is the spring constant. Inserting this proportionality into the previous equation yields Eq. (8). To the same level of approximation as cV in Eq. (9) we find for cP

where γ is the Grüneisen parameter, i.e., γ=V∂P/∂E|V=∂P/∂T|V∂T/∂E|V=VαPB/CV. Note that (10) simply follows via CP = CV + VTBα^2^ in conjunction with (9). Finally it is worth pointing out that (9) and (10) can be obtained based on the ‘law’ of corresponding states. For instance the free enthalpy is expressed in terms of the ratio T/ωD, where ωD is the Debye frequency, G = Go(P) + ωDf(ωD/T)^68^. From this thermodynamic potential we may calculate the heat capacities using ∂T ln ωD|V = −αP/2. This essentially extends the approximate validity of (9) and (10) to some extend.

Equation (9) is a limiting case of a phonon theory of liquid thermodynamics developed by Trachenko et al.^36,50,70^. According to its key idea, the liquid acquires solid-like behavior at frequencies exceeding a certain inverse relaxation time, the Frenkel frequency ωF, where the liquid may support shear waves (Recent work has shown that there is a gap in momentum space rather than a frequency gap^49^, which however does not affect the thermodynamic functions derived previously). The authors develop an energy expression of the liquid combining contributions from longitudinal phonons, transverse (shear) phonons beyond ωF, and diffusion below ωF. Expanding their result, i.e., Eq. (5) in the above reference, in the high temperature limit, neglecting again the temperature dependence of αP, yields

where

The Frenkel frequency, ωF, in ref.^36^ is identified with the inverse relaxation time, τM, of the Maxwell model for viscoelasticity, i.e., ωF ≈ 1/τM and τM = η/G∞. Here η is the viscosity and G∞ is the shear modulus at infinite frequency. In the ‘solid-limit’, i.e., ωF = 0, Eq. (11) reduces to Eq. (9). In the liquid state, decreasing ωF leads to an increase of cV. Note that every translational mode contributes kB/2 to the heat capacity, whereas a fully excited one-dimensional oscillator contributes kB. Thus, the conversion from diffusive to oscillatory motion tends to increase of cV. Increasing the temperature, however, leads to a decrease of ωF and therefore decreases cV. Notice also that dr/dT > 0. The factor (1 + αPT) will counteract this to some extend, but experiments typically yield a decreasing cVexp^36^.

At this point we may ask how the addition of nanoparticles will affect the specific heat capacity, i.e., cV as well as cP. As discussed previously^6^ it is not sufficient to use mass weighted averages, e.g.,

Here wnp is the mass fraction nanoparticles and the respective specific heat capacities of the two components are cP,liq and cP,np. Instead it appears necessary to consider modified expressions.

One such expression is the interacting mesolayer model developed in ref.^6^. The model assumes that the center of each particle coincides with the center of a spherical shell of radius R + Δ, where R is the particle radius (in general Δ ≫ R). This modified version of Eq. (13) is given by

with

and

ρ¯=ρnpwnp+ρliq(1−wnp) (Note that Eqs (15) and (17) are the corrected versions of the mistyped Eqs (15) and (17) in ref.^6^). Here κ is the specific heat capacity in the mesolayer divided by the specific heat of the neat liquid. Note that cP(wnp = 0) = cp,liq and cP(wnp = 1) = cP,np, i.e., Eq. (14) agrees with (13) in these limits. In between, however, Eq. (14) deviates from the linear interpolation. Increasing wnp from zero in the interacting mesolayer model yields a cP above the linear interpolation, when κmax > 1. According to the underlying idea, the heat capacity in a liquid shell of thickness Δ is enhanced due to the presence of the central particle. A further increase of wnp leads to a statistical overlap of the mesolayers surrounding each particle, which means that Y decreases from one to zero and κ changes from κ = κmax to κ = κmin. Notice that max and min refer to the che case when the mesolayers are strongly overlapping, which may affect the heat capacity inside the mesolayer. The model, however, does not specify the microscopic origin of these changes. Nevertheless, in ref.^6^ it is shown, that it fits the experimental data quite well. Here we employ Eq. (14) to cP as well as to cV. At this point we return to possible microscopic reasons for the change of the specific heat capacity in the presence of nanoparticles.

It is well known^71–73^ that adding particles to a liquid causes its viscosity to increase. This hydrodynamic ‘reinforcement’, which indeed is a long-range effect, should thus decrease ωF. This is because ωF ∝ 1/η, resulting in an increase of cV (as well as cP). Nevertheless, it is known that the addition of particles to an elastic matrix, usually an elastomer matrix, will increase its shear modulus^74^. Overall there is no net effect on τM = η/G∞. However, the Maxwell model is an oversimplification and therefore this may not hold in real system governed by relaxation time distributions rather than a single one.

It is interesting to mention that polymer melts and the liquids discussed here are quite similar. In particular, the addition of nanoparticles does affect cP in a polymer melt beyond what one expects on the basis of Eq. (13). This remains true even if the polymer is vulcanized, i.e physically cross linked. In ref.^65^ the authors observe a pronounced decrease of cP when small amounts carbon nanotubes are added, i.e., 1% wt. and 5% wt., both below and above the glass transition.

In summary, the effect of nanoparticles on the heat capacity, as well as on other physicochemical properties, of liquids is not colligative, i.e., it does not merely depend on the nanoparticle weight fraction. Even qualitatively it depends on particle type and size, possibly morphology. In other respects it is basic and almost universal, i.e., it affects the liquid and the solid phase alike and it is observed over a wide range of base fluids of small molecules as well as polymers. The above theoretical framework allows us to focus on three aspects - the generation of additional shear modes entering through the Frenkel frequency, anharmonicity of the intermolecular potentials, and particle (network) induced liquid structure.

Acknowledgements

The authors are grateful to the Deutsche Forschungsgemeinschaft for its support of this work through grant HE 1545/18-1.

Author Contributions

Sven Engelmann - computer simulations and model development Reinhard Hentschke - study supervision, model development and drafting of the manuscript.

Data Availability

The datasets generated during and analysed during the current study are available from the corresponding author on reasonable request.

Competing Interests

The authors declare no competing interests.

Footnotes

Contributor Information

Sven Engelmann, Email: s.engelmann@uni-wuppertal.de.

Reinhard Hentschke, Email: hentschk@uni-wuppertal.de.

References

Associated Data

Data Availability Statement

The datasets generated during and analysed during the current study are available from the corresponding author on reasonable request.