Figure 6.3 shows one way of representing planar Couette flow in a periodic system. In the Figure we only employ 2 particles per unit cell. In an actual computer simulation this number typically ranges from about one hundred to possibly several tens of thousands. As the particles move under Newton's equations of motion they feel the interatomic forces exerted by the particles within the unit cell and by the image particles whose positions are determined by the instantaneous lattice vectors of the periodic array of cells. The motion of the image cells defines the strain rate, γ≡∂u_x∕∂y for the flow. The motion of the cell images is such that their individual origins move with an x-velocity which is proportional to the y-coordinate of the particular cell origin.

(6.23)

If the Reynolds number is sufficiently small and turbulence does not occur, we expect that the motion of image particles above and below any given cell will, in time, induce a linear streaming velocity u(r), on each of the particles within the cell.

If during the course of time, a given particle moves out of a cell it will be replaced by its periodic image. If the particle moves through a y-face of a cell (that is, through the planes y=0 or y=L) the replacing image particle will not have the same laboratory velocity, nor necessarily the same x-coordinate. This movement of particles into and out of the primitive cell promotes the generation of a stable linear streaming velocity profile.

Although there are jump discontinuities in both the laboratory coordinates and the laboratory velocities of particles between cells there is no way in which the particles can actually sense the boundaries of any given cell. They are merely bookkeeping devices. The system is spatially homogeneous. As we shall see those components of particle velocity and position which are discontinuous have NO thermodynamic meaning.

We have depicted the Lees Edwards boundary conditions in the so-called sliding brick representation. There is a completely equivalent deforming cube representation that one can use if one prefers (see Figure 6.4). We will mainly use the language of the sliding brick representation - our choice is completely arbitrary however.

We will now consider the motion of particles under Lees Edwards boundary conditions in more detail. Consider a simulation cube of side L, located so that the streaming velocity at the cube origin is zero (that is the cube 0<{x,y,z}<L. The laboratory velocity of a particle i is then the sum of two parts; a peculiar or thermal velocity c_i, and a streaming velocity u(r_i), so

(6.24)

Imagine that at t=0 we have the usual periodic replication of the simulation cube where the boundary condition is

(6.25)

(with the modulus of a vector defined to be the vector of the moduli of the elements). As the streaming velocity is a function of y only, we need to consider explicitly boundary crossings in the y direction. At t=0, r_i has images at r'_i at r_i+jL, and r''_i at r_i+jL. After time t the positions of particle i and these two images are given by

(6.26)

where c_i and y_i (and their images) are functions of time. Now by definition the peculiar velocities of a particle and all of its periodic images are equal, c_i=c'_i=c''_i, so that

(6.27)

Similarly we can show that

(6.28)

If r_i(t) moves out the bottom of the simulation cube, it is replaced by the image particle at r'_i(t)

(6.29)

or if r_i(t) moves out of the top of the simulation cube, it is replaced by the image particle at r''_i(t)

(6.30)

The change in the laboratory velocity of a particle is given by the time derivative of equations (6.29) and (6.30). These rules for imaging particles and their velocities are shown schematically in Figure 6.4.

There is a major difficulty with the boundary driven algorithm. The way in which the boundaries induce a shearing motion to the particles takes time to occur, approximately given by the sound traversal time for the primitive cell. This is the minimum time taken for the particles to realise that the shear is taking place. The boundary driven method as described above, therefore cannot be used to study time dependent flows. The most elegant solution to this problem introduces the SLLOD algorithm. We will defer a discussion of thermostats and the evaluation of thermodynamic properties until after we have discussed the SLLOD algorithm.

The Boundary Driven shear flow algorithm has a number of disadvantages, the principle one being its lack of contact with response theory. We will now describe two synthetic field algorithms for simulating any form of flow deformation. Historically the first fictitious force method proposed for viscous flow calculations was the DOLLS tensor method (Hoover et.al, 1980). This method can be derived from the DOLLS tensor Hamiltonian,

(6.31)

It generates the following equations of motion

(6.32)

These equations of motion must be implemented with compatible periodic boundary conditions. If the strain rate tensor has only one nonzero element and it is off-diagonal, the deformation is planar Couette flow and Lees-Edwards boundary conditions must be used. If the strain rate tensor is isotropic then the flow is dilational and the appropriate variation of Lees-Edwards boundaries must be used. Other flow geometries can also be simulated using these equations.

One can see from the first of the equations (6.32), that since is obviously a laboratory velocity, the momenta p_i are peculiar with respect to the low Reynolds number streaming velocity u(r)=r⋅∇u. We call this streaming velocity profile the zero wavevector profile. If the Reynolds number is sufficiently high for turbulence to occur, the p_i are peculiar only with respect to the zero wavevector profile. They will not be peculiar with respect to any possible turbulent velocity profiles.

From (6.32) the dissipation is easily shown to be

(6.33)

where P is the instantaneous pressure tensor (3.152), whose kinetic component is given in terms of the peculiar momenta p_i. Since the DOLLS tensor equations of motion are derivable from a Hamiltonian, the AIΓ condition is clearly satisfied and we see immediately from equations (6.33) and (5.73), that in the linear regime, close to equilibrium, the shear and bulk viscosities will be related to equilibrium fluctuations via the Green-Kubo formula (T.6.3). This proves that the DOLLS tensor algorithm is correct for the limiting linear regime. The linear response of the pressure tensor is therefore,

(6.34)

The DOLLS tensor method has now been replaced by the SLLOD algorithm (Evans and Morriss,1984b). The only difference between the SLLOD algorithm and the DOLLS tensor equations of motion involves the equation of motion for the momenta. The Cartesian components that couple to the strain rate tensor are transposed. Unlike the DOLLS tensor equations, the SLLOD equations of motion cannot be derived from a Hamiltonian.

(6.35)

It is easy to see that the dissipation function for the SLLOD algorithm is precisely the same as for the DOLLS tensor equations of motion. In spite of the absence of a generating Hamiltonian, the SLLOD equations also satisfy AIΓ. This means that the linear response for both systems is identical and is given by (6.34). By taking the limit γ→0, followed by the limit t→∞, we see that the linear shear viscosity can be calculated from a nonequilibrium simulation, evolving under either the SLLOD or the DOLLS tensor equations of motion. With, ∇u=ji(∂u_x∕∂y), and calculating the ratio of stress to strain rate we calculate,

(6.36)

From (6.34) we see that the susceptibility is precisely the Green-Kubo expression for the shear viscosity (Table 6.1). Because the linear response of the SLLOD and DOLLS tensor algorithms are related to equilibrium fluctuations by the Green-Kubo relations, these algorithms can be used to calculate the reaction of systems to time-varying strain rates. If the shear rate is a sinusoidal function of time, then the Fourier transform of the susceptibility gives the complex, frequency-dependent shear viscosity measured in viscoelasticity (§2.4 & §4.3).

If the strain rate tensor is isotropic then the equations of motion describe adiabatic dilation of the system. If this dilation rate is sinusoidal then the limiting small field bulk viscosity can be calculated by monitoring the amplitude and phase of the pressure response and extrapolating both the amplitude and frequency to zero (Hoover et.al.1980). It is again easy to see from (6.3.13) that the susceptibility for the dilation induced pressure change, is precisely the Green-Kubo transform of the time dependent equilibrium fluctuations in the hydrostatic pressure (Table 6.1).

Although the DOLLS tensor and SLLOD algorithms have the same dissipation and give the correct linear behaviour, the DOLLS tensor algorithm begins to yield incorrect results at quadratic order in the strain rate. These errors show up first as errors in the normal stress differences. For irrotational flows (∇u=(∇u)^T) so the SLLOD and DOLLS tensor methods are identical, as can easily be seen from their equations of motion.

We will now show that the SLLOD algorithm gives an exact description of shear flow arbitrarily far from equilibrium. This method is also correct in the high Reynolds number regime in which laminar flow is unstable. Consider superimposing a linear velocity profile on a canonical ensemble of N-particle systems. This will generate the local equilibrium distribution function for Couette flow, f₁

(6.37)

Macroscopically such an ensemble is described by a linear streaming velocity profile,

(6.38)

so that the second rank strain rate tensor, ∇u, has only one nonzero element, (∇u)_yx=γ. The local equilibrium distribution function is not the same as the steady state distribution. This is easily seen when we realise that the shear stress evaluated for f₁, is zero. The local distribution function is no more than a canonical distribution with a superimposed linear velocity profile. No molecular relaxation has yet taken place.

If we allow this relaxation to take place by advancing time using Newton's equations (possibly supplemented with a thermostat) the system will go on shearing forever. This is because the linear velocity profile of the local distribution generates a zero wavevector transverse momentum current. As we saw in § 3.8, the zero wavevector momentum densities is conserved. The transverse momentum current will persist forever, at least for an infinite system.

Now let us see what happens under the SLLOD equations of motion (6.34), when the strain rate tensor is given by (6.38). Differentiating the first equation, then substituting for ṗ_i using the second equation gives,

(6.39)

If the strain rate γ is switched on at time zero, and remains steady thereafter,

(6.40)

Thus is a delta function at t=0. Now consider subjecting a canonical ensemble to these transformed SLLOD equations of motion, (6.39). If we integrate the velocity of particle i, over an infinitesimal time interval about zero. We see that,

(6.41)

So at time 0⁺ the x-velocity of every particle is incremented by an amount proportional to the product of the strain rate times its y coordinate. At time 0⁺, the other components of the velocity and positions of the particles are unaltered because there are no delta function singularities in their equations of motion. Applying (6.41) to a canonical ensemble of systems will clearly generate the local equilibrium distribution for planar Couette flow.

The application of SLLOD dynamics to the canonical ensemble is thus seen to be equivalent to applying Newton's equations to the local distribution function. The SLLOD equations of motion have therefore succeeded in transforming the boundary condition expressed in the form of the local distribution function into the form of a smooth mechanical force which appears as a mechanical perturbation in the equations of motion. This property is unique to SLLOD dynamics. It is not satisfied by the DOLLS tensor equations of motion for example. Since one cannot really call into question, the validity of the application of Newtonian dynamics to the local distribution as a correct description of Couette flow we are lead to the conclusion that the adiabatic application of SLLOD dynamics to the canonical ensemble gives an exact description of Couette flow.

Knowing that the SLLOD equations are exact, and that they generate Green-Kubo expressions for the shear and bulk viscosities, provides a proof of the validity of the Green-Kubo expressions themselves. The SLLOD transformation of a thermal transport process into a mechanical one, provides us with a direct route to the Green-Kubo relations for the viscosity coefficients. From equation (6.35) we see that we already have these relations for both the shear and bulk viscosity coefficients. We also see that these expression are identical to those we derived in Chapter 4, using the generalised Langevin equation. It is clear that the present derivation is simpler and gives greater physical insight into the processes involved.

Compared to the boundary driven methods, the advantages of using the SLLOD algorithm in computer simulations are many. Under periodic boundaries the SLLOD momenta which are peculiar with respect to the zero wavevector velocity field, and are continuous functions of time and space. This is not so for the laboratory velocities v_i. The internal energy and the pressure tensor of the system are more simply expressed in terms of SLLOD momenta rather than laboratory momenta. The internal energy E is given as,

(6.42)

while the ensemble averaged pressure tensor is,

(6.43)

For simulations of viscoelasticity special measures have to be taken in the boundary driven algorithm to ensure that the time varying strain rate is actually what you expect it to be. In the SLLOD method no special techniques are required for simulations of time dependent flows. One simply has to solve the equations of motion with a time dependent strain rate and ensure that the periodic boundary conditions are precisely consistent with the strain derived by integrating the imposed strain rate γ(t).

Since the SLLOD momenta are peculiar with respect to the zero wavevector velocity profile, the obvious way of thermostatting the algorithm is to use the equations,

(6.44)

The thermostatting multiplier α, is calculated in the usual way by ensuring that .

(6.45)

The temperature is assumed to be related to the peculiar kinetic energy. These equations assume that a linear velocity profile is stable. However as we have mentioned a number of times the linear velocity profile is only stable at low Reynolds number, (Re=ρmγL²∕η).

In Figure 6.7 we show the shear viscosity of 2048 WCA particles as a function of strain rate. The fluid is close to the Lennard-Jones triple point. The reduced temperature and density are 0.722 and 0.8442 respectively. The simulations were carried out using the Gaussian isokinetic SLLOD algorithm. We see that there is a substantial change in the viscosity with shear rate. Evidently WCA fluids are shear thinning in that the viscosity decreases with increasing strain rate. It turns out that this is common to all simple fluids for all thermodynamic state points. Shear thinning is also a widely observed phenomenon in the rheology of complex molecular fluids.

The imposed shear causes a major change in the microscopic fluid structure. This is manifest in all the thermodynamic properties of the system changing with shear rate. In Figure 6.8 we see the internal energy of the fluid plotted as a function of strain rate. For reduced strain rates in the range 0-1.5, we see that both the shear viscosity and the internal energy change by approximately 50% compared to their equilibrium values. Furthermore the viscosity coefficient appears to vary as the square root of the strain rate while the energy appears to change with the 1.5 power of the strain rate. Over the range of strain rates studied, the maximum deviation from the functional forms is 2.5% for the viscosity, and 0.1% for the internal energy. There has been much recent discussion of the relation of these apparently non-analytic dependences to mode-coupling theory (see, Yamada and Kawasaki, 1973; Kawasaki and Gunton, 1973; Ernst et. al., 1978; Evans , 1983; Kirkpatrick, 1984, van Beijeren, 1984 and deSchepper et. al., 1986). It is clear that the final resolution of this matter is still a long way off.

One of the most interesting and subtle rheological effects concerns the diagonal elements of the pressure tensor. For Newtonian fluids (ie fluids characterised by a strain rate independent and frequency independent viscosity) the diagonal elements are equal to each other and to their equilibrium values. Far from equilibrium, this is not true. We define normal stress coefficients, η₀, η_-, (the so-called out-of-plane and in-plane normal stress coefficients) as,

(6.46)

(6.47)

Figure 6.9 shows how these coefficients vary as a function of γ^1/2 for the WCA fluid.The out-of-plane coefficient is far larger than the in-plane coefficient, except at very small strain rates where both coefficients go to zero (ie the fluid becomes Newtonian). These coefficients are very difficult to compute accurately. They require both larger and longer simulations to achieve an accuracy that is comparable to that for the shear viscosity. In terms of the macroscopic hydrodynamics of Non-Newtonian fluids, these normal stress differences are responsible for a wide variety of interesting phenomena (eg the Weissenberg effect see Rainwater et. al. (1985 a,b)).

If one allows the strain rate to be a sinusoidal function of time and one extrapolates the system response to zero amplitude, one can calculate the linear viscoelastic response of a fluid. Figure 6.10 shows complex frequency dependent shear viscosity for the Lennard-Jones fluid (Evans, 1980), at its triple point.

If one compares Figure 6.10 with the Maxwell model for viscoelasticity, Figure 2.4, one sees a qualitative similarity with the low frequency response being viscous and the high frequency response being elastic. The shape of the two sets of curves is however quite different. This is particularly so at low frequencies. An analysis of the low frequency data shows that it is consistent with a nonanalytic square root dependence upon frequency.

(6.48)

where , , are the real and imaginary parts of the viscosity coefficient. Since the frequency dependent viscosity is the Fourier-Laplace transform of the memory function (2.76), we can use the Tauberian theorems (Doetsch, 1961), to show that if (6.48) represents the asymptotic low frequency behaviour of the frequency dependent viscosity, then the memory function must have the form,

(6.49)

This time dependence is again consistent with the time dependence predicted by mode-coupling theory (Pomeau and Resibois,1975). However as was the case for the strain rate dependence the amplitude of the effect shown in Figure 6.10, is orders of magnitude larger than theoretical predictions. This matter is also the subject of much current research and investigation.

Similar enhanced long time tails have been observed subsequently in Green-Kubo calculations for the shear viscosity (Erpenbeck and Wood, 1981). Whatever the final explanation for these enhanced long time tails, they are a ubiquitous feature of viscous processes at high densities. They have been observed in the wavevector dependent viscosity (Evans, 1982a) and in shear flow of 4-dimensional fluids (Evans, 1984). The situation for two dimensional liquids is apparently even more complex (Evans and Morriss 1983a and Morriss and Evans 1989).