As we have shown in Chapter 2, the transport coefficients appearing in the entropy production of a simple, one component fluid are the bulk and shear viscosities and the thermal conductivity. We have already devoted considerable attention to the shear viscosity, because of its importance as a practical transport coefficient and as a test property for nonequilibrium statistical mechanics. The determination of the thermal conductivity by homogeneous nonequilibrium molecular dynamics techniques is also an interesting and subtle subject, leading to fundamental questions regarding the relationship between the synthetic heat field appearing in nonequilibrium molecular dynamics simulations and the temperature gradient, and the existence of “heat waves”. When binary and multicomponent systems are considered, new transport coefficients associated with diffusion and the Soret and Dufour effects appear, and additional questions regarding the application of linear response theory to equations of motion with phase-space compression arise.

The thermal conductivity can be determined by equilibrium molecular dynamics simulations using the Green–Kubo formula, by homogeneous molecular dynamics methods that introduce a synthetic “heat field”, or by inhomogeneous molecular dynamics methods that examine the heat flow between thermal reservoirs. Here, we will focus on the homogeneous nonequilibrium molecular dynamics methods, because the equilibrium methods have been described in detail elsewhere, and the reservoir methods suffer from the same problems as the “boundary-driven” methods for the determination of the shear viscosity. Reservoir methods work by directly inducing a temperature gradient which must be large due to the poor signal to noise ratio of the molecular dynamics methods. This large temperature gradient also inevitably leads to inhomogeneity of the density. In this case the value of the thermal conductivity obtained from the simulation is the average over a range of temperatures and densities. This is an even greater limitation if we wish to compute the thermal conductivity near the critical point or near a phase transition. Similar considerations apply to the determination of the diffusion coefficient of the binary fluid, and the Soret and Dufour coefficients. Fortunately, in all cases, there exist homogeneous nonequilibrium molecular dynamics algorithms for the determination of these transport coefficients. The equations of motion used in these algorithms generate the desired fluxes in homogeneous systems. In each case, there is a statisticalmechanical proof that the equations of motion lead to the correct determination of the corresponding linear transport coefficient.

A detailed microscopic description of nonequilibrium classical mechanical systems can be achieved through the methods of nonequilibrium statistical mechanics. Beginning with the development of the Green–Kubo relations in the 1950s, modern nonequilibrium statistical mechanics has now reached a high level of sophistication, culminating in time-dependent nonlinear response theory and the derivation of fluctuation theorems. Several different approaches to nonequilibrium statistical mechanics that have slightly different fundamental perspectives have been developed. Some examples of these differing approaches can be seen in the books by McLennan [37], Eu [38], Zwanzig [39], Zubarev [40] and Gaspard [41]. The methodology we will use closely follows Evans and Morriss [2]. It evolved in close connection with the development of many of the nonequilibrium molecular dynamics simulation algorithms discussed in this book, and it is particularly suitable for our discussion. Readers seeking a broader introduction to nonequilibrium statistical mechanics are advised to consult the references provided.

Fundamentals of Classical Mechanics

Equations of Motion

Before proceeding with our discussion of nonequilibrium statistical mechanics, we will introduce the main concepts of classical mechanics. This will allow us to establish our notation and some basic ideas before moving on to the more complex material that follows. The simplest and most familiar formulation of classical mechanics for a single particle is based on Newton's second law of mechanics,

where m is the particle's mass, r is its position and ris its velocity. The total force F T may be a function of the particle's position, velocity and may also be an explicit function of time. For a system of interacting particles, the equation of motion for particle i may be a function of the positions and velocities of all other particles in the system, as well as an explicit function of time. The Lagrangian and Hamiltonian formalisms are richer and more elegant formulations of classical mechanics.

Generalised Hydrodynamics

The realisation that classical Navier-Stokes hydrodynamics would fail at molecular length-scales is not new. It has long been expected that generalised hydrodynamics could play an important part in the prediction of transport properties of inhomogeneous fluids. In 1983 Alley and Alder [321] noted that generalised hydrodynamic models could be very useful to the development of molecular-level predictive tools. In generalised hydrodynamics, the transport coefficients are no longer regarded as constants, nor are they simply material properties of a fluid at the local thermodynamic state point. Transport now becomes a fully nonlocal property of fluids, and likewise all the transport “coefficients” now become nonlocal in both time and space. In other words, the transport “coefficients” are now replaced by kernels and the governing constitutive equations are integral functions (convolutions) over both space and time. The kernels themselves now have both a wavelength and frequency dependence. In fact, significant progress in the development of generalised hydrodynamics was made in the 1970s by Akcasu and Daniels [369] and Ailawadi et al. [370]. The books by Boon and Yip [371], Eu [372] and Hansen and McDonald [61] also treat this subject in considerable detail, and we refer readers to these references for a much more thorough theoretical foundation.

Our purpose is to show how generalised hydrodyamics can be used as a predictive tool for inhomogeneous fluids. In particular, we will show how classical Navier-Stokes hydrodynamics breaks down at molecular length scales when the spatial variation in the driving thermodynamic force (e.g. the strain rate) is significant over the range of molecular interactions. Such conditions can occur not only for highly confined fluids, but also for fluids under shock conditions [373–375], where moment and gradient expansions originally derived from gas kinetic theory are commonly employed [372, 376]. Nonlocal effects are also observed in shear-banding phenomena [377], micellar solutions [378], Brownian suspensions of rigid fibres [379] and jammed glassy systems [380] and may well play a significant role in turbulence.

We will first demonstrate where breakdown in classical theory occurs in both confined fluids and systems in which the curvature in the strain rate exceeds the width of the relevant transport kernel. We then express the governing constitutive equation in its generalised form and demonstrate several techniques to compute the transport kernels. Once the kernels are computed, it is a straightforward matter to use them to predict properties of interest.

We live in a world out of equilibrium – a nonequilibrium world. We are surrounded by phenomena occurring in nature, in industrial and technological processes and in controlled experiments that we can only understand with the aid of a theoretical framework that encompasses nonequilibrium processes. Our understanding of these phenomena is largely based on a macroscopic theory that starts with the balance equations for the densities of mass, momentum, energy and other macroscopic quantities. To solve these equations, it is necessary to introduce relationships based on experiments that relate the observable properties of materials to the variables that define their macroscopic state. These relationships may describe equilibrium or locally equilibrium states of the material and in this case they are called equations of state. But we also need other relationships that relate the fluxes of properties to the property gradients that drive them. These are called constitutive or transport equations. The main subject of this book is the study of these transport equations and the material properties, such as the transport coefficients that account for the differences in the behaviour of different substances, using molecular dynamics simulation methods.

The molecular dynamics (MD) simulation method was developed soon after the Monte Carlo (MC) method, for the purpose of studying relaxation and transport phenomena [9]. Both MC and MD employed periodic boundary conditions, in which the system of interest is assumed to be replicated periodically in all directions, to limit (but not totally eliminate) the effects of the finite system size. At first, applications of this new technique focused on the structure, dynamics and equations of state of equilibrium systems [10–12]. The development in the 1950s of the Green–Kubo formalism, relating linear transport coefficients to equilibrium fluctuations in the corresponding fluxes [13, 14], made it possible to use equilibrium simulations to study nonequilibrium properties. However these methods, based on the computation of time correlation functions, were difficult to apply to all of the transport properties except self-diffusion due to their large computational requirements in comparison to the computing power available at that time. In addition, they could only address transport processes in the linear regime, i.e. where the flux is directly proportional to the thermodynamic driving force. These factors motivated the development of nonequilibrium molecular dynamics (NEMD) methods.

In this chapter, we introduce homogeneous nonequilibrium molecular dynamics simulation techniques by discussing the theoretical background to the SLLOD equations of motion. When these equations of motion are used in conjunction with compatible periodic boundary conditions and a homogeneous thermostat, they provide a very robust, reliable and well-understood method for studying fluids subjected to homogeneous flows. Here, we introduce the SLLOD equations of motion for the simple case of atomic fluids. This provides the groundwork for our discussion of methods for simulating homogeneous flows of molecular fluids in Chapter 8.

The SLLOD Equations of Motion

Background

To conduct microscopic simulations of flows driven by boundaries, mimicking real physical systems (e.g. Couette or elongational flows) we must explicitly include the walls. This inevitably induces density inhomogeneities into the fluid. If one is interested in nano-confined flow, then this is an appropriate simulation strategy since spatial inhomogeneity needs to be explicitly included in the simulation. However, if one is concerned with computing bulk properties such as mass, momentum and heat transport coefficients that we do not want to be distorted by surface effects, then the explicit use of boundaries is inappropriate.

An alternative to using atomistic wall boundaries is to generate flow through a suitable implementation of periodic boundary conditions. The first and most popular method of inducing flow through the periodic boundary conditions employs the so called Lees-Edwards boundary conditions [15] to generate planar shear flow. In such a scheme, a simulation box is replicated in all directions by periodic images.

So far we have largely considered NEMD simulation methods for fluids undergoing homogeneous flows. The exceptions to this were the “method of planes” and volume averaging techniques, in which expressions were derived that allow us to compute densities and fluxes for inhomogeneous fluids for various geometries (see Chapter 4). In this chapter we will use some of these expressions to compute relevant properties of highly confined fluids under several different flow conditions. In addition, we will show how to implement appropriate equations of motion to faithfully model fluids subject to spatially periodic fields and fluids under extreme confinement. It is with the latter application in mind that NEMD techniques are particularly important in the field of nanofluidics.

For over 150 years classical Navier-Stokes hydrodynamics [315] has been a wonderfully successful theoretical tool for predicting the properties of fluids and gases under a large variety of conditions. Its success extends from describing the dynamics of galactic motion [316], the aerodynamics of flight [317], the hydrodynamics of substances from liquid water to dense polymer melts [28], and right down to the flows of fluids on the microscale [318, 319]. It has even been shown to be accurate down to nanoscale dimensions [96], as long as certain conditions are maintained. In a numerical study on the Lennard-Jones fluid, it was first clearly demonstrated that, for an atomic fluid confined by atomistic walls, the Navier-Stokes equations were valid down to confinement spaces as low as around 5 to 10 atomic diameters [96, 320]. Below this spacing the fluid becomes highly inhomogeneous in space so the assumptions of constant density, constant viscosity, etc. break down, as do the Navier-Stokes equations.

At such small length scales another significant problem arises: the transport properties of fluids become nonlocal in nature. Although this effect has implicitly been built into the theory of generalised hydrodynamics [321] (see also Chapter 11) it has only recently been validated when the spatial extent of variations in the velocity gradient of the fluid are of the order of the width of the viscosity kernel [322]. The kernel itself is a nonlocal material property of the system with both wavevector and frequency dependence [321, 323] and has been accurately computed and parameterised for the Lennard- Jones fluid [324].