Measuring Heat Flux Beyond Fourier's law

TL;DR

This study uses nonequilibrium molecular dynamics to analyze how shear flow influences heat flux, revealing additional shear-dependent components and validating theoretical predictions about heat transfer beyond Fourier's law.

Contribution

The paper introduces a detailed molecular dynamics approach to quantify shear-induced heat flux components, extending Fourier's law to include shear effects in thermal conductivity.

Findings

Heat flux in flow direction varies linearly with strain rate.

Additional shear-dependent heat flux components are quantified.

Results are consistent with theoretical predictions and experimentally testable.

Abstract

We use nonequilibrium molecular dynamics (NEMD) to explore the effect of shear flow on heat flux. By simulating a simple fluid in a channel bounded by tethered atoms, the heat flux is computed for two systems: a temperature driven one with no flow and a wall driven, Couette flow system. The results for the temperature driven system give the Fourier's law thermal conductivity, which is shown to agree well with experiments. Through comparison of the two systems, we quantify the additional components of the heat flux parallel and normal to the walls due to shear flow. To compute the heat flux in the flow direction, the Irving-Kirkwood equations are integrated over a volume, giving the so-called volume average form, and they are also manipulated to get expressions for the surface averaged and method of planes forms. The method of planes and volume average forms are shown to give equivalent results for the heat flux when using small volumes. The heat flux in the flow direction is obtained consistently over a range of simulations, and it is shown to vary linearly with strain rate, as predicted by theory. The additional strain rate dependent component of the heat flux normal to the wall is obtained by fitting the strain rate dependence of the heat flux to the expected form. As a result, the additional terms in the thermal conductivity tensor quantified in this work should be experimentally testable.

Tables

Table 1. Table 1 : Comparison of λ 𝜆 \lambda predicted using parameterization obtained from results in system 1 at measured T l subscript 𝑇 𝑙 T_{l} and T w ​ a ​ l ​ l subscript 𝑇 𝑤 𝑎 𝑙 𝑙 T_{wall} , to the λ 𝜆 \lambda value in system 2 using the zero strain rate intercept.

${\rho }_l$	$T_l$	System 2	Eq. (52)	Eq. (52)
		intercept	$\lambda ({\rho }_l,T_l)$	$\lambda ({\rho }_l,T_{wall})$
$\mathrm{\hspace{0.33em}0.65}$	$\mathrm{\hspace{0.33em}0.9}$	3.50	4.00	3.37
$\mathrm{\hspace{0.33em}0.75}$	$\mathrm{\hspace{0.33em}0.95}$	5.04	5.60	4.89
$\mathrm{\hspace{0.33em}0.85}$	$\mathrm{\hspace{0.33em}1.0}$	7.03	7.65	6.77

Equations

$${{\mathbf{J}}_q}^{approx}=-{\mathbf{\lambda }}_{\text{eff}}\cdot \nabla T+\zeta \frac{\partial \dot{\gamma }}{\partial y}\mathbf{i}+2\xi \dot{\gamma }\frac{\partial \dot{\gamma }}{\partial y}\mathbf{j}.$$

$${\boldsymbol{\lambda_{\textrm{eff}}}}=\left[{\begin{array}[]{*{20}{c}}{\lambda+3{\lambda_{2}}{{\dot{\gamma}}^{2}}}&{-{\lambda_{1}}\dot{\gamma}}&0\\ {-{\lambda_{1}}\dot{\gamma}}&{\lambda+3{\lambda_{2}}{{\dot{\gamma}}^{2}}}&0\\ 0&0&{\lambda+{\lambda_{2}}{{\dot{\gamma}}^{2}}}\end{array}}\right].$$

$$\rho \left(\mathbf{r},t\right)$$

$$\rho \mathbf{v}\left(\mathbf{r},t\right)$$

$$\rho e\left(\mathbf{r},t\right)$$

$$\frac{\partial \left(\rho \mathbf{v}\right)}{\partial t}$$

$$\frac{\partial \left(\rho e\right)}{\partial t}$$

$$\mathbf{P}\left(\mathbf{r},t\right)=\sum _i^N{m}_i{\mathbf{c}}_i{\mathbf{c}}_i\delta \left(\mathbf{r}-{\mathbf{r}}_i\right)$$

$$-\frac{1}{2}\sum _i^N\sum _{j\ne i}^N{\mathbf{r}}_{ij}{\mathbf{F}}_{ij}{\int }_0^1\delta \left(\mathbf{r}-{\mathbf{r}}_i-s{\mathbf{r}}_{ij}\right)ds$$

$${\mathbf{J}}_q\left(\mathbf{r},t\right)=\sum _i^N{u}_i{\mathbf{c}}_i\delta \left(\mathbf{r}-{\mathbf{r}}_i\right)$$

$$-\frac{1}{2}\sum _i^N\sum _{j\ne i}^N{\mathbf{r}}_{ij}{\mathbf{F}}_{ij}\cdot \left[{\mathbf{v}}_i-\mathbf{v}\left(\mathbf{r},t\right)\right]{\int }_0^1\delta \left(\mathbf{r}-{\mathbf{r}}_i-s{\mathbf{r}}_{ij}\right)ds$$

$${\rho }_V$$

$${\left(\rho v\right)}_V$$

$${\left(\rho e\right)}_V$$

$${\vartheta }_i={\int }_{\Delta V}\delta \left(\mathbf{r}-{\mathbf{r}}_i\left(t\right)\right)d\mathbf{r}$$

$$={\int }_{-\infty }^{\infty }\vartheta \left(\frac{\mathbf{r}-{\mathbf{r}}_m}{\Delta \mathbf{r}}\right)\delta \left(\mathbf{r}-{\mathbf{r}}_i\left(t\right)\right)d\mathbf{r}=\vartheta \left(\frac{{\mathbf{r}}_i-{\mathbf{r}}_m}{\Delta \mathbf{r}}\right)$$

$$\begin{gathered} \vartheta \left(\frac{\mathbf{r}-{\mathbf{r}}_m}{\Delta \mathbf{r}}\right)= \\ \Pi \left(\frac{x-x_m}{\Delta x}\right)\Pi \left(\frac{y-y_m}{\Delta y}\right)\Pi \left(\frac{z-z_m}{\Delta z}\right) \end{gathered}$$

$$\Pi \left(\frac{x-x_m}{\Delta x}\right)=\left[H\left(x-x_{-}\right)-H\left(x-x_{+}\right)\right]$$

$${\int }_{\mathrm{\Delta V}}P{}^K(\text{r},t)\mathrm{dr}=\sum _{i=1}^N{m}_i{\text{c}}_i{\text{c}}_i{\vartheta }_i$$

$$\bar{\mathbf{v}}\left({\mathbf{r}}_m,t\right)\equiv \frac{{\int }_{\mathrm{\Delta V}}\rho v\mathrm{dr}}{{\int }_{\mathrm{\Delta V}}\rho \mathrm{dr}}=\frac{{\sum }_{i=1}^N{m}_i{\text{v}}_i{\vartheta }_i}{{\sum }_{i=1}^N{m}_i{\vartheta }_i}$$

$${\int }_{\mathrm{\Delta V}}P{}^{\varphi }(\text{r},t)\mathrm{dr}=-\frac{1}{2}\sum _{i,j}^N{\mathbf{r}}_{ij}{\mathbf{F}}_{ij}{\int }_{\mathrm{\Delta V}}{\int }_0^1\delta (\mathbf{r}-{\mathbf{r}}_{\mathbf{s}})\mathbf{ds}\mathrm{dr},$$

$${\int }_{\Delta V}{\int }_0^1\delta \left(\mathbf{r}-{\mathbf{r}}_s\right)dsd\mathbf{r}$$

$$={\int }_0^1{\int }_{-\infty }^{\infty }\vartheta \left(\frac{\mathbf{r}-{\mathbf{r}}_m}{\Delta \mathbf{r}}\right)\delta \left(\mathbf{r}-{\mathbf{r}}_s\right)d\mathbf{r}ds$$

$$={\int }_0^1\vartheta \left(\frac{{\mathbf{r}}_s-{\mathbf{r}}_m}{\Delta \mathbf{r}}\right)ds.$$

$${\ell }_{ij}\equiv {\int }_0^1{\vartheta }_{s}ds,$$

$$P{}_V(\text{r},t)=\frac{1}{\mathrm{\Delta V}}\left[\sum _{i=1}^N{m}_i{\text{c}}_i{\text{c}}_i{\vartheta }_i-\frac{1}{2}\sum _{i,j}^N{\mathbf{r}}_{ij}{\mathbf{F}}_{ij}{\ell }_{ij}\right].$$

$${\mathbf{J}}_q{}^K=\frac{1}{\mathrm{\Delta V}}\sum _{i=1}^N{u}_i{\text{c}}_i{\vartheta }_i.$$

$${\int }_{\mathrm{\Delta V}}v(\text{r},t){\int }_0^1\delta (\text{r}-{\text{r}}_s)ds\mathrm{dr}$$

$$={\int }_0^{1}v({\text{r}}_s,t){\vartheta }_{s}ds\approx \bar{\mathbf{v}}\left({\mathbf{r}}_m,t\right).$$

$${\mathbf{J}}_q{}_V=\frac{1}{\mathrm{\Delta V}}\left[\sum _{i=1}^N{u}_i{\text{c}}_i{\vartheta }_i-\frac{1}{2}\sum _{i,j}^N{\mathbf{r}}_{ij}{\mathbf{F}}_{ij}\cdot \left({\text{v}}_i-\bar{\mathbf{v}}\left({\mathbf{r}}_m,t\right)\right){\ell }_{ij}\right].$$

$${\mathbf{J}}_{qV}^K\left({\mathbf{r}}_m,t\right)=\frac{1}{\Delta V}\left[\sum _{i=1}^N{e}_i{\mathbf{v}}_i{\vartheta }_i-\bar{\mathbf{v}}\left({\mathbf{r}}_m,t\right)\sum _{i=1}^N{e}_i{\vartheta }_i\right].$$

$${\mathbf{J}}_{qV}^{\varphi }\left({\mathbf{r}}_m,t\right)=$$

$$-\frac{1}{\Delta V}\frac{1}{2}\left[\sum _{i,j}^N{\mathbf{r}}_{ij}{\mathbf{F}}_{ij}\cdot {\mathbf{v}}_i{\ell }_{ij}-\left(\sum _{i,j}^N{\mathbf{r}}_{ij}{\mathbf{F}}_{ij}{\ell }_{ij}\right)\cdot \bar{\mathbf{v}}\left({\mathbf{r}}_m\right)\right].$$

$${\oint }_A{\mathbf{P}}^T\cdot d\mathbf{A}={\int }_{\Delta V}\nabla \cdot \mathbf{P}d\mathbf{r}.$$

$${\int }_{\Delta V}\nabla \cdot {\mathbf{P}}^{K}d\mathbf{r}$$

$$=-\sum _{i=1}^N{m}_i{\mathbf{c}}_i{\mathbf{c}}_i\cdot \frac{\partial }{\partial {\mathbf{r}}_i}{\vartheta }_i,$$

$$\frac{\partial }{\partial {\mathbf{r}}_i}{\vartheta }_i=\frac{\partial }{\partial {\mathbf{r}}_i}\vartheta \left(\frac{{\mathbf{r}}_i-{\mathbf{r}}_m}{\Delta \mathbf{r}}\right)=$$

Full text

Abstract

We use nonequilibrium molecular dynamics (NEMD) to explore the effect of shear flow on heat flux. By simulating a simple fluid in a channel bounded by tethered atoms, the heat flux is computed for two systems: a temperature driven one with no flow and a wall driven, Couette flow system. The results for the temperature driven system give the Fourier’s law thermal conductivity, which is shown to agree well with experiments. Through comparison of the two systems, we quantify the additional components of the heat flux parallel and normal to the walls due to shear flow. To compute the heat flux in the flow direction, the Irving-Kirkwood equations are integrated over a volume, giving the so-called volume average form, and they are also manipulated to get expressions for the surface averaged and method of planes forms. The method of planes and volume average forms are shown to give equivalent results for the heat flux when using small volumes. The heat flux in the flow direction is obtained consistently over a range of simulations, and it is shown to vary linearly with strain rate, as predicted by theory. The additional strain rate dependent component of the heat flux normal to the wall is obtained by fitting the strain rate dependence of the heat flux to the expected form. As a result, the additional terms in the thermal conductivity tensor quantified in this work should be experimentally testable.

I Introduction

In his 1878 work, Joseph Fourier proposed the existence of a single coefficient of thermal conductivity to describe heat conduction in an infinite solid (Fourier and Freeman, 1878). The simple linear relationship between heat flux and temperature gradient is known as Fourier’s law, an empirical observation of great importance in many areas of science and engineering. However, there are cases when a single coefficient of thermal conductivity is not sufficient to model the full range of physics. Examples include systems under extremely high shear rates as observed in tribological applications, nanoscopic systems where interfacial effects dominate and highly inhomogeneous systems where local constitutive equations break down. In these cases, a detailed understanding of heat flux at the small scale is essential. Deviations from Fourier’s law are potentially significant in nanochannels, MEMS and NEMS, nanoscale cpu components and complex materials such as aerogels with nano-scale pores, but they are exceptionally difficult to measure. It is in this context that molecular dynamics (MD) simulations are uniquely placed to provide insight.

Reliable methods for obtaining the Fourier’s law thermal conductivity are well known, including the Green-Kubo (Green, 1954; Kubo, 1957), Evans heat flow algorithm (Evans, 1982) and sinusoidal transverse field (STF) methods (Gosling et al., 1973). The thermal conductivity can also be obtained by comparing temperature profiles obtained in inhomogeneous nonequilibrium molecular dynamics (NEMD) simulations with predictions of constitutive laws (Todd and Daivis, 2017).

When the thermodynamic forces are large, the fluxes may depend on them nonlinearly. Under these conditions, Curie’s principle no longer applies and the heat flux may depend on nonlinearly coupled combinations of the temperature gradient and the velocity gradient (Daivis and Evans, 1993). If the velocity gradient varies rapidly with position, nonlocal effects also become important and the second spatial derivative of the velocity field may also be required as a thermodynamic force (Baranyai et al., 1992; Daivis and Coelho, 2000). Both of these effects lead to a component of the heat flux vector parallel to the streaming velocity in addition to a modification of the usual perpendicular heat flow due to Fourier’s Law. These results can be summarised in a constitutive equation of the form given by Daivis and Coelho (Daivis and Coelho, 2000). For a system in steady planar shear with a velocity gradient of the form $\nabla{\mathbf{v}}=\dot{\gamma}\left(y\right){\mathbf{ji}}$, a truncated Taylor expansion of the heat flux vector in powers of the thermodynamic forces $\nabla{\mathbf{v}}$, $\nabla\nabla{\mathbf{v}}$ and $\nabla T$ gives

[TABLE]

In a flow with a uniform velocity gradient field, such as homogeneous planar shear flow, ${\partial\dot{\gamma}}/{\partial y}$ is zero and only the first term is relevant.

The thermal conductivity tensor for a homogeneously shearing fluid with the fluid velocity in the $x$ direction and the gradient in the $y$ direction, i.e. with $\nabla{\mathbf{v}}=\dot{\gamma}\left(y\right){\mathbf{ji}}$, is given by

[TABLE]

The three thermal conductivity coefficients are: $\lambda$, the linear Fourier’s law thermal conductivity, $\lambda_{2}$ the coefficient of a nonlinear shear rate dependent heat flux in the direction of the temperature gradient ($y$) and $\lambda_{1}$ the coefficient for the shear induced heat flux in the flow direction, parallel to the walls ($x$).

Note that the off diagonal terms are chosen to be negative to give positive values of $\lambda_{1}$, a different convention than used previously (Daivis and Evans, 1993). The linear dependence on shear rate for the off-diagonal and quadratic dependence for the diagonal components of the thermal conductivity tensor has been confirmed for a model of liquid butane by Daivis and Evans (1995).

The tensor character of the generalised thermal conductivity reflects the anisotropy of the fluid induced by the shearing velocity gradient field. While shear rate dependent contributions to the heat flux should in principle exist, it has been suggested that it would be very challenging to measure them experimentally (Todd and Evans, 1995). In fact very few measurements of anisotropic thermal conductivity due to shear flow exist Venerus et al. (2004).

Kinetic theory, using Bhatnagar-Gross-Krook and Grad 13-moment style schemes has also confirmed the tensor character of the thermal conductivity for a fluid in planar shear flow (Garzó, 1994), as well as the existence of a secondary heat flux parallel to the streamlines, driven purely by the gradient of the shear rate, in dilute gases under planar Poiseuille flow (Risso and Cordero, 1997). Using the STF method, the shear-rate gradient induced heat flux was shown by Baranyai et al. (1992) to result in as much as a 25% departure from the purely linear solution for the temperature profile. Heat flux in Poiseuille flow of a dense fluid was also obtained from NEMD methods by Todd and Evans (1997). The magnitude of the second coefficient was shown in shear roller flow as a function of wavenumber (Monaghan and Morriss, 1997) for 2D flow. Han and Lee (2004) also proposed a local form of the method of planes formalism Todd et al. (1995a, b); Daivis et al. (1996) and captured the parallel component of the heat flux in an NEMD simulation. The strain rate coupling was also explored by Menzel et al. (2017) using NEMD for a polymer solution which indicates a very small value of the strain rate coupling coefficient $\xi$ with a magnitude of approximately one in reduced units, in agreement with literature values.

In this work, we explore the heat flux in planar shear flow using inhomogeneous NEMD simulations. In order to accurately compute the contribution of the parallel component, we use a large molecular system with approximately a quarter of a million atoms over a million timesteps. To isolate the shear rate dependent contribution to the heat flux from other contributions, we compare two systems; one which applies only a temperature gradient by maintaining different temperatures at the two channel walls with no strain rate and a second system with planar shear flow which generates its own temperature gradient and equal temperatures at the thermostatted walls. By running a parameter study over a range of densities in both systems, we isolate the shear rate dependent components of the thermal conductivity as a function of shear rate at each density.

We consider wall driven planar shear flow and obtain the mass density, streaming velocity, pressure tensor and heat flux vector using the formulation of Irving and Kirkwood (1950), integrated over either a cuboidal volume in space or a small rectangular surface. Spatial integration or averaging is important as it allows terms involving the Dirac delta functions in the Irving-Kirkwood equations to be numerically computed (Smith et al., 2012). The fluxes can then be written in one of two forms; $i)$ the volume average (VA) form (Cormier et al., 2001; Hardy, 1982; Lutsko, 1988) and $ii)$ a surface averaged form, which is a generalisation of the method of planes (MOP) (Todd et al., 1995a, b; Daivis et al., 1996; Heinz et al., 2005) providing all components of the pressure tensor and heat flux vector on the localised surfaces of a volume (Han and Lee, 2004).

This paper is organised as follows: in Section II, the formulae for obtaining the heat flux are derived in both volume averaged and surface averaged forms. Next, the extended constitutive equations relevant to the current systems are outlined. In Section III the details of the two NEMD simulations, conducted on systems 1 and 2, are outlined. The results and discussion are presented next in Section IV. Finally, this work finishes with a summary and conclusions in Section V.

II Theory

Expressions for the pressure tensor and heat flux vector can be derived by following the Irving and Kirkwood (1950) procedure. They typically begin (Evans and Morriss, 1990) with definitions of the local instantaneous mass density $\rho$, momentum density $\rho{\mathbf{v}}$ and total energy density $\rho e$,

[TABLE]

where ${m_{i}}$, ${{\mathbf{v}}_{i}}$ and ${e_{i}}$ are the mass, velocity and total energy of atom $i$. The delta function (strictly speaking a generalised function) represents the atomic localisation of each property. The densities are substituted into the momentum and energy balance equations,

[TABLE]

leading to explicit expressions (Evans and Morriss, 1990; Todd and Daivis, 2017) for the local instantaneous pressure tensor $\mathbf{P}\left({{\mathbf{r}},t}\right)$ and heat flux vector ${{\mathbf{J}}_{q}}\left({{\mathbf{r}},t}\right)$,

[TABLE]

where $u_{i}$ is the internal energy of particle $i$, ${\mathbf{c}}_{i}={{{\mathbf{v}}_{i}}-{\mathbf{v}}\left({{\mathbf{r}_{i}},t}\right)}$ is its thermal or peculiar velocity and ${\mathbf{r}}_{ij}={\mathbf{r}}_{j}-{\mathbf{r}}_{i}$. The Irving-Kirkwood expressions for the local pressure tensor and heat flux vector are well established, but they are purely formal. They contain delta functions which must be integrated in order to generate numerical values for the fluxes. In this section, volume averaged and surface averaged expressions for the heat flux vector are derived by integrating the Irving and Kirkwood (1950) fluxes over either the volume or the surfaces of a three-dimensional cuboid in Cartesian coordinates. When this cuboid is sufficiently small we obtain good approximations to the local values of the pressure tensor and heat flux vector. This presentation is kept brief as the majority has been given previously (Smith et al., 2012; Heyes et al., 2014).

II.1 Volume average fluxes

To obtain the volume averaged form of the pressure tensor and heat flux vector, we integrate both sides of the flux expressions in space. The integration volume $\Delta V$ can be regarded as a control volume in the fluid dynamics sense (Hirsch, 2007; Potter and Wiggert, 2002). Integrating the densities and dividing by the volume gives the average densities

[TABLE]

within volume $\Delta V$. Integrating the delta function gives a value of $1$ for each particle within $\operatorname*{\Delta{\it{V}}}$ and zero for each particle outside that volume (Smith et al., 2012). This is represented by the spatial selector function $\vartheta_{i}$ which can be written as

[TABLE]

The $\vartheta$ function is a three dimensional generalisation of the rectangle function defined by Bracewell (1999),

[TABLE]

where the one dimensional rectangle function is defined by

[TABLE]

and $H\left({x}\right)$ is the Heaviside step function. Similar relations hold for the rectangle functions in the $y$ and $z$ directions. In the definition of the rectangle function, ${x_{m}}$ is the midpoint of the rectangle $x_{m}=(x_{-}-x_{+})/2$, and ${\Delta x}$ is its width, and it is assumed that ${x_{+}}>{x_{-}}$.

Although we are mainly concerned with the heat flux vector in this paper, it is worth considering the simpler case of the pressure tensor first. Because the kinetic part of the pressure tensor is essentially a density, its volume average is easily evaluated in this form,

[TABLE]

where $\operatorname*{{\textbf{c}}_{\mathnormal{i}}}$ is the peculiar velocity $\operatorname*{{\textbf{c}}_{\mathnormal{i}}}\operatorname*{\equiv}\operatorname*{\textbf{v}_{\mathnormal{i}}}-\operatorname*{\textbf{v}}(\textbf{r}_{i},t)$, and the average velocity $\operatorname*{\textbf{v}}$ can be evaluated from the locally averaged momentum and mass densities

[TABLE]

and $\textbf{v}_{i}$ is the total velocity of particle $i$, i.e. the sum of the thermal and streaming velocities, with the streaming velocity at the location of the particle $\textbf{r}_{i}$ while the velocity obtained from Eq. (16) is the average for a volume with centre at ${\mathbf{r}}_{m}$.

Next, let us consider the configurational part of the pressure tensor,

[TABLE]

where ${{\bf{r}}_{\lambda}}={{\bf{r}}_{i}}+s{{\bf{r}}_{ij}}$ is a point on the line between molecules $i$ and $j$ and we use the subscript $i,j$ on the summation to denote sum over all $i$ and all $j$ with $i\neq j$. Swapping the order of integration, inserting the three dimensional rectangle function to limit the range of integration and using the sifting property of the delta function, we obtain

[TABLE]

The integrand on the right hand side is denoted by $\vartheta_{s}$, a function which is one if the point ${\bf{r}}_{s}$ is inside $\Delta V$ and zero otherwise, with the whole right hand side typically expressed as,

[TABLE]

where the integral is equal to the fraction of the line between particles $i$ and $j$ that is enclosed by the volume $\operatorname*{\Delta{\it{V}}}$. The average of the local pressure over a small volume $\operatorname*{\Delta{\it{V}}}$ is $\operatorname*{\textbf{P}}_{V}=\left(1/\operatorname*{\Delta{\it{V}}}\right)\int_{\operatorname*{\Delta{\it{V}}}}\operatorname*{\textbf{P}}(\textbf{r},t)\operatorname*{\it{d}\textbf{r}}$. Inserting the microscopic expressions, the volume averaged pressure tensor is then,

[TABLE]

To obtain the volume averaged form of the heat flux vector, we follow a similar procedure. Similarly to the pressure tensor, the average of the kinetic part of the local heat flux vector Eq. (II) over a small volume $\operatorname*{\Delta{\it{V}}}$ is,

[TABLE]

To compute the volume average of the configurational part of Eq. (II), we change the order of integration, introduce the three dimensional rectangle function and use the properties of the delta function to find

[TABLE]

where ${\mathbf{\bar{v}}}\left({{{\mathbf{r}}_{m}}},t\right)$ is the volume averaged velocity at $\mathbf{r}_{m}$, discussed in the appendix section A. This assigns the velocity contributions based on interaction length inside a volume, as discussed further in the appendix. Instead of evaluating the velocity integrated along the line, the velocity is approximated as the average value in the same volume in which we have evaluated the pressure. For simplicity, as we are using volume averages for other quantities, we assume the fluid velocity is also constant inside $\Delta V$, so we can write

[TABLE]

Error will be introduced into the computation of the heat flux vector if this relation is used when $\operatorname*{\textbf{v}}$ is not constant, but the error can be minimised by reducing the dimensions of $\Delta V$.

For practical computation it is convenient to split the heat flux into kinetic and configurational parts, so the total volume average heat flux is $\int_{\operatorname*{\Delta{\it{V}}}}\operatorname*{{\bf{J}}_{\it{q}}}\operatorname*{\it{d}\textbf{r}}=\int_{\operatorname*{\Delta{\it{V}}}}\left(\operatorname*{{\bf{J}}_{\it{q}}}^{\!\!\!K}+\operatorname*{{\bf{J}}_{\it{q}}}^{\!\!\!\phi}\right)\operatorname*{\it{d}\textbf{r}}$, that is, using the sum of the volume averaged kinetic and configurational components. Then the volume averaged form of the kinetic part of the heat flux vector suitable for computation is

[TABLE]

If the streaming velocity is approximately constant within $\Delta V$, the configurational component of the heat flux vector should be well approximated by

[TABLE]

In this form, the average velocities, average pressure and other quantities can be accumulated during the simulation until satisfactory statistics are achieved and then both components of the heat flux vector can be calculated at the end of the run.

II.2 Surface averaged fluxes

When the fluxes are averaged over a finite surface, we obtain generalisations of the method of planes results of Todd et al. (1995a, b) as given in Han and Lee (2004) and Heinz et al. (2005). A mathematical expression for fluxes over a volume surface in terms of functionals was derived previously for the pressure tensor (Smith et al., 2012) and is extended here for the heat flux vector.

The integral of a tensor over a closed surface is related to the volume integral of its divergence through the divergence theorem. For the pressure tensor, which is a second rank tensor, this is

[TABLE]

This relationship gives us a convenient way to obtain the surface averaged fluxes from the local expressions for the pressure tensor, Eq. (II), and heat flux vector, Eq. (II). Applying the divergence theorem to the kinetic part of the pressure tensor, we find

[TABLE]

where the delta function identity $\frac{\partial}{{\partial x}}\delta\left({x-y}\right)=-\frac{\partial}{{\partial y}}\delta\left({x-y}\right)$ has been used.

The derivative of the control volume function $\vartheta_{i}$ defined in Eq. (II.1) and Eq. (14) is obtained by well known manipulations of the Heaviside functions,

[TABLE]

where ${\mathbf{i}}$, ${\mathbf{j}}$ and ${\mathbf{k}}$ are unit vectors in the $x$, $y$ and $z$ directions. The physical interpretation of this function is that it selects atoms that cross the surface. The $\mathbf{i}$ component is the $x$ surfaces of the control volume, a square in space located at the top $x_{+}$ or bottom $x_{-}$ surface bounded by $[y_{-},y_{+}]$ and $[z_{-},z_{+}]$. A similar interpretation can be applied for the $y$ and $z$ direction with the j and k components. Using the divergence theorem, we therefore have

[TABLE]

The left hand side, a nine-component tensor dotted with each of the six surface vectors, gives us three Cartesian components for each of the six faces of the rectangular box, making 18 terms in total. Each term of this expression can be matched to a corresponding term on the right hand side. For example, the term on the top $x$ surface of the rectangular box, $P_{xx}\Delta{A_{x+}}$ is matched to the $x_{+}$ face of the i component, i.e.

[TABLE]

In general, for the $\alpha\beta$ component of the surface averaged pressure tensor, with $\alpha_{\pm}$ denoting either top $\alpha_{+}$ or bottom $\alpha_{-}$ faces of the rectangular box, we have

[TABLE]

To evaluate the time averages of these components of the pressure tensor, we only need to apply the delta function identity,

[TABLE]

which expresses the Dirac delta function as the sum of a Dirac delta for each of the function’s roots, where the index $n$ denotes the time $t_{n}$ at which crossings occur. In the time averaging process, we integrate this function over time as the system evolves, and this will isolate and count molecular crossings of the $\alpha_{+}$ surface,

[TABLE]

with the sifting property moving the crossing time, $t_{n}$, into the rectangle functions which checks if the crossing is within the area $\Delta\beta\Delta\gamma$. Similar expressions to Eq. (32) and Eq. (33) can be applied for the $\alpha_{-}$ surface.

To evaluate the expression for the surface average of the configurational part of the local pressure tensor, we substitute it into the volume integral of the divergence theorem, to obtain

[TABLE]

which simplifies to

[TABLE]

This must be equal to the surface integral of the divergence theorem and can be written as

[TABLE]

where

[TABLE]

The final equation for practical computation of the $\alpha\beta$ component of the pressure tensor at the $\alpha_{+}$ surface of the rectangular box is

[TABLE]

where $\beta_{\alpha+}\equiv\beta_{+}-\beta_{i}-\frac{\beta_{ij}}{\alpha_{ij}}(\alpha_{+}-\alpha_{i})$. The two rectangle functions select only those $ij$ interactions whose $\mathbf{r}_{ij}$ vector intersects the $\alpha_{+}$ surface within the area $\Delta\beta\Delta\gamma$ and the sgn function selects only those interactions having one particle on each side of the $\alpha_{+}$ surface. In the limit of infinite planar surfaces, the last two rectangle functions are always equal to 1, since $\Delta\beta$ and $\Delta\gamma$ both become infinitely large, and the expression reduces to the method of planes result. A similar treatment can be applied to the $\alpha_{-}$ surface.

The corresponding version of the divergence theorem for the heat flux vector is

[TABLE]

The expression for the surface averaged kinetic part of the heat flux vector is analogous to the corresponding part of the pressure tensor,

[TABLE]

The configurational part of the heat flux vector is more complicated. Separating the particle velocity and fluid streaming velocity terms of the configurational part of Eq. (II) we have

[TABLE]

The first term is analogous to the configurational part of the pressure tensor and results in

[TABLE]

with $S_{ij}$ again defined by Eq. (37). The second term in Eq. (40), containing the fluid streaming velocity, is

[TABLE]

which can be integrated to give

[TABLE]

where ${\mathbf{v}}\left({{{\mathbf{r}}_{x+}}}\right)={\mathbf{v}}\left({{{\mathbf{r}}_{i}}+\frac{{{x_{+}}-{x_{i}}}}{{{x_{ij}}}}{{\mathbf{r}}_{i}}}\right)$ is the fluid velocity at the point where the $\mathbf{r}_{ij}$ vector intersects the $x_{+}$ plane and similarly for ${\mathbf{v}}\left({{{\mathbf{r}}_{x-}}}\right)$, etc.

Inserting Eq. (41) and Eq. (43) into Eq. (40) and using the divergence theorem for the left hand side, we find that the final result for practical computation of the surface heat flux vector is

[TABLE]

This result shows that the relative velocity that appears in the expression for the surface averaged configurational part of the heat flux vector is the particle velocity minus the streaming velocity at the point where the $\mathbf{r}_{ij}$ vector intersects the $\alpha_{+}$ plane. This is again the plane peculiar velocity discussed by Todd et al. (1995b) Again, a similar expression exists for the $\alpha_{-}$ plane.

As with the volume average form, to calculate the surface averaged heat flux in a simulation, we decompose $\oint_{A}\operatorname*{{\bf{J}}_{\it{q}}}\cdot d\textbf{A}=\oint_{A}\left(\operatorname*{{\bf{J}}_{\it{q}}}^{\!\!\!K}+\operatorname*{{\bf{J}}_{\it{q}}}^{\!\!\!\!\phi}\right)\cdot d\textbf{A}$ into kinetic and configurational parts, obtaining the kinetic component from

[TABLE]

and the configurational component from

[TABLE]

where velocity $\bar{\mathbf{v}}$ is again from definition Eq. (16) extrapolated to the surface of the cell. This can again be computed in two parts that are combined at the end of the run. We have shown only the $x$ component, with similar forms for the $y$ and $z$ component heat fluxes on the other surfaces of the volume.

III Simulations

In this work, we use two different types of molecular simulation in order to isolate the three coefficients of thermal conductivity given in Eq. (2). Both cases simulate a wall bounded channel. In system 1, we apply a temperature gradient to obtain the Fourier’s law thermal conductivity coefficient $\lambda$ in the absence of a shear rate, recording heat flux over a range of densities. These heat flux results are compared to experimental results over a similar range of densities. In system 2, we simulate wall driven Couette flow over the same range of densities where we observed good experimental agreement for system 1. We obtain the coefficients of strain rate dependent thermal conductivity both parallel to the flow, $\lambda_{1}$, and normal to the flow, $\lambda_{2}$. In this section, details of the molecular simulation methodology common to both systems are first outlined before the key features of systems 1 and 2 are discussed.

III.1 Setup

In this work, we use a simple truncated and shifted Lennard-Jones potential (the so-called Weeks-Chandler-Anderson (WCA) potential Weeks et al. (1971)) in order to test the theoretical predictions:

[TABLE]

More realistic potentials and fluids have been avoided here because they could introduce further complexity due to rotational, non-Newtonian and coupled transport effects which might obscure the phenomena that are the focus of this work. The Lennard-Jones potential also has the advantage that it gives good agreement with experiments using liquid Argon. The equations of motion for the fluid are solved using the Leapfrog-Verlet scheme. The WCA cut off distance, $r_{c}=2^{\frac{1}{6}}$, was used, an assumption justified by previous observations that many transport properties are dominated by repulsive interactions (Todd and Daivis, 2017). Both systems have a timestep of $\Delta t=0.005$. Each simulation ran for a total of $1,000,000$ steps with samples taken every $25$ timesteps and results were written to disk every 25 samples (i.e. an output every $625$ timesteps). The system was initialised to an integer number of FCC units with the liquid region created by randomly removing atoms until the desired liquid density, $\rho_{l}$, was obtained. The walls remained as an FCC crystal lattice with density of $\rho_{w}=1.0$ and a thickness of $4$ atomic units ($L_{wall}=4$) from the top and bottom boundaries in the $y$ direction. The walls were tethered using the anharmonic potential from Petravic and Harrowell (2006) with spring constants for second, fourth and sixth power of displacement terms $k_{2}=0.0,k_{4}=5\times 10^{3}$ and $k_{6}=5\times 10^{6}$ respectively. This potential was chosen by Petravic and Harrowell (2006) to allow sufficient motion for momentum and energy interchange while preventing wall atoms from moving far enough to compromise heat flux or shear stress definitions. The system temperature was set to $T=1.2$ initially with the temperature of both walls controlled using a Nosé-Hoover thermostat applied only to the outer $L_{therm}=2$ reduced units of atoms in $y$, i.e.  half the wall is thermostatted in $y$ with the $2$ units closest to the liquid unthermostatted. The anharmonic tethering is purported to have better thermostatting characteristics, with a harmonic potential shown to impact the canonical distribution achieved by the Nosé-Hoover thermostat (Petravic and Harrowell, 2005).

The wall atoms have the following equations of motion:

[TABLE]

Here $Q_{\psi}$ is a coefficient which determines the Nosé-Hoover thermostat’s heat bath inertia, namely how strongly the wall temperature will tend to set point $T_{0}$. $\textbf{r}_{i_{0}}=\textbf{r}_{i}-\textbf{r}_{0}$, is atom $i$’s position relative to lattice site $\textbf{r}_{0}$ (the atom’s tethering location) where it slides at the wall speed $v_{w}$ in the streamwise $x$ direction with unit vector ${\bf{\hat{n}}}_{x}=[1,0,0]$. The two walls are thermostatted independently, i.e.  they are attached to separate heat baths where the strength of thermostatting is chosen to be proportional to $N_{therm}$, the number of atoms thermostatted in each wall, with $Q_{\psi}\operatorname*{\equiv}5\times 10^{-4}N_{therm}$.

Before averages are taken, the simulation is run for an initial period of time to reach steady state, which for system 1 is indicated by a time-independent approximately linear temperature profile, while in system 2 it is indicated by time-independent, approximately linear velocity and parabolic temperature profiles.

III.1.1 System 1

In system 1, we used separately thermostatted walls and densities varying from $\rho_{l}=0.4$ to $0.9$ with $T_{bot}=0.8$ at the bottom and $T_{top}=1.2$ at the top, giving a liquid temperature of $T_{l}=1.0$ on average. A range of average temperatures were also studied at fixed densities $\rho_{l}=\{0.65,0.75,0.85\}$, varying from $T_{bot}=0.5$ and $T_{top}=0.9$ to $T_{bot}=1.0$ and $T_{top}=1.4$ so that the temperature gradient remained constant. The system size is $L_{x}=31.75$, $L_{y}=47.62$ and $L_{z}=31.75$ which results in a liquid height of $39.62$, which for the case with $\rho_{l}=0.8$ has $N=37,902$ atoms, but varying for other densities. The domain is split into $22$ by $66$ by $22$ averaging bins in $x,y$ and $z$ respectively.

III.1.2 System 2

In system 2, as depicted in Fig. 1, both walls are thermostatted to the same temperature $T_{wall}=0.7$ with the walls atoms counter sliding at velocities set to $v_{w}=\pm 1$. The wall temperature was chosen to give a liquid temperature of $T_{l}=1$ on average. In order to give higher shear rates, a smaller channel height, $L_{y}$, was used than in the temperature study of system 1. As a result, the shear dependent terms in equations Eq. (53) and Eq. (54) are enhanced but with height chosen to give a liquid region of $15.81$, which is the smallest size where bulk behaviour has been shown to broadly match hydrodynamics (Travis and Gubbins, 2000). The domain size is expanded in other directions to improve statistics, with $L_{x}=139.69$, $L_{y}=23.81$ and $L_{z}=93.66$ so the number of atoms modelled is $N=259,359$ for the case with a density of $\rho_{l}=0.8$, again varying for other densities. The domain is split into $12$ by $256$ by $8$ averaging bins in $x,y$ and $z$ respectively.

IV Results and Discussion

In this Section, the results from the molecular simulations are presented and discussed, starting with system 1 in section IV.1 and then system 2 in section IV.2. Results from system 1 are presented first in order to obtain the thermal conductivity coefficient in the absence of shear.

IV.1 System 1 – Temperature Gradient

As outlined in Section III.1.1, a temperature gradient was induced by thermostatting each of the atomistic walls separately to different temperatures. The computed heat fluxes as a function of channel position are presented in Fig. 2, where we have shown the configurational component, $J_{q,y}^{\phi}$, kinetic component, $J_{q,y}^{K}$, and the resulting total heat flux calculated from the sum of both, $J_{q,y}=J_{q,y}^{K}+J_{q,y}^{\phi}$. The resulting heat flux is calculated for the channel using using Eq. (24)/ Eq. (45) for VA/MOP kinetic contributions and Eq. (25)/ Eq. (46) for VA/MOP configurational components. Fig. 2 shows results from both the VA and MOP approaches to calculating the heat flux, split into kinetic and configurational parts, with good agreement between all computations. The near-wall discrepancy between VA and MOP is attributed to density oscillations which are measured differently by the two methods, either as atoms inside a volume or crossings over a plane.

We obtain the temperature gradient $\partial T/\partial y$ from the derivative of a linear fit to temperature in the inner liquid region of the channel, starting $2.5$ reduced units away from the wall to avoid density oscillations and the Kapitza resistance at the surface. The heat flux, $J_{q,y}$, measured for the same inner liquid region, is then divided by temperature gradient to obtain the liquid’s thermal conductivity directly, $\lambda=-J_{q,y}/(\partial T/\partial y)$. Only the thermal conductivity of Fourier’s law $\lambda$, is relevant as the strain rate of system 1 is zero, so Eq. (1) and Eq. (2) simplify to

[TABLE]

This process is applied over the range of different liquid densities from $\rho_{l}=0.4$ to $\rho_{l}=0.95$.

The resulting thermal conductivity as a function of density is presented in Fig. 3 and compared to experimental results for liquid Argon from Roder et al. (1987) with both axes expressed in reduced units.

Good agreement is observed over the experimental range with discrepancies only apparent for $\rho_{l}$ below $0.6$ and above $0.85$. This is attributed to the limitations of the WCA model to reproduce experimental results for liquid Argon at more extreme densities. As a result, we limited the study to densities, $\rho_{l}$, in the range $[0.6,0.85]$ in the study of planar shear flow in system 2.

The study over a range of density was performed at a constant mean temperature of $T_{l}=1.0$, however the Fourier’s law thermal conductivity is also weakly dependent on temperature. In system 2, the mean liquid temperature was only controlled indirectly through the wall thermostats and can vary as shear rate changes. To quantify this, we performed a study between $T=0.7$ and $T=1.2$ at three of the densities which agree with experimental results $\rho=\{0.65,0.75,0.85\}$, shown in Fig 4. The experimental data from Roder et al. (1987) at comparable temperature and density is limited but all relevant values are included on Fig 4 and agree fairly well. The simulation results in Fig 4 are fitted using a linear function for each density. This temperature correction to the Fourier’s law thermal conductivity can then be included into a general fit, $\lambda_{f}$,

[TABLE]

where $T_{l}$ is the mean measured liquid temperature and $m_{T}=m_{T}(\rho_{l})$ is a gradient fitted to data at each of the three density values shown in Fig 4. Using second order polynomial fits to the data in Fig. 3 gives $\lambda_{f}^{\prime}(\rho_{l})=21.3\rho_{l}^{2}-14.2\rho_{l}+3.92$ and fitting a second order curve to the varying gradients obtained from Fig 4 yields $m_{T}(\rho_{l})=13.0\rho_{l}^{2}-17.0\rho_{l}+6.63$ which allows $\lambda_{f}(\rho_{l},T_{l})$ to be approximated for a given temperature and density.

IV.2 System 2 – Velocity Gradient

Having obtained the thermal conductivity in the absence of shear, we move on to system 2 as outlined in Section III.1.2, i.e. planar shear flow. We ran simulations over a range of strain rates by changing system size, $L_{y}=\{23.8,47.5,71.3,95.2\}$, as well as varying over the range of densities that matched to experimental results in system 1, $\rho_{l}=\{0.6,0.65,0.7,0.75,0.8,0.85\}$. The challenge in varying density and channel heights is that the system temperature and pressure also change in an uncontrolled manner, outlined in Fig 5, which summarises all the cases run in this section. The lines of constant density and constant applied strain rate are shown on Fig 5, highlighting that despite only varying either density or strain rate, there results a corresponding change in temperature and pressure. The Fourier’s law thermal conductivity coefficient, $\lambda$, at each temperature and density is estimated from Eq. (52) and shown in Fig 5 by the colouring of the points. This variation is observed despite thermostatting both walls to a temperature $T_{wall}=0.7$ in all cases.

The method for obtaining the components of heat flux is identical to the one used in system 1, based on Eq. (24)/ Eq. (45) for kinetic contributions and Eq. (25)/ Eq. (46) for configurational. However, in the presence of shear flow, we have additional components of heat flux. The $x$ component of heat flux, $J_{q,x}$, as a function of channel position is shown in Fig. 6 (a) and $y$ component of heat flux $J_{q,y}$ is shown in Fig. 6 (b). These are decomposed into the various contributions as before: kinetic heat flux, $\operatorname*{{\bf{J}}_{\it{q}}}^{\!\!K}$, and configurational flux, $\operatorname*{{\bf{J}}_{\it{q}}}^{\!\!\phi}$ as well as their sum $\operatorname*{{\bf{J}}_{\it{q}}}{}$. All simulations were run until the thermal and hydrodynamic steady state has been reached before properties were accumulated.

Both the VA and MOP approaches are used to calculate the heat flux and good agreement is observed between the two measurement techniques, as shown in Fig. 6. Previous studies have shown that the VA and MOP stresses are equivalent (Heyes et al., 2011) in the small bin limit and this observation is extended here to the energy equation and heat flux. The MOP form of stress is known to have worse statistics than VA (Hadjiconstantinou et al., 2003; Smith et al., 2017), a problem which is even more acute in the calculation of heat flux. Short averaging periods ($25$ records collected over $625$ timesteps) were found to result in a systematic error in the $J_{q,x}$ component of the MOP heat flux (when compared to the VA value). This may be due to very long correlations between $\operatorname*{\textbf{P}}$ and $\operatorname*{\textbf{v}}$. This error was found to disappear when averages from the entire one million timesteps simulation are used. In order to ensure agreement between MOP and VA approaches, averages of velocity, $\operatorname*{\textbf{v}}$, and the pressure tensor, $\operatorname*{\textbf{P}}$, are then taken over the entire simulation run length to calculate the correct stress heating term, $\operatorname*{\textbf{P}}\cdot\operatorname*{\textbf{v}}$. Interpolation of density and velocity measurements to the appropriate cell surface are also used to improve agreement between VA and MOP (the so-called ‘plane peculiar’ velocity problem as noted before and first pointed out by Todd et al. Todd et al. (1995b)). These are required as density and velocity measurements are obtained for atoms in a volume, so they represent the volume average about the centre of that cell, inducing a systematic error of order $\Delta y/2$ when calculating quantities at the cell surfaces. Having obtained the profiles for $J_{q,x}$ and $J_{q,y}$ in the sheared channel of system 2, we turn our attention to the calculation of the thermal conductivity coefficients $\lambda_{1}$ and $\lambda_{2}$.

For the case of planar shear flow with constant strain rate, Eq. (1) and Eq. (2) can be simplified to,

[TABLE]

where we have collected the strain rate dependent terms into single coefficients $\Lambda_{1}\operatorname*{\equiv}\lambda_{1}\dot{\gamma}$ and $\Lambda_{2}\operatorname*{\equiv}3\lambda_{2}\dot{\gamma}^{2}$. In the following section, we aim to use the values for $J_{q,x}$ and $J_{q,y}$ measured in the molecular simulation to obtain the two unknown thermal conductivity coefficients.

In order to obtain the coefficients $\Lambda_{1}$ and $\Lambda_{2}$, the expressions for $\partial T/\partial y$ and $\dot{\gamma}=\partial v/\partial y$ must be obtained. The near-wall density oscillations in the channel of system 2 affect the velocity and temperature results and prove especially problematic in calculating the numerical derivatives required for $\dot{\gamma}=\partial v/\partial y$ and $\partial T/\partial y$. In addition to density oscillations, numerical derivatives are also susceptible to noise in the data, which persist even for the large numbers of samples collected in these studies. To avoid these problems, least squares fits to the parabolic temperature profile and linear velocity functions are used and the resulting functional forms are differentiated instead. The profiles of velocity and temperature are shown in Fig. 7 (a) with fits to linear $v^{LS}=a_{1}y$ and quadratic, $T^{LS}=b_{1}y^{2}+b_{2}$, functions respectively, where $\{a_{1},b_{1},b_{2}\}$ are fitting parameters. The expressions, $dT^{LS}/dy$ and $\partial v^{LS}/\partial y=\dot{\gamma}^{LS}$, shown in Fig. 7 (b) as lines are obtained from the analytical derivative of the least squares fits. In doing this, we have decomposed the temperature profile, $T=T^{LS}+T^{\prime}$, into a purely quadratic part $T^{LS}$, as well as a separate density and noise dependent part $T^{\prime}$. Although the parabolic fit is good in Fig 7 (a), departures due to both noise and density stacking are evident in the derivative of Fig 7 (b). The derivative of the component of temperature due to density oscillations and noise can be expressed mathematically as $\partial T^{\prime}/\partial y=\partial T/\partial y-\partial T^{LS}/\partial y$. We replace $\partial T^{\prime}/\partial y$ with the less noisy gradient of density, $\partial T^{\prime}/\partial y\to\partial\rho/\partial y$, to improve the quality of the fit discussed in Appendix B.

The equations used in the fitting shown in Fig 7 $c)$ and $d)$ are therefore,

[TABLE]

where we have introduced two new coefficients $d_{x}$ and $d_{y}$ to be fitted. The importance of using $\partial\rho/\partial y$ as a proxy for $\partial T^{\prime}/\partial y$ is highlighted in Appendix B, Fig. 10, where different starting locations for the fitting profile can be seen to give very different measurement for heat flux. We have ignored nonlocal effects in the interface, which in principle should be included. However, the inclusion of nonlocal density coupling in shear flow is known to be a complex undertaking (Dalton et al., 2015; Glavatskiy et al., 2015) which is not warranted in this instance, and the approach used here provides a practical and convenient way to improve fitting in the presence of near-wall density stacking.

The approximate expressions for fluxes, Eq. (55) and Eq. (56), are fitted to the results for $J_{q,x}$, as shown in Fig. 7 (c) and $J_{q,y}$ in Fig. 7 (d), with these fits giving values for $\Lambda_{1}$ and $\lambda+\Lambda_{2}$ respectively. It seems reasonable to expect that we should be able get the coefficient of wall normal heat flux, $\Lambda_{2}$, by subtracting the coefficient of heat flux $\lambda$ obtained from Eq. (52) at the measured densities and temperatures. However, due to the highly coupled nature of the temperature and density highlighted in Fig 5, as well as variation in the channel itself, this was found to provide poor results. The values for $\Lambda_{1}$ and $\Lambda_{2}$ are found to be of similar magnitude to each other, $\sim 0.5$, but an order of magnitude smaller than the Fourier’s law coefficient $\lambda$, which is of order $\sim 5$. It is for this reason that these coefficients are difficult to obtain and temperature dependence becomes important. Further discussion of errors is provided in Appendix B. However, as the two heat flux contributions $\Lambda_{1}$ and $\Lambda_{2}$ are known functions of strain rate, we can use this dependence as another route to the strain rate dependent thermal conductivity coefficients.

To explore a range of strain rates, we compared coefficients from $L_{y}=23.81$ to three extra domain sizes, $L_{y}=47.5$, $L_{y}=71.3$, $L_{y}=95.2$, run with the same wall sliding velocity $v_{w}=\pm 1$, each at three densities $\rho_{l}=0.65$, $\rho_{l}=0.75$ and $\rho_{l}=0.85$. In addition to varying system size, varying wall sliding velocities was also explored at $L_{y}=23.81$ for $v_{w}=\{1.0,2.0\}$ at densities $\rho_{l}=\{0.65,0.75,0.85\}$. However, for $v_{w}=2.0$ we observe large pressure heating with temperature as high as $T\approx 4$ and, as a result, only varying system size was used in the strain rate analysis.

A plot of $\Lambda_{1}$ and $\Lambda_{2}$ against strain rate is presented in Fig. 8. The parallel heat flux term $\Lambda_{1}=\lambda_{1}\dot{\gamma}$ is linearly proportional to strain rate with a zero intercept, allowing us to use the fits of Fig. 8 to get a direct estimate of $\lambda_{1}$. The $\Lambda_{2}$ term is obtained by plotting the results for $\Lambda_{2}+\lambda$ from the fits in Fig. 7 $d)$, against the strain rate and fitting to $\lambda+3\dot{\gamma}^{2}\lambda_{2}$. The results are consistent with the expected quadratic strain rate dependence when $\lambda$ and $\lambda_{2}$ are allowed to vary freely in the fit. Subtracting the fitted intercept ensures a value of zero at the origin as required and gives us the results shown in Figure 8. The zero strain rate intercepts provide estimates of $\lambda$ for each density and the fit to the strain rate dependence provides a direct estimate of $\lambda_{2}$. The fits to results for all the densities studied appear to be similar, suggesting that the $\lambda_{2}$ coefficient is not very sensitive to density within the range of densities studied here.

The values of $\lambda$ measured in system 2 are compared to the predictions using Eq. (52) based on the liquid’s density and temperature in Table 1. As the wall is thermostatted to $T_{wall}=0.7$, in the absence of shear the liquid would be expected to equilibrate to the same temperature as the walls. The zero-strain intercept should therefore be the predicted $\lambda$ for this equilibrium case at $T_{wall}$, which is also shown in Table 1. The system 2 prediction of $\lambda$ is seen to be close to the one expected at the wall temperature, within $5\%$, but with a slightly greater value in all cases. The use of an average temperature over the liquid region or wall temperature is crude as both temperature and $\lambda$ would be expected to vary inside the channel. A more detailed analysis could take into account the varying thermal conductivity coefficient as a function of $y$, using the approach of Todd and Evans (1997).

A summary of results for the three coefficients of thermal conductivity $\lambda$, $\lambda_{1}$ and $\lambda_{2}$ is presented in Fig. 9. We use the strain rate $\dot{\gamma}^{LS}=\partial v^{LS}/\partial y$ to work out parallel heat thermal conductivity coefficients by simply dividing the VA and MOP values by the strain rate, $\lambda_{1}=\Lambda_{1}/\dot{\gamma}^{LS}$ at each density. As a cross check, we also use the fits obtained in Fig 8 over the range of strain rate studies, shown as crosses on Fig 9 with excellent agreement. The Fourier’s law thermal conductivity coefficients $\lambda$ are also presented in Fig. 9 for all studied density. In each case, the value is corrected based on the average liquid temperature using Eq. (52). The average liquid temperature measured in each system is as summarised in overview Fig. 5. Only $\lambda$ values for the highest strain rate case, $\dot{\gamma}=0.084$, are shown, other results are similar after correction.

Finally, the three $\lambda_{2}$ values obtained from the fits to strain rate of Fig 8 are also presented in Fig. 9. A full discussion of possible sources of error is included in Appendix B. Although the parallel heat flux has been previously observed in NEMD simulations of a dense fluid undergoing Poiseuille flow Han and Lee (2004), this is the first time that the values of both $\lambda_{1}$ and $\lambda_{2}$ have been evaluated directly from the heat fluxes. We further note that Han and Lee (2004) observed a stronger heat flux in the parallel direction, in contrast to our much weaker heat flux observed for planar Couette flow. In future studies it will be interesting to repeat these measurements for Poiseuille flow to see whether we obtain similarly strong parallel heat fluxes, and measure the resulting values of $\lambda_{1}$ and $\lambda_{2}$.

The measured values of $\lambda_{1}$ in Fig 9 are found to be of similar magnitude to that found by Daivis and Evans (1993) using a Green-Kubo like method, where a value of $\lambda_{1}\approx 2.2$ was obtained for a LJ atomic system with a density of $\rho_{l}=0.8$ and a temperature of 2.0. Daivis and Evans (1993) also found that $\lambda_{2}$ was approximately zero, in contrast to the results presented here. The reason for this is not presently understood, but the different state point, LJ potential cutoff and range of shear rates studied may all be partially responsible.

Given the good agreement for $\lambda$ compared to liquid Argon experimental data, the presented measurements of $\lambda_{1}$ and $\lambda_{2}$ should be testable against laboratory experiments. The high shear rates required and small magnitudes of this secondary heat flux would make such measurements extremely challenging for this system. It is for this reason that molecular simulation remains an invaluable tool for studies of this kind.

V Conclusions

In this work we have used the Irving and Kirkwood (1950) equations in integrated form to obtain expressions for heat flux, both as an average quantity in a volume (Volume Average, VA) and defined on surfaces (Method of Planes, MOP). By obtaining the heat flux on the various surfaces, the MOP approach is extended to provide heat flux parallel to the flow. The MOP and VA definitions are shown to give equivalent results for measured heat flux when properly averaged. From a non-equilibrium molecular dynamics (NEMD) simulation using a channel with a temperature gradient, the coefficient of heat flux from Fourier’s law is obtained and shown to have good agreement to experimental results. The Fourier’s law coefficient is parameterised over a range of densities and temperatur es. Boundary driven planar Couette flow is then used to explore the additional terms due to strain rate in the direction parallel and normal to the wall, as predicted by non-equilibrium thermodynamics.

An estimate of the second coefficient of thermal conductivity in the wall normal direction is obtained by fitting a quadratic in strain rate. This estimate of second coefficient is obtained at three different densities and appears to be broadly independent of density, within the range of densities studied in this work. The zero strain intercept provides a separate estimate of the Fourier thermal conductivity coefficient, which is found to be slightly lower than the parameterised form obtained from system 1. In addition, the results for the parallel component of heat flux appear to be very convincing. As a function of strain rate, it is shown to vary linearly as predicted by theory, allowing us to fit over a wide range of strain rates providing this coefficient as a function of density. The consistency of this coefficient from all simulations suggests this value could be used as a numerical prediction testable through experiments.

Appendix A Additional Derivations

In the Appendix, Section A has a discussion of the alternative form of the Irving and Kirkwood (1950) operator as an integral along a line, with the result that the velocity expressions can be defined along this line, as an average in a volume or as an average over a surface. In Section B, the various sources of error in this work are set out and the attempts to mitigate them discussed, including the use of density gradients in the fitting.

The pressure from Irving and Kirkwood (1950) is expressed as the difference between two Dirac delta functions. To derive an expression for the configurational part of the local, instantaneous pressure tensor, we will need to express the difference between two delta functions as a divergence. This can be done as follows. Using the fundamental theorem of calculus, we can write

[TABLE]

However,

[TABLE]

so we can write

[TABLE]

which is the required form as an integral along the line between the atoms. The volume integral of this expression is as given in Eq. (19), which occurs in quantities which depend on inter-particle interactions such as configurational pressure and heat flux. The heat flux includes velocity and it is necessary to write the averaged fluid velocity integrated along the line,

[TABLE]

Then the full expression for the volume averaged heat flux vector is

[TABLE]

Note this includes a contribution from the difference between the velocity $\operatorname*{\textbf{v}_{\mathnormal{i}}}$ of particle $i$ and the fluid velocity averaged over all points on the inter-molecular interaction line $\mathbf{r}_{ij}$ that are inside $\Delta V$, ${\mathbf{\bar{v}}}\left({{{\mathbf{r}}_{ij}}}\right)$. If $v$ is constant inside $\Delta V$, then the $\ell_{ij}$ on the top and bottom cancel and we are left with $\bf{v}$ at ${\bf{r}}_{m}$ as assumed in equations (23) and (25).

In order to get the expressions for the heat flux, we can convert the kinetic part of the heat flux vector in Eq. (23) into a form where the instantaneous streaming velocity is not required. We can write

[TABLE]

where the last line follows from

[TABLE]

In the main text, we use the velocity averaged along a line and in a volume when we consider heat flux in a volume; while for the surface or plane based measures, velocity is defined on the surface as follows,

[TABLE]

where this is the velocity on the a plane located at $x^{+}$ obtained from the average flux of atoms (Todd et al., 1995a) and the density, obtained as an average over a cell centred on $\bf{r}$, as $\sum m_{i}\vartheta_{i}$ should be extrapolated to the location of the cell surface. An alternative approach would be to collect density at the location of the surface,

[TABLE]

where the density at a plane can be computed using the method of Daivis et al. (1996). In practice, both approaches are used and compared, giving similar results in this work as the statistics collected are in steady-state temperature and shear driven flows. However, in more general cases such as turbulent flow (Smith, 2015) this distinction may become important.

Appendix B Uncertainty Analysis

We discuss the sources of possible error in this study, which include: (i) temperature variation, (ii) finite system size effects, (iii) departure from classical Navier-Stokes hydrodynamics and (iv) statistical errors.

First, consider (i) temperature due to shear heating and the thermostatting of only the wall atoms. Different density fluids have different average temperatures, with denser fluids generally hotter. In the worst case of $\rho_{l}=0.85$ and a channel of height $L_{y}=23.81$ in system 2, the mean fluid temperature increases to $T=1.25$, which is compared to the value of $T=1.0$ at the same density when obtaining $\lambda$ in system 1. The difference in heat flux due to the difference in temperature, estimated using NIST isochoric data at $\rho_{l}=0.85$ for liquid Argon (Perkins et al., 1991), would give a difference of $\sim 2.5\%$ in the value of $\lambda$, which could imply an almost $60\%$ error in the calculation of $\lambda_{2}$. This is a direct consequence of the relative magnitudes of the $\lambda$ and $\lambda_{2}$ coefficients, making isolation of these higher order heat flux terms difficult. It is also unclear if $\lambda_{2}$ would be expected to have the same functional dependence on temperature as $\lambda$. It should be noted that strain rate dependence is very difficult to decouple from temperature dependence without thermostatting the entire system; an approach which should be avoided as previous work has demonstrated that wall-only thermostatting is essential in order is preserve the correct system dynamics (Bernardi et al., 2010a, b). The use of temperature correction based on system 1 has been employed in this work in an attempt to solve this dependence.

Another possible source of error, (ii), is finite size effects. In the work of Hyżorek and Tretiakov (2016), the effect of channel height on measured heat flux shows that confinement effects cause a departure from the experimental Fourier’s law. This is observed for channels smaller than about $20$ in reduced units, especially at higher densities. Although Hyżorek and Tretiakov (2016) use a longer cutoff radius $r_{c}=5$ and an infinite wall-fluid interaction length (factors which potentially increase the impact of nano-confinement) these trends suggest that the value of $\lambda$ could be reduced by the smaller channel used in system 2 as compared to system 1. This would, in turn, suggest that the obtained $\lambda_{2}$ in this work may increase in error for high densities. The departures observed by Hyżorek and Tretiakov (2016) are of the order $1\%$ when applied to get $\lambda_{2}$ for low density systems ($\rho_{l}=0.4$), increasing to $10\%$ for the highest densities they considered ($\rho_{l}=0.7$). Extrapolation using a fit to digitised data from Hyżorek and Tretiakov (2016) suggest this may rise as high as $25\%$ in the highest density, $\rho_{l}=0.85$, presented in this work. Again, it is not clear if the $\lambda_{1}$ and $\lambda_{2}$ coefficients will be similarly impacted by nano-confinement. The use of larger systems, as large as $L_{y}=95.2$ reduced units in the strain rate study, help reduce this effect. However, large channels necessarily have low strain rates and so poor signal to noise ratio. Fitting over a range of strain rates is used to get values of $\lambda_{1}$ and $\lambda_{2}$.

The uncertainty induced due to (iii) departures from classical hydrodynamics include density oscillations, non-zero slip lengths, non-Newtonian effects and departure from expected velocity and temperature profiles. Indeed, a departure from a linear velocity profile is observed in the high density large channel cases (greater than $40$ reduced units), an effect especially apparent near the wall. As a result, the measured strain rate is taken not simply as $v_{w}/L_{y}$ but obtained as the average profile of velocity over the domain. This is the reason that the strain rates in Fig. 8 are different for the three different densities. No wetting parameter is used in this work but the higher density walls will potentially cause slip and Kapitza-like jumps (Kim et al., 2008). The impact of slip-length, both in temperature and velocity profiles, has been reduced by fitting only to the liquid region, i.e.  applying fitting to averaging bins located a distance from the wall. Slip is linked to density stacking effects discussed previously and so the inclusion of density in our fits also help to reduce this effect, as discussed in point (iv).

These measurements are also subject to (iv) large levels of noise with very weak signals, shown for example in Fig. 6 (a). These errors have been reduced as far as possible by collecting statistics for large atomistic systems over long runs. In addition, the use of least square fits to obtain derivatives of temperature and velocity improve statistics. However, the density stacking near the wall poses a problem for fitting least squares solutions to obtain coefficients of heat flux. To highlight the problem, consider Fig 10 $a)$ which shows the fitting parameter (i.e.  gradient of the line) returned when we start the fit at different bins in the channel.

In this case, a value of zero on the $x$-axis means a fit to all $170$ liquid bins, while $10$ for example excludes the $10$ bins either side of the channel closest to the wall. This means we would fit a line to only $150$ innermost ones. By including the part of the channel with density peaks, it is clear that we obtain completely different coefficients, by almost an order of magnitude in Fig 10 $a)$. This is a particular problem in this work given the very small magnitudes of the $\lambda_{1}$ and $\lambda_{2}$ coefficients compared to variations in density in this highly confined channel. As a result, it is almost impossible to work out a meaningful location to start fitting and using only the inner-most part of the channel gives very poor statistics. Instead, we improve the fit by replacing the term for the oscillating component of the temperature field $\partial T^{\prime}/\partial y$ with a term proportional to the gradient of density itself $\partial\rho/\partial y$, as the gradient of density correlates very well with the gradient of temperature but has much less noise as shown in Fig 11. Todd and Evans (1997) have shown that this is a good approximation near the centre of the channel when the pressure is constant, which is what we can expect when local thermodynamic equilibrium holds. In the interfacial region, the thermodynamics is more complex but the density and temperature gradients remain highly correlated.

By including the density derivative term in Eq. (55) and Eq. (56), the results as a function of fitting location become much less susceptible to density peaks, as shown in Fig 10 $b)$. The starting point for the fitting is chosen to be bin $11$, shown by the black line on Fig 10 $b)$ as this appears to be consistent with the average gradient in the channel and includes the majority of bins in the fit.