J. CHEM. SOC. FARADAY TRANS., 1994, 90(13), 1931-1940 Self-diff usion and Viscoelasticity of Dense Hard-sphere Colloids David M, Heyes and Paul J. Mitchell Department of Chemistry, University of Surrey, Guildford, UK GU2 5XH Brownian dynamics (BD) simulation is used to calculate the viscoelasticity of model near-hard-sphere colloidal liquids using a continuous potential r-“ interaction between the model colloidal particles. The exponent n was varied between 6 and 72. The real and reciprocal components of the complex shear viscosity, q’ and q”, were computed via time-correlation functions under no-shear conditions using a Green-Kubo formula. Also, oscil-latory shear non-equilibrium BD was employed at finite strain amplitude in the linear-response regime. We find that the normalised stress autocorrelation function can be approximated very well by a two-parameter stretched exponential over the complete volume-fraction range.The parameters used to specify the stretched exponential and also the viscosities and long-time self-diffusion coefficients are quite sensitive to the value of n at a chosen volume fraction. The Newtonian viscosity decreases and the long-time self-diffusion coefficient increases with the softness of the interaction, in agreement with experiment. The value n = 36 gives best agreement with the experimental data, and therefore appears to be a good ‘effective’ interaction which we suggest includes the time-averaged effects of the many-body hydrodynamics to some extent. The state dependence of the derived spectrum of relaxation times is determined.As for experimental systems, the complex viscosity scales with the dimensionless (‘ longest’) relaxation time, Doz,/a2, where a is the radius of the particle and Do is the self-diffusion coefficient in the zero-density limit. Also in the intermediate frequency regime 20 <a2w/Do < 200 we find that both the real and imaginary parts of the complex shear vis- cosity decay as ca. o-’l2in agreement with experiment. The rheology of stabilised dispersions has been widely studied by experiment, simulation and statistical mechanics theory. An area of particular interest recently is their visco- elastic behaviour which is conveniently characterised by a complex viscosity composed of real and imaginary parts, q’ the in-phase and q” the out-of-phase response, respectively. The linear viscoelasticity for stabilised dispersions has been measured experimentally by van der Werff et al.’ in the volume-fraction range 0.10 < 4 < 0.60.They discovered an inverse-square root dependence of both q’ -q’(c0) [where q’(c0) is the infinite oscillation frequency viscosity] and q” in an intermediate frequency regime ca. 20 < a20/D, < 200 where a is the radius of the particle and Do is the self- diffusion coefficient in the zero-density limit. Cichocki and Felderhof2 and de Schepper et aL3 derived analytic expres- sions for the complex viscosities which also have a high-frequency ca. o-’/* limiting behaviour. Here we report the results of BD simulations of model sta- bilked dispersions that explore the viscoelastic response of model colloidal liquids consisting of spherical particles, as a function of the volume fraction of solids and the extent of repulsiveness of the interaction potential.This is the first time that this has been performed for a simple generic colloidal interaction using a well established representation of colloidal dynamics (the Langevin equations of motion without many- body hydrodynamics) which is often used as a reference rep- resentation of the colloidal liquid state. Wide ranges of solids fraction and particle softness are considered. The motivation is to discover how this simple model system (with a minimum of parameters and assumptions) behaves, before extending the model in the future to include hydrodynamic forces.On the technical side we propose two new methods to cal- culate the linear dynamic viscosities for model colloidal systems. In the first method, the temporal evolution of the shear stress fluctuations in an unsheared sample are used to calculate a time-correlation function, C,(t), which is numeri- cally identical to the shear stress relaxation function mea- sured in linear stepstrain rheology experiments. The Fourier transform of C,(t) gives the complex moduli and shear vis- cosities. This approach has the advantage of applying no shear rate or strain to the system, so the response function is guaranteed to be in the linear-response regime, which can be time consuming to establish by the direct application of a shear strain profile to the system in experiment or simulation.The second route to the dynamic viscosities is to use non- equilibrium BD. The model dispersion is subjected to a homogeneous oscillating shear flow field with an exp( -iot) time variation, in analogous fashion to the operation of oscil- latory shear rheometers where the liquid is excited into a non-equilibrium state by an externally applied strain field. Brownian Dynamics Simulations The BD method is the same as we have employed in previous work4 and calculates the trajectories of N model colloidal particles from N coupled irreversible equations of motion. If each colloidal particle has mass m, index i and is at position denoted by ri within the (cubic) BD cell, then the position update algorithm is based on, mL’i = Fi + Ri -5i.i (1) where F is an effective interaction force calculated from Fi = -v c V(Jri-rjl) (2) j# i where V is an effective non-bonding chemical interaction between colloidal particles i and j, which is assumed to be pair-wise additive.R is the Langevin random force and 5 is the friction coeficient. The timescale for momentum relax- ation of the colloidal particle, called the Brownian relaxation time, is zB= mt-’ = m/3nrsqs, where qs is the viscosity of the suspeqding medium. A finite difference integration of eqn. (1) is used to evolve the assembly of particles through time and space, r,(t + At) = r,(t) + [F,(t) + R,(t, At)] x At/( + P(t)Atr, (3) where At is the time-step.The last term in eqn. (3) allows for the inclusion of a time-dependent linear shear-flow field in 1932 the suspending medium. In the present model we have omitted many-body hydrodynamics, to discover the predic- tive ability of a simple 'reference' model of a colloidal system. It is interesting to discover the properties of this (still much used) description of a colloidal system for the volume-fraction and interaction-potential dependence of the complex vis- cosity. The colloidal particles interact through an inverse- power potential, V(r)= E(CT/T)" (4) where o is the equivalent hard-core diameter of the model colloid molecule and r is the separation between the centres of two model particles.We set E = k, T and 6 < n < 72. This interaction would represent a stabilised colloidal particle and for n greater than ca. 36 is sufficiently hard to be equivalent to a hard-sphere system for many purposes. The reduced number density of particles, p = No3/V and the solids volume fraction, 4 = 7rp/6. The frequency, a,is made dimensionless by multiplying by a characteristic structural relaxation time, z,, which is the time it takes a colloidal particle at infinite dilution to diffuse a distance a = 012, Z, = 3TCCT3z7s/4kBT = a2/D, (5) For colloidal particles of diameter in excess of 0.1 pm we have, tB< t,, so we can choose a time-step h such that T, 4 h 4 t,. The time-step in the simulation, h, is chosen with h = 622D,, where 6, is the standard deviation of the random displacement.The value of 6, is chosen as large as possible within the bounds of algorithm stability and accuracy (determined empirically from a trial series of simulations), The values of z, and h depend on the choice of Z, (an input parameter), although the ratio of h/z, (the efficiency of trav- elling through configuration phase space) is independent of the value of zB. We typically chose the value zB = 0.316 x 10-3~(rn/~)1/2and 6, = 0.0090. The value of hz, < in these potential-based reduced units is comparable to the values chosen by other workers, (e.g.ref. 6 and 7). For the inverse-power potentials, the interaction energy, pressure and mechanical properties are simply related.The average interaction energy per particle, u, is given by rN where rij= ri -rj. Then the osmotic pressure is given by, P = np(u)/3 (7) and the infinite-frequency shear rigidity modulus in the zero strain amplitude limit, G, = G'(w .+ m), is given by G, = (n2 -3n)p(u)/15 (8) making use of the formula of Zwanzig and Mountain8 Simi- larly for the infinite-frequency bulk rigidity modulus in the zero-strain amplitude limit, K, = K'(o-, co),is given by, K, = (n2 + 3n)p(u)/9 (9) The components of the stress tensor, CT,are needed to compute the dynamic rheology uaS= (raijrpij/rij)Vij (10)N i=l j=i+' Eqn. (10)only includes the contribution to the stress from the direct interactions between the colloidal particles. It does not include the solventsolloidal particle or solvent-mediated hydrodynamic forces between the colloidal particles.The stress in this free-draining level of approximation is treated at J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 a pair-wise additive level. For a more realistic representation including hydrodynamic interactions, the relationship between particle stress and macroscopic properties is much more complex (see e.g. ref. 9). In eqn. (lo), the total stress has a kinetic term pk, TI, where I is the identity matrix. On the timescale of the motion of the Brownian particle, they have reached thermal equilibrium owing to the large number of collisions with the solvent molecules. The magnitude of the kinetic component of the stress is negligible compared with the potential term [eqn.(lo)]. Therefore we do not need to consider this component. Also, as we are mainly interested in the off-diagonal elements of the stress tensor, we do not need to include it because the kinetic contribution to the stress is a diagonal matrix on the rheological timescale. Self-diffusion We also compute the self-diffusion coefficient from the 'local' slope of the mean-square displacement, W(t). For an arbi- trary particle, index 1, we have, W)= ml(t) -r1(0)l2> which is averaged over all particles and dWD(t) = -dt The short-time diffusion coefficient, Ds is D(t)as D(t + 0) and the long-time diffusion coefficient, DL is D(t + m). Let us define the relaxation function, CD(t), D(t) -DL cD(t) = Ds -DL (14) It follows from the Smoluchowski equation that the dimen- sionless function CD(t)can be expressed as CD(t)= rP,(u)exp( -ut) du (15) with a spectral density normalised to, [PD(u) du = rPD(u) d ln(u) = 1 (16) A two-pole approximation to the spectral density valid as p --+ 0 derived by Cichocki and Felderhof" if we neglect hydrodynamic interactions is (17) We make a comparison between this analytic [i.e.eqn. (17) substituted in eqn. (15)] and the BD simulated CD(t)below. Shear-stress Time Correlation Function A linear-response expansion of the position Langevin equa- tion used in this work'' leads to a Green-Kubo expression for the linear shear viscosity in terms of the shear-stress time autocorrelation function, C,(t)defined as where indicates an average over time origins in eqn.(-0.) (18). This method was first used by Levesque et who aplied it to molecular liquids. The infinite-frequency linear shear modulus is given by G, = C,(O). In an unsheared system, the stress fluctuations of all the off-diagonal elements of the stress tensor are equivalent. Therefore, we improved the statistics of C,(t) by considering ox,, ox2,and oy2separa-tely in eqn. (18); and then averaging the three at the end of J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 the simulation. (The stress tensor is symmetric so gap= osol.) The function, C,(t),is exactly the same function as the stress- relaxation function derived from step-in-strain experiments taken in the linear strain limit.The time correlation functions have to extend typically for ca. 20000 time-steps to ensure decay of the function to zero. In order to reduce the com- puter memory requirements, the correlation function was constructed in a piece-wise fashion from three separate corre- lation functions with time origins started (and with a resolution of) every 1, 10 and 100 time-steps. These three correlation functions extended for progressively longer in time, and non-overlapping pieces were merged for the purpose of subsequent analysis and presentation. The number of entries in the histogram used to calculate the time correlation function decreases as time increases. Nevertheless the statistics are reasonable for the ca.500000 time-steps covered for each production simulation. The present BD model’ only incorporates the thermody- namic interactions between the colloidal particles and ignores the many-body hydrodynamic solvent-mediated forces. A measure of the contribution of many-body hydrodynamics to the total viscosity is given by the experimental value of the viscosity in the second Newtonian plateau, qm, which is entirely hydrodynamic in origin. If we assume that the hydro- dynamic contribution to the viscosity is equal to q, at all shear rates, then the Green-Kubo formula in the present model gives the difference between the Newtonian viscosity, qo (the zero-shear-rate limit) and qm. qo is related to C,(r) through vo = dt (19)v, + p) where we need to take a value for q,,, from another source.This is not as serious as it may first appear because often the main interest, certainly for dynamic rheology, is the deviation in the behaviour of the sample from the qocvalue. To be rig- orous, adding the high-frequency or high Peclet limiting vis- cosity in this fashion should be accompanied by a correction to the Stokes drag coefficient as the short-time diffusion coef- ficent is strongly affected by hydrodynamic interactions.’ It is convenient to render the colloid liquid’s viscosity in dimen- sionless units, by dividing them through by the host liquid’s viscosity, to form the relative viscosity, qro = qo/qsand qrL = q,/qs. The complex dynamic viscosity, q*(u)= ((0)+ iq”(o) (20) is related to the dynamic shear modulus, G*(o) through G*(o) = wq*(w).We have = q‘(c0) + FI i C,(t)exp(-icot) dt 1 (21)~’(0) r and q”(w)= 9I i rCs(t)exp(-iot) dt I (22) where q‘(m)= qm which is hydrodynamic in origin.Cs(t)can be represented by a superposition of exponential relaxation functions, C,(t*) = G, J:’P,(~)exp( -t*/z) dt (23) with the normalisation condition j2 P,(z) = 1. Then ~*(c;o)= f(m) + G, where H(z)= zP,(z). Cichocki and Felderhof” give a limiting form for P, as p +0 in the absence of hydrodynamic inter- actions, 1 9U3l2 n (1 -4~)~Ps(u)= -+ 9u(l -u)’ where t = 2t*U2/9D0. If we substitute eqn. (25) in eqn. (23) we have an analytic prediction for C,(t). Applied Oscillatory Shear The second method we use to calculate the viscoelasticity is to apply an oscillating shear strain field to the colloidal par- ticles.The contents of the BD cell are sheared homoge-neously with a time-dependent strain, y(t) over n cycles, y(t) = ’Jo cos ot (26) where yo is the strain amplitude. The analytic expressions for the dynamic viscosities are q’ = -i r2naIw a,,(t‘) dt’ cos(ot’) nT0‘J 0 and ql = -a,,(t’) dt’ sin(wt’) n=yo‘J0 The number of cycles in a simulation of a specified number of time-steps and oscillatory frequency was adjusted to be a whole number by scaling down the time-step. The number of cycles varied with the value of oz,.The number was adjusted empirically to produce reasonable statistics for the average quantities. For example, for COT, > 100 we chose n = 200 and for oz,< 0.1 we used n = 10.Limitations on the availability of computer time necessitated a smaller number of cycles at low frequency, as each cycle takes an increasing number of time-steps as the frequency decreases. Therefore, the sta-tistical uncertainty of the cycle averages was larger as the fre- quency was lowered. In the following section we use the time-correlation method and direct application of a shear flow to explore the linear viscoelasicity of these model colloidal liquids over a wide volume-fraction range. We focus on the effect of the softness of the interaction potential, which is adjusted via the value of the potential exponent, n, in eqn. (4).The softness in experimental systems can be adjusted by varying the relative diameter of the (hard) core and the stabilising barrier (shell).13 In order to cover as much of parameter space as possible, we investigated 4 over essentially the whole liquid volume fraction range for n = 36.Also at the higher 4 values we explored the effects of varying the value of n between n = 6 and n = 72. Results and Discussion A series of equilibrium BD simulations was carried out at a range of volume fractions using N = 108, 256 and 500 model colloidal particles in the simulation cell. Simulation details and the computed properties of the system are given in Table 1. Above a volume fraction of 4 x 0.4 there is an increasing system-size dependence, as expected because configurational phase space becomes more structured with increasing 4.Only those configurations consistent with a finite periodic system are allowed, which produces an increasingly more unrepresentative average with increasing 4. Generally the internal energy (and other derived properties) and the vis- cosities decrease with increasing N. The self-diffusion coefi- cients increase with N. J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 Table 1 Thermodynamic, mechanical and Newtonian transport coefficients by Green-Kubo n 4 tsiml'r u KD ~ 108 36 0.075 455 0.0322 27.98 0.880 0.365 0.016 0.038 108 36 0.115 143 0.0563 48.97 0.799 0.979 0.044 0.073 108 36 0.150 444 0.08 19 71.35 0.749 1.857 0.069 0.115 108 36 0.250 113 0.1892 166.3 0.581 7.155 0.295 0.366 108 36 0.350 63 0.3893 346.3 0.402 20.6 1 0.995 1.17 256 36 0.400 330 0.5443 485.7 0.29 1 32.93 1.90 2.26 108 18 0.427 195 1.3359 272.5 0.279 19.61 1.90 3.34 -108 36 0.427 397 0.6572 590.1 0.253 42.45 3.05 -256 36 0.427 377 0.6574 589.9 0.235 42.46 2.40 -500 36 0.427 170 0.6599 592.3 0.230 42.62 2.95 500 6 0.450 74 3.7603 71.6 0.455 3.88 0.54 4.79 -500 12 0.450 180 2.1908 192.9 0.285 13.56 2.16 -500 18 0.450 202 1.5434 3 17.3 0.227 23.88 2.75 -500 36 0.450 207 0.7750 700.0 0.193 52.75 3.97 -500 72 0.450 1228 0.37 18 1323 0.177 105.83 >2.70 500 6 0.472 159 4.0706 78.84 0.43 1 4.40 0.61 7.00 -500 12 0.472 175 2.4609 218.91 0.254 15.97 2.54 -108 18 0.472 60 1.7518 362.33 0.195 28.43 3.40 -500 18 0.472 164 1.7631 365.00 0.196 28.61 3.6 -108 36 0.472 153 0.8939 81 1.7 0.154 63.82 7.10 -256 36 0.472 151 0.8945 811.8 0.1 50 63.86 7.05 500 36 0.472 217 0.90 16 819.1 0.157 64.37 5.7 --500 72 0.472 46 0.4392 1567.1 <0.136 131.12 >6.69 ~___ KD is the value of qro -q,, from the Krieger-Dougherty form~lae.'~ The mean-square force in the a-Cartesian direction is (F:).The thermodynamic and mechanical quantities are accurate to ca. 1%, while the viscosity and self-diffusion coefficients have a larger error which, at the highest solids fractions, is ca. 10%. Only the interaction component of the thermodynamic and mechanical quantities are given. tsim is the length of the production simulation. We first consider the variation of the self-diffusion coeffi- increases, the particles do no have as far to diffuse before they cient with volume fraction and interaction softness.In Fig. 1 interact strongly with their first coordination shell, which is shown the volume-fraction dependence of the diffusion gives rise to the start of the 'long-time' diffusion regime. In coefficient relaxation function, CD(t)for the n = 36 potential. Fig. 2 we explore the effect of softness of the interactions on The function decays more rapidly with time with increasing the form of CD(t),for the q5 = 0.450 state and varying volume fraction. At the lowest volume fraction considered, 6 d n d 36. Although broadly similar, there is an indication 4 = 0.075, there is very good agreement between the pre- of a more rapid decay with time for the higher n.The Iong- diction of eqn. (15) and the simulated curve up to times time self-diffusion coefficient becomes larger with increasing <0.2a2/D, at least. The agreement is not so good at longer softness of the potential (decreasing n) as shown in Fig. 3 for times, although it must be borne in mind that there is greater 4 = 0.472 states. The soft potential facilitates cooperative statistical uncertainty in the simulation results at low volume motion of the particles and consequently enhances the long- fractions. The simulated curve is a little above the theoretical time self-diffusion coefficient. The slower decay of the softer prediction, which we note is itself an approximation. The (smaller n) potentials especially for r > 0 leads to greater simulation CD(t)decays more rapidly with time as the solids coupling between the trajectories of the particles and hence a fraction increases, which we attribute to the decreasing dis- larger mobility.The particle self-diffusion coefficient tance between nearest neighbours. As number density decreases as the particles aproach closer to the hard-sphere limit (i.e.as n + co). 1.4 1.2 1.o 1.0 c 0.8 n 0.8Bcv 0 0.6 c, 0.4 0.2 0.0 -0.2 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 tD,a-2 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7Fig. 1 CD(t)defined in eqn. (14) as a function of volume fraction, C$ tD,a-2for n = 36. 4, from top to bottom, 0.075, 0.115, 0.150, 0.250, 0.350, 0.427, 0.450, 0.472, 0.490 and 0.527. (-) is the C$ -0 prediction Fig.2 As for Fig. 1, except that CD(t) as a function of n for given by eqn. (15). 4 = 0.450 are shown. n, from top to bottom, 6, 12, 18 and 36 J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 1.1 1 I I I I I 1.0 c 4 0.9 0.8 0.7 5c 0.6 a‘ 0.5 0.4 0.3 0.2 ’ 0.1 4 1 I I I I 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 tD,ir2 Fig. 3 D(t)/D, for 4 = 0.472 for n, from top to bottom, 6, 12, 18, 36 and 72 We now turn to the stress relaxation function C,(t),which unlike CD(t)represents a property of the whole system rather than being a single-particle function. Fig. 4 shows the C,(t) for n = 36 at a range of volume fractions. The function becomes more long-lived with increasing 4, in contrast to cD(t)which decays more rapidly with increasing volume frac- tion.Therefore there is a qualitatively different density depen- dence between cD(t) and C,(t).In general, and particularly at low 4, CJt)decays significantly more rapidly with time than the corresponding CD(t).This is reasonable as the shear-stress pair interaction decays rapidly with pair separation. There- fore the particles do not have to move as far to cause a decay of the same magnitude for C,(t) when compared with cD(t). The function C&) essentially measures the extent of displace- ment of the particles, whereas C,(t) also measures this dis- placement but, in some sense, multiplied by the rapidly decaying pair-interaction function, V(r). The displacement trend inherent in C,(t) is ‘scaled’ by the interaction potential to form C,(t).The great difference in timescales between cD(t)and C,(t), produced by the simulations at low volume fractions, is not reflected in the corresponding analytic approximations to these two decay functions based on eqn. (15) and eqn. (23, respectively. The distinct difference in timescales does not 1.0 I I 1 I I I 0.8 i\ 0.6 h ;c 0.4 Ll 0.2 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 tD,a-2 Fig. 4 C,(t) for 4, from top to bottom, 0.075, 0.115, 0.150, 0.250, 0.350, 0.427 and 0.500. N = 108 in each case. (-) is the 4 -+ 0 pre- diction given by eqn. (25) used in eqn. (23). 1935 appear to be accounted for in the analytic approximations. The analytic model is essentially a zero 4 expansion.Using eqn. (25), the analytic C,(t) decays significantly more slowly than the lowest volume fraction BD computed function (for 4 = 0.075). The discrepancy is somewhat reduced as volume fraction increases but the volume fraction is not a feature of the theory. The C,(t) can be represented very well by a so-called ‘fractional’ or ‘stretched’ exponential [except at very short times t/z, < 0.2 x where this analytic form underesti- mates the simulation C,(t)]. The stretched exponential has the analytic form Two examples of least-squares fits to the BD C,(t) using the analytic form of eqn. (29) are given in Fig. 5 for two extreme values of n at 4 = 0.472. The figure shows the BD and fitted C,(t) for n = 6, the softest potential considered and n = 72, the most steeply repulsive interaction used.Bearing in mind that there are only two disposable parameters the level of agreement between the fitted and BD function is very good in both cases. This stretched exponential function is useful for post-simulation analysis of the viscoelastic behaviour of the system, because of its comparative analytic simplicity. In Table 2 we present the values of z’ and for the different systems. We also give in the table the mean relaxation time z = z‘T(l/p)/p for the normalised function The simulations show that the mean relaxation time increases with volume fraction in agreement with the experimental data on sterically stabilised silica particles by van der Werff et a/.’ In Fig.6 we show the function C,(t) for a range of softness parameters (n) at 4 = 0.450. The function decays more rapidly as n increases i.e. the potential becomes more repulsive. The interaction becomes more short-ranged with increasing n so that the particles do not have to move apart as far before their contribution to the stress is significantly reduced. This causes C,(t) to decay more rapidly for the high n, steeper potentials. Both C,(t) (Fig. 2) and C,(t) (Fig. 6) manifest the same qualitative variation with n. In the 4 +0 limit the stress relaxation function C,(t) appears to be relatively insensitive to the value of 4 (Fig. 4 1.o 0.8 0.6 h%0.4 c, 0.2 0.0 I I I I I-0.2 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 tD,a-2 Fig.5 Normalised stress relaxation autocorrelation function from simulation, (-) and (---), and the stretched exponential least- squares fit, (-----) and (. . .), for the states (a) n = 6, 4 = 0.472, = 0.653 and T’ = 0.0414, and (b) n = 72, 4 = 0.472, 0= 0.238 and Z‘ = 0.0009066. N = 500. 1936 Table 2 Parameters of the stretched exponential least-squares fit to the simulation C,(t)from eqn. (29) 256 36 0.400 0.331 0.00372 0.0228 108 18 0.427 0.45 1 0.01 7 194 0.0423 108 36 0.427 0.330 0.00420 0.0292 500 6 0.450 0.632 0.0427 0.0602 500 12 0.450 0.473 0.0258 0.0578 500 18 0.450 0.433 0.0186 0.0502 500 36 0.450 0.315 0.00465 0.0341 500 72 0.450 0.29 1 0.000830 0.00841 500 6 0.472 0.653 0.04 14 0.0563 500 12 0.472 0.485 0.0306 0.0648 108 18 0.472 0.425 0.0195 0.0553 500 18 0.472 0.416 0.0191 0.0568 108 36 0.472 0.287 0.00369 0.0398 256 36 0.472 0.283 O.Oo409 0.0465 500 36 0.472 0.309 0.00472 0.0373 500 72 0.472 0.238 0.0009 1 0.0198 The mean relaxation time is defined in eqn. (30).and Table 2). In this region the parameters for the stretched exponential representation of C,(t) have values, = 0.360 & 0.005 and z’ = 0.0035 f0.001 for 4 < 0.25 and n = 36. There is a gradual decrease in the value of p and increase in z’ with increasing volume fraction. The functional form of eqn. (29) can be written alternatively in the relaxation time spectrum form of eqn. (23). The distribution of relax- ation times Ps(~)or equivalently H(z)for several representa- tive choices of p and z’ are presented in Fig.7 using the numerical procedure developed in ref. 15 specifically for the stretched exponential. The figure shows that the stretched exponential produces a relaxation-time distribution function which is broad at the low z end but has a relatively sharp cutoff at the high z end, which becomes sharper as j?+ 1. (At fl = 1 we recover a process with a single relaxation time, so the relaxation spectrum takes the form of a S function.) There is a slower decay of H(z) in the z +0 limit. The spectrum of relaxation times becomes broader as /I decreases with increasing 4. If each of these relaxation times is identified with a particular microscopic mechanism for stress release, then this trend indicates that there are many more significant dynamical processes that come into existence at the higher volume fractions.The growth in population of relaxation times is especially at the lower z end. The rapid decrease in C,(t)occurs in the time domain of the transition between the short-time and long-time diffusion coefficients. 1.0 I I 1 I I I I 0.8 0.6 hc % 0.4 u 0.2 0.0 -0.2 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 rD,a-2 Fig. 6 Normalised stress relaxation autocorrelation function from simulation for n, from top to bottom, 6, 12, 18 and 36 at 4 = 0.450 and N = 500 J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 -1 n h c) v3 -7 -6 -5 -4 -3 -2 -1 log(7&la2) Fig.7 Stretched exponential relaxation time spectrum distribution function for the 4 = 0.450 state values: (a) = 0.291 and z‘ = 0.000840, (b) p = 0.473 and 7’= 0.00258 and (c) fi = 0.632 and Z’= 0.04266 In Table 1 we give the integrated shear viscosity from eqn. (19). The states at # = 0.472 and above manifest a strong N dependence with systems in excess of N = 500 considered essential. The # = 0.472 state for example has an N depen-dence with the N = 500 state point qro -qrm value being 5.7 as opposed to 7.0 for N = 256. These high volume fractions and relatively small system sizes have a noticeable N depen-dence for the viscosity. The trend is for the shear viscosity to increase for the smaller N systems. This is reasonable as the trajectories of the particles through phase space are more frustrated by the constraints imposed by the periodic bound- ary conditions.Pathways to lower stressed states have higher activation energies than for larger N systems, and this causes the viscosity to be higher. The relaxation mechanisms are slower, and they have longer stress relaxation times and therefore larger viscosities. The N dependence of the viscosity decreases with volume fraction. For example, Table 1 reveals that for n = 36 the 108 to 500 values of qrO-q,, at $ = 0.427 are 3.05 to 2.95, respectively, compared with 7.10 to 5.7 at 4 = 0.472. Therefore there is a small N dependence at $ = 0.427 but a more significant effect at 4 = 0.472. Over the complete volume-fraction range the following analytic expressions fit the experimental relative viscosity data of near hard-sphere dispersions (summarised in the Krieger-Dougherty expression^'^) quite well.The expressions are, and, qroo= (1 -4/0.71)-2 Table 1 shows that the Krieger-Dougherty expression for qrO -qrm, give values that are quite close to the simulation results for n = 36 (obtained by numerical integration of the time-correlation functions). This is noteworthy as the model does not include many-body hydrodynamics in the equations of motion, which is often considered to be important. Clearly these interactions are present in the real system. Nevertheless, these data indicate that at these high volume fractions, the simple Langevin dynamics embodied in the computer algo- rithm used here, coupled with the continuous interaction V(r),in some ‘mean-field’ sense reproduces well the rheology of the real systems.More significantly, perhaps, is that it sug- gests that the contribution of the many-body hydrodynamics to the non-Newtonian viscosity is relatively insensitive to J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 shear rate. At high volume fraction there is probably appre- ciable screening of the hydrodynamic interactions, so that their effective range is likely to be only several particle diam- eters at most. At lower volume fractions, this screening is likely to be less effective, so that the current hydrodynamics- free approach should be a poorer representation of the system and a mean-field hydrodynamics treatment less ade- quate.This is corroborated by the viscosities given in Table 1. Best percentage agreement with experiment using n = 36 is for the high 4 states, and it becomes less satisfactory as #) +0. In Table 1 we see that the interaction energy, mean-square force, elastic moduli and viscosity increase with n. This trend for thermodynamic and mechanical properties is a feature of inverse power potentials, for which it has recently been shown that u cc n-’.16 The behaviour of the viscosity is con- sistent with the experimental evidence of Mewis et all3 who showed that for the same effective hard-sphere particle diam- eter the Newtonian viscosity (and also the so-called infinite shear rate viscosity) increased with the particle ‘hardness ’.These experimental results are consistent with the trends exhibited by our model systems. As n increases, then the New- tonian viscosity increases. The experimental colloid particles had a varying stabilising layer thickness. Larger colloidal particles with the same stabilising shell thickness are harder when distance is scaled by the its mean diameter. Mewis et al. noticed that larger particles at the same mean volume frac- tion exhibited higher relative Newtonian viscosities.’ The simulated viscosities decrease with decreasing n. We now consider the oscillatory shear simulations carried out using the non-equilibrium BD method. We are interested in the linear response region here, so it is important to deter- mine the maximum strain amplitude that can be used (which will depend on frequency) to remain, to a good approx- imation, within the linear response regime.In Fig. 8 we explore the strain amplitude dependence of the storage modulus, G’, for a 4 = 0.472 state. Simulations were carried out over a range of yo between 0.01 and 1.5 at a series of fixed frequencies, ma2/Do= 10100 and 1000 which covers the practical range of interest. Fig. 8 shows that the response is linear for yo < 0.1 at all frequencies, and there is a reduction in the dynamic modulus or shear strain ‘softening’ for higher amplitudes. There does not appear to be a significant fre- quency dependence to the value of yo above which there is a significant departure from linear response.A value of yo = 50 I -kl 30 --0” 25 0+ + + ++++ ++ +20:; 1 15 -+ O0 0 10 t 1 Fig. 8 Strain amplitude, yo, dependence of the storage modulus of an N = 108 and 4 = 0.472 system at the frequencies, oaZ/D,: 0,10;+, 100 and 0,lOOO 1937 0.02 was chosen for the oscillatory shear simulations reported below to ensure linear response. It is interesting to explore the time-average response of a range of physical properties to the applied oscillatory shear as the strain amplitude increase into the non-linear response region. Fig. 9 shows the amplitude dependence of the average interaction energy per particle, u, for the 4 = 0.472 and n = 36 state point. We note that in the strain amplitude regime where there is shear softening of G’ (ie.yo > 0.1) there is an associated dramatic increase in u, which becomes larger in magnitude with higher frequency. This is reasonable, as an increase in strain amplitude and frequency will both indepen- dently cause a departure from the equilibrium state. We expect their combined effect to be roughly additive. The larger yo and o then the further from equilibrium on average is the state, which causes an increase in the internal energy u. In Fig. 1qa) we show the corresponding plot of the first normal pressure difference, P,, -P,, where P = --b. This manifests a quite distinct dependence on yo when compared with u. At the two highest frequencies considered oa2/Do= 10 and 100 there is a negative dip in the function between 0.5 < yo < 1.0 before going positive and rising steeply with yo at yo x 1 (more steeply than u).We note that the region of rapid ascent is delayed to larger amplitudes than in the case of internal energy. The second normal pressure difference P,, -P,, and its dependence with yo at three fixed frequencies is given in Fig. lqb). These data follow very much the same trend as the first normal pressure differences of Fig. lqa). The complex viscosities obtained from the stretched expo- nential fits to the time-correlation functions using eqn. (21) and eqn. (22) are statistically indistinguishable from those obtained using the direct application of an oscillating shear flow field for oa2/Do< 1OOO. At higher frequencies systematic differences appear which arise from the failure at very short times of the fitted function to agree with the simulated C,(t).The frequency domain of oa2/Do> loo0 is, however, not of much practical interest. The real and imaginary parts of the relative complex viscosity for two widely separated volume fractions are given in Fig. 11. The qf and qf’ for n = 36 shift to lower frequencies as 4 increases, in agreement with the experimental results of van der Werff et al.’ In Fig. 12 we present the real and imaginary parts of the relative complex viscosity for two values of n at the same volume fraction, 4 = 0.472. As the interaction becomes softer (i.e. n decreases) the peak in the imaginary part of the viscosity shifts to higher frequency and becomes sharper (associated with an increase in p, shown in Table 2). The corresponding real part also I I3.0 1 2.8 2.6 2.4 2.2 2.0 U 1.8 1.6 1.4 1.2 -2 -1 0 1 log Yo Fig.9 Variation of the time average value of u with strain ampli- tude at frequencies as in Fig. 8 for the 4 = 0.472, n = 36 states 1938 12 10 8 26 I h kk 4 0 o..., 1-2-‘ -2 -1 0 1 log Yo 35 I I I I 0 30 (b)25 t i0 0 0 0 I I-10 ‘ -2 -1 0 1 log Yo Fig. 10 Variation of (a)the time-average value of P, -P,, and (b) P,, -P,, with strain amplitude at frequencies as in Fig. 8 moves to higher frequency. Although we are unaware of any experimental studies that would confirm or deny this trend, it is intuitively reasonable in the same sense as we rationalised the softness dependence of the Newtonian viscosity.The softer particles (but maintaining the same effective hard- 1.o 0.9 > z 0.8b.$ 0.7 0.6 c -Y x 0.5 8E-v-‘e10.150.4 0 E* 0.3 .-> -x 0.2 ?2 0.1 0.0 -3 -2 -1 0 1 23 4 0.00 0.05 0.10 0.15 0.20 log oa2/Do (w a2/D,)-’I2 J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 00 0 0 0 0 0 ~ O+++&+0.4 +, 1 -3 -2 -1 0 1 23 4 log oa2/Do Fig. 12 As for Fig. 11, except we compare two systems with differ- ent n (0, +, n = 6; 17, x, n = 36) at the same solids fraction, 4 = 0.472 sphere diameter, a) possess interactions that facilitate coop- erative motion.Therefore it is necessary to go to a higher frequency to cause the same extent of viscoelastic response as for harder (i.e. higher n) particles. Workstation movies of the simulation cell under shear distortion in the linear response regime do not reveal any restructuring (e.g. long-range order) during oscillatory shear, if the affine shear distortion is extracted. The dynamic viscosities are presumably therefore primarily influenced by the competing relative timescales of the mobility of the particles within their essentially equi- librium structures, and that of the oscillatory shear (ca. 2n/o). Between a20/D, = 20 and 200 Fig. 13 shows that q’ decays as o-”’.Both q’ and q” can be approximated in this regime by the following analytic forms q; = q:, + A’(a*m/Do)-”2 (33) and 7:’ = A”(a2m/Do)-l’* (34) This behaviour is in agreement with the results of van der Werff et al.’ who measured the viscoelastic behaviour of 1 1 1 n 8 0.20 Fig.11 Real (upper curves) and imaginary (lower curves) part of the Fig. 13 Normalised real part of the relative complex shear viscosity relative complex shear viscosity as a function of frequency for the as a function of (aZw/D,)-’/2for n = 36. 0,4 = 0.350, calculated by n = 36 model colloidal systems at the volume fractions: 4 = 0.150 Fourier transformation of the correlation function. +, 4 = 0.350; 0, (0,+), N = 108 and q5 = 0.450 (0,x), N = 500 0.427 and x ,0.472, calculated by the non-equilibrium BD method. J.CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 sterically stabilised silica spheres. It is important to note that at very high frequency, there is always a departure from this relationship. This is essential in practice as a q' decaying as at higher frequency implies a G' that grows as 01/2 which is physically unreasonable as this is equivalent to G, -,co for a -,00. This is understandable, at least for a system without hydrodynamics, because this limit represents an affine deformation of the microstructure by shear, which is impossible for rigid spheres at contact. The dimensionless constants A' and A" in eqn. (33) and eqn. (34), respectively, have values 1.8 & 0.1, 3.6 & 0.2 and 5.5 & 0.5 at the volume fractions # = 0.35, 0.427 and 0.472, respectively, for n = 36.The experimental values at these volume fractions are approximately double in magnitude. For 4 = 0.35, 0.427 and 0.472 these parameters have the values 3.5 & 0.5, 9.0 f1.0 and 30 f5. The slope of the linear regime of the complex viscosity with LO-'/'is sensitive to the softness parameter n. The dependence of the imaginary component of the complex viscosity with n is presented in Fig. 14. As the value of n diminishes the dimensionless constants A' and A" decrease. In this case, at 4 = 0.472 there are A" = 1.3, 3.7 and 7.8 for n = 6, 12 and 72, respectively. In Fig. 15 we show the imagin-ary component of the complex viscosity with 4 for n = 36 8T 0.20 + :-I0 ox 0 8"$ 0.15 0 +A 0 X r-+ O O *exe: 0.10 0.05 10.00 0.20 0.00 0.05 0.10 0.15 (w aZ/D,)-!2 Fig.14 Normalised imaginary part of the relative complex shear viscosity for the model colloidal systems as a function of (u~~/D,,-''~at 4 = 0.472 for n: 0,6; +, 12; 0,18; x, 36 and A, 72, using the correlation function method 0.35 , 1 a 0.30 -a + A 0.25 h 8 0.20 n ++ h v 3 r-U v z .' 0.15 h 0 z +0 E-0.10 X 0.05 0.00 0.00 0.05 0.10 0.15 0.20 (w a2/D,)-'Iz Fig. 15 As for Fig. 14 except the dependence on 4 for n = 36 is presented. 4: 0,0.075; +, 0.115; 0,0.150; x, 0.250; A, 0.350; *, 0.400;0,0.427; +, 0.450; 0,0.472 and x, 0.490. 1939 using the correlation function route. We note that the linear regime for o-'I2 is compressed into a higher frequency band as volume fraction decreases. Van der Werff et a/.also show that the complex viscosity scales with a reduced frequency, oa2zl/Do where tl is the 'longest' relaxation time in their fit to the data. If we use the van der Werff values for a2zl/D, then the present simulation data also scales onto the same curve which is also statistically indistinguishable from their experimental data (see Fig. 16). The particle size and volume fraction therefore can be incorp-orated in this dimensionless scaling parameter. The values of a2z,/Doare approximately an order of magnitude larger (but roughly proportional to) the values of a2z/Do which we obtained from our simulations. The Stokes-Einstein relationship is for this system D,q,(O)/D, = 1.The ratio DLqr(O)/D0is essentially linearly increasing with volume fraction up to 6 z 0.45, before lev-elling off. This ratio also shows an approximately linear dependence with n-l, as shown in Fig. 17 for two of the higher solid fractions. This indicates that as the softness 1.o ' """""""d, ' 1 1 0.9 t >..z 0.8 8 F.p 0.7 W B+0.6 t 0 0.5 0 -fi 0.4 0 0.3 .-5 -t;; 0.2 ?? 0.1 0.0 -3 -2 -1 0 1 2 3 log wrl Fig. 16 Normalised (a) real and (b) imaginary part of the relative complex shear viscosity as a function of frequency for three volume fractions 4: 0, x, 0.150; +, A, 0.427 and 0,*, 0.490 us. m1,where Dot,/u2 = 0.08, 0.32 and 0.45, respectively. The 'bars' on the figure represent the limits of the band of experimental points from van der Werff et al.' a" 3.4--.- h 3.2 - ca'' 3.0 2.8 - 0 2.6 - 2.4 - + + 0 2.2 I ' L 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 1/n Fig. 17 DLq,(0)/Do11s. n for the simulation data. 4: 0,0.450 and +,0.472. increases (i.e. n-increases) the self-diffusion coefficient increases more rapidly than the relative viscosity decreases. Conclusions We have shown that BD computer simulation can be used to investigate the viscoelastic behaviour of model near-hard- sphere colloidal liquids up to volume fractions of ca. 0.50 using two complementary techniques (equilibrium and non- equilibrium). We have explored the influence of solids volume fraction and extent of softness of the interparticle interactions on the self-diffusion and complex viscosity.The linear stress relaxation function (stress autocorrelation function) is approximated well by a stretched exponential for all states considered. The values of the parameters in the stretched exponential depend on solid fraction and softness of the interaction potential between the particles. Despite the formal absence of many-body hydrodynamics, the simple model adopted accounts well for the experimentally observed trends in viscoelastic behaviour as solids fraction and particle soft- ness are varied. Newtonian and complex shear viscosities are also numerically quite close to the experimental values at the same volume fractions if we choose a potential index, n, close to 36.This suggests that, for the future, it could be possible to develop an essentially ‘mean-field’ treatment of many-body hydrodynamics which is suitable for simulation at high solids fraction, and which is not much more computationally inten- sive than the current model (in marked contrast with existing approximate schemes for including many-body hydrody- namics e.g. ref. 17). P.J.M. thanks the SERC and ECC International Ltd for a research fellowship (grant number GR/H80644). Computa-J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 tions were carried out on the CONVEX C3 at the University of London Computer Centre, Personal DEC Stations and a Silicon Graphics Indigo XZ workstation in the Chemistry Department at the University of Surrey. The referees are thanked for helpful comments. References 1 J. C. van der Werff, C. G. de Kruif, C. Blom and J. Mellema, Phys. Rev. A, 1989,39,795. 2 B. Cichocki and B. U. Felderhof, Phys. Rev. A, 1991,43,5405. 3 I. M. de Schepper, H. E. Smorenburg and E. G. D. Cohen, Phys. Rev. Lett., 1993,70, 2178. 4 J. R. Melrose and D. M. Heyes, J. Chem. Phys., 1993,98,5873. 5 D. M. Heyes and J. R. Melrose, J. Non-Newtonian Fluid Mech., 1993,46,1. 6 K. J. Gaylor, I. K. Snook, W. van Megen and R. 0. Watts, Chem. Phys., 1979,43,233. 7 H. Lowen, J-P. Hansen and J-N. ROUX,Phys. Rev. A, 1991, 44, 1169. 8 R. Zwanzig and R. W. Mountain, J. Chem. Phys., 1965,43,4464. 9 J. F. Brady, J. Chem. Phys., 1993,99,567. 10 B. Cichocki and B. U. Felderhof, Physicu A, 1993,98,423. 11 W. Hess and R. Klein, Adv. Phys., 1983,32, 173. 12 D. Levesque, L. Verlet and J. Kurkijarvi, Phys. Rev. A, 1973, 7, 1690. 13 J. Mewis, W. J. Frith, T. A. Strivens and W. B. Russel, AIChE J., 1989,35,415. 14 W. B. Russel, D. A. Saville and W. R. Schowalter, Colloidal Dis-persions, Cambridge University Press, Cambridge, 1989, p. 466. 15 C. P. Lindsey and G. D. Patterson, J. Chem. Phys., 1980, 73, 3348. 16 D. M. Heyes and P. J. Aston, J. Chem. Phys., 1994,100,2149. 17 G. Bossis and J. F. Brady, J. Chem. Phys., 1989,91,1866. Paper 4/00386A; Received 21st January, 1994