Computer simulations have been carried out for two different systems (Evans and Morriss, 1988). Two statistical mechanical systems were studied. The first was a system of 72 soft disks, (Φ=4ε(σ∕r)¹²), in two dimensions at a reduced density, ρ^*=ρσ² = 0.6928, a reduced temperature, T^* = kT∕ε = 1, and for a range of reduced strain rates, γ ^*= γ(m∕ε)^1/2σ = ∂u_x∕∂y (m∕ε)^1/2σ. The second system was studied more extensively. It consisted of 256 WCA particles. The system was three dimensional and the density was set to ρ^*=ρσ³ = 0.8442 while the temperature was T^* = kT∕ε = 0.722 (ie the Lennard-Jones triple point state).

In each simulation the direct NEMD average of the shear stress, pressure, normal stress difference and thermostat multiplier α, were calculated along with their associated transient correlation functions using typically 60,000 nonequilibrium starting states. For the three dimensional system each nonequilibrium trajectory was run for a reduced time of 1.5 (600 timesteps). Each 60,000 starting state simulation consisted of a total of 54 million timesteps made up of 2 x 15,000 x 600 timesteps at equilibrium and 4 x 15,000 x 600 perturbed nonequilibrium timesteps. The trajectory mappings described in §7.4 were used to generate the 4 starting states for the nonequilibrium trajectories.

In Figure 7.1 we present the results obtained for the reduced shear stress P_xy ^*=P_xy(σ²∕ε), in the 2 dimensional soft disk system. The imposed reduced strain rate is unity. The values of the shear stress calculated from the transient correlation function expression, (P_xy(T)), agree within error bars, with those calculated directly, (P_xy(D)). The errors associated with the direct average are less than the size of the plotting symbols whereas the error in the integral of the transient correlation function is approximately ±2.5% at the longest times. Although the agreement between the direct simulation results and the TTCF prediction is very good it must be remembered that the total response for the shear stress is the sum of a large linear effect which could be correctly predicted by the Green-Kubo formula and a smaller (~25%) nonlinear effect. Thus the statistical agreement regarding the TTCF prediction of the intrinsically nonlinear component of the total response is therefore approximately 10%.

The shear-induced increase in pressure with increasing strain rate (shear dilatancy) is an intrinsically nonlinear effect and is not observed in Newtonian fluids. The Green-Kubo formulae predict that there is no coupling of the pressure and the shear stress because the equilibrium correlation function, <Δp(t)P_xy(0)>, is exactly zero at all times. In Figure 7.2 we present the direct and transient correlation function values of the difference between the pressure p^*=p(σ²∕ε) and its equilibrium value, p₀ ^*, (Δp^*=p^*-p₀ ^*). The agreement between the direct average, and the value obtained from the transient correlation function expression at γ ^* = 1.0 is impressive. It is important to note that the agreement between theory and simulation shown in Figure 7.2, is a test of the predictions of the theory for an entirely nonlinear effect. It is a more convincing check on the validity of the TTCF formalism than are the results for the shear stress because there is no underlying linear effect.

The results for the x-y element of the pressure tensor in the three dimensional WCA system are given in Figure 7.3. Again the agreement between the TTCF prediction (T), and the Direct simulation (D), is excellent. We also show the long time steady state stress computed by conventional NEMD (denoted, SS). It is clear that our time limit for the integration of the Transient Time Correlation Functions is sufficient to obtain convergence of the integrals (i.e. to ensure relaxation to the nonequilibrium steady state). We also show the Green-Kubo prediction for the stress (GK). A comparison of the linear and nonlinear responses shows that the intrinsically nonlinear response is only generated at comparatively late times. The response is essentially linear until the stress overshoot time (t^*~0.3). The figure also shows that the total nonlinear response converges far more rapidly than does the linear GK response. The linear GK response has obviously not relaxed to its steady state limiting value at a t^* value of 1.5. This is presumably because of long time tail effects which predict that the linear response relaxes very slowly as t^-1/2, at long times.

In Figure 7.4 we show the corresponding results for shear dilatancy in three dimensions. Again the TTCF predictions are in statistical agreement with the results from direct simulation. We also show the steady state pressure shift obtained using conventional NEMD. Again it is apparent that t^* = 1.5 is sufficient to obtain convergence of the TTCF integral. Although it is not easy to see in the figure, the initial slope of the pressure response is zero.

This contrasts with the initial slope of the shear stress response which is G_∞. This is in agreement with the predictions of the transient time correlation formalism made in §7.3. Figures 7.1, 7.3 clearly show that at short time the stress is controlled by linear response mechanisms. It takes time for the nonlinearities to develop but paradoxically perhaps, convergence to the steady state asymptotic values is ultimately much faster in the nonlinear, large field regime.

Comparing the statistical uncertainties of the transient correlation and direct NEMD results shows that at reduced strain rates of unity conventional NEMD is clearly the most efficient means of establishing the steady state response. For example under precisely the same conditions: after 54 million timesteps the TTCF expression for P_xy is accurate to ± 0.05, but the directly averaged transient response is accurate to ±0.001. Because time is not wasted in establishing the steady state from each of 60,000 time origins, conventional steady state NEMD needs only 120 thousand timesteps to obtain an uncertainty of ± 0.0017. If we assume that errors are inversely proportional to the square root of the run length, then the relative uncertainties for a 54 million timestep run would be ± 0.05, ± 0.001 and 0.00008 for the TTCF, the directly averaged transient response and for conventional NEMD, respectively. Steady state NEMD is about 600 times more accurate than TTCF for the same number of timesteps. On the other hand, the transient correlation method has a computational efficiency which is similar to that of the equilibrium Green-Kubo method. For TTCFs time origins cannot be taken more frequently than the time interval over which the TTCFs are calculated. An advantage of the TTCF formalism is that it models the rheological problem of stress growth(Bird et. al., 1977), not simply steady shear flow, and we can observe the associated effects such as stress overshoot, and the time development of normal stress differences.

Figure 7.5 shows the transient responses for the normal stress differences, P_yy-P_zz and P_xx-P_yy, for the three dimensional WCA system at a reduced strain rate of unity. The normal stress differences are clearly more subtle than either the shear stress or the hydrostatic pressure. Whereas the latter two functions seem to exhibit a simple overshoot before relaxing to their final steady state values, the normal stress differences show two maxima before achieving their steady state values (indicated SS, in the figure). As before it is apparent that t^* = 1.5, is sufficient time for an essentially complete relaxation to the steady state.

Over the years a number of numerical comparisons have been made between the Green-Kubo expressions and the results of NEMD simulations. The work we have just described takes this comparison one step further. It compares NEMD simulation results with the thermostatted, nonlinear generalisation of the Green-Kubo formulae. It provides convincing numerical evidence for the usefulness and correctness of the Transient Time Correlation Function formalism. The TTCF formalism is the natural thermostatted, nonlinear generalisation of the Green-Kubo relations.