Treating interactions with polar and solvent sites alike, the interaction parameter is The correction of eq A5 for the surface contacts per polymer sphere is now the difference of eqs A12 and A13 and or, with eqs A2 and A3, With eqs A9 and A10 reduction of eq A8 yields Equations A5 and A15 now yield for the corrected contact energy of the two-part solution (copolymer sphere and remainder) El E'" kT
--+
-
-x@'[(Nn + nsN2 + U + ufe)pn(l - $)21
+
4;)
(A161
In eq A16 the first term in the brackets represents the FloryHuggins result for the homogeneous polymer solution, eq 3. Therefore, the second term accounts for the extra intramolecular contact energy of the single polymer sphere. Since there is such a contribution for each polymer molecule the corrected contact energy for the copolymer solution becomes E = -x@'[(Nn + ns)$2 + U + uf)pNn(l - $)2] + O ( P ) kT (A171
We add Ef" to E l in eqs to account for the (1 - a)n,q sphere/ solution contacts. Since the total number of contacts between lattice sites is invariant, we must now subtract from El the amount Ecp which is the contact energy of (1 - u)nc sphere sites among themselves and of (1 - a)$ bulk solution sites among themselves. Methods similar to those used in deriving eq A12 now yield
+
Dropping terms of order No,dividing by M = Nn 4,and using = Nn/(Nn + n,) converts eq A17 into eq 5 of the main text, We note that the derivation of the Flory-Huggins expression in eq 3 yields terms proportional to $, similar to those in eq A12. Such terms have been suppressed because they have no thermodynamic significance.
Brownian Dynamics Simulations of an Order-Disorder Transition in Sheared Sterically Stabiked Colloidal Suspensions Angeliki Artemis Rigos* Department of Chemistry, Merrimack College, North Andover, Massachusetts 01845
and Gerald Wilemski* Physical Sciences Inc., 20 New England Business Center, Andover, Massachusetts 01810 (Received: September 24, 1991)
The shear thinning behavior of a sterically stabilized nonaqueous colloidal suspension was investigated using nonequilibrium Brownian dynamics simulationsof systems with 108 and 256 particles. At a volume fraction of 0.4, the suspension is thixotropic: it has a reversible shear thinning transition from a disordered state to an ordered, lamellar state with triangularly packed strings of particles. The time scale for the transition is set by the free particle diffusion constant. For the smaller system, the transition occurs gradually with increasing shear rate. For the larger system, the transition is sharp and discontinuous shear thinning is found.
Introduction Colloidal suspensions are attractive candidates to study for unusual shear rate dependent phenomena. Although they are simpler in structure and dynamics than polymer solutions and melts, they exhibit much of the same rheological behavior. Moreover, structural transitions occur at moderate shear rates due to the large particle sizes and low Reynolds number flows 0022-3654/92/2096-3981$03.00/0
involved. This is in contrast to molecular liauids where extremelv large shear rates, well outside the expehen'tally accessible rang;, would be required toObserve ampab1e effects found in "PUkr simulations using nonequilibrium molecular (NEMD)*'" (1) Erpenbeck, J. J. Phys. Rev. Let?. 1984, 52, 1333. (2) Woodcock, L.V. Phys. Reo. Lerr. 1985.51, 1513; Chem. Phys. Lerr. 1984,iii,455.
0 1992 American Chemical Society
3982 The Journal of Physical Chemistry, Vol. 96, No. 10, 1992 Shear-induced structural changes in concentrated suspensions have been the focus of much recent experimental work.6-’2 Discontinuous shear thinning7 and evidence for incompletely ordered particle “stringsns are just two of the intriguing new findings. While these and other experimental results” are not yet fully understood, they resemble behavior found in recent computer simulations using NEMD1-5s’3and nonequilibrium Brownian dynamics (NEBD).I4-” NEBD simulations may be a useful supplement to theoretical efforts1s-21in understanding the structure and behavior of sheared colloidal suspensions since they can provide unambiguous results for well-defmed model systems. In contrast to NEMD, the NEBD method is inherently designed to simulate the dynamics of particles suspended in a fluid, and it does not need thermostating which sometimes leads to artifacts in NEMD simulation^.^^.^^ Over the past several years, several investigators have used NEBD techniques to explore the behavior of sheared colloidal suspensions. Bossis and BradyZez7have developed to a high degree the treatment of hydrodynamic interactions among colloidal particles in their Stokesian dynamics simulations. They have found that monolayer suspensions of hard particles shear thicken at very high Peclet numbers due to the formation of large hydrodynamically interacting clusters, while at low Peclet numbers shear thinning occurs because of the disappearance of the direct Brownian contribution to the viscosity.27 They have also observed a layered structure in highly concentrated, strongly interacting systems,24but layering was not induced or affected by increasing the shear rate. Heyes14 has compared the rheology of molecular fluids and colloidal suspensions using both NEMD and NEBD in 3-D. With NEBD and both a Lennard-Jones (LJ) ( N = 108 particles) and a repulsive LJ (RLJ) ( N = 32) potential, he observes continuous shear thinning. Although it is difficult to see the direct effect of system size and potential on his results, he proposes that the attractive tail of the LJ potential has no noticeable effect. H e concludes that shear thinning in both the NEMD and NEBD simulations is mainly influenced by the effective hard-sphere diameter which is determined by the short-range repulsive force. Heyes’ more recent NEBD simulation^,'^ again in 3-D with N = 32 and 108 particles using a soft-sphere potential and without hydrodynamic interactions, produce a shear-induced continuous (3)Hess, S. Int. J . Thermophys. 1985,6,657. (4)Heyes, D. M.; Morriss, G. P.; Evans, D. J. J . Chem. Phys. 1985,83, 4760. (5) Heyes, D. M. J. Chem. Soc., Faraday Trans. 2 1986,82, 1365. (6)Ackerson, B. J.; Hayter, J. B.; Clark, N. A.; Cotter, L. J . Chem. Phys. 1986. 84. 2344. (7) Chen, L. B.; Zukoski, C. F. Phys. Reu. Lett. 1990,65,44;J . Chem. Soc., Faraday Trans. 1990,86,2629. (8)Ackerson, B. J.; Pusey, P. N. Phys. Reu. Lett. 1988,61,1033. (9)Johnson, S. J.; de Kruif, C. G.; May, R. P. J . Chem. Phys. 1988,89, 5909. (10)Ackerson, B. J.; van der Werff, J.; de k i f , C. G. Phys. Reu. A 1988, 37,4819. (11) van der Werff, J. C.; Ackerson, B. J.; May, R. P.; de Kruif, C. G. Physica A 1990,165,375. (12)de Kruif, C. G.; van der Werff, J. C.; Johnson, S . J.; May, R. P. Phys. Fluids A 1990,2, 1545. (13)Stevens, M. J.; Robbins, M. 0.;Belak, J. F. Phys. Reu. Lerr. 1991, 66,3004. (14)Heyes, D. M. J . Non-Newtonian Fluid Mech. 1988,27,47. (15)Heyes, D. M. Phys. Lett. A 1988,132,399. (16)Xue,W.; Grest, G. S. Phys. Reu. Lett. 1990,64,419. (17)Wilemski, G. J . Stat. Phys. 1991, 62, 1239. (18)Hess, S. Phys. Reu. A 1980,22,2844. (19)Ronis, D. Phys. Reu. Lett. 1984,52,473;Phys. Reo. A 1984,29,1453; 1986 34, 1472. (20)Koo, H.-M.; Hess, S . Physica 1987,145A. 361. (21)Dhont, J. K. G. Physica 1987,146A,541;J . Fluid Mech. 1989,204, 421. (22)Evans, D. J.; Morriss, G . P. Phys. Rev. Lett. 1986,56, 2172. (23)Hanley, H. J. M.; Morris, G. P.; Welberry, T. R.; Evans, D. J. Physica 1988. 149A,406. (24)Bossis, G.; Brady, J. F. J . Chem. Phys. 1984,80, 5141. (25)Brady, J.; Bossis, G. J . Fluid Mech. 1985,155, 105. (26)Phillips, R. J.; Brady, J. F.; Bossis, G. Phys. Fluids 1988,31,3462. (27)Bossis, G.; Brady, J . F. J . Chem. Phys. 1989,91,1866.
Rigos and Wilemski
1.05
1.02
1.08
rid
1.1
12
1
rld
F g r e 1. Potential used in the Brownian dynamics simulations. For r / d < 1.007, the steep attractive part of the potential is not shown. The potential maximum is 154kT at r / d = 1.0055.
phase change of the dispersed particles into strings forming a distorted triangular lattice in the plane perpendicular to the streaming direction. Recently, using NEBD, Xue and GrestI6 simulated 3-D aqueous suspensions of polyballs subjected to oscillatory shear flow without hydrodynamic interactions. For concentrated systems of 500 and 504 particles near the melting point, they observed shear-induced transitions to a string-ordered phase with the strings forming a triangular lattice in the plane perpendicular to the flow. They also found that the diffusion coefficients perpendicular to the flow direction sharply declined a t the transition shear rate, which suggests that their systems should also exhibit discontinuous shear thinning. However, they did not explicitly confirm discontinuous shear thinning by calculating the suspension stress versus shear rate. Finally, in Wilemski’sl’ previous NEBD study of 3-D aqueous and nonaqueous suspensions of 108 particles without hydrodynamic interactions, the effects of the colloidal interparticle potential OR shear viscosity were explored at particle volume fractions 4 of 0.2 and 0.4. The aqueous systems, which had solidlike order at rest, were found to exhibit apparent discontinuous shear thinning accompanied by a structural transition from a disordered to a lamellar state. More recent results by Cook and WilemskiZBfor a larger, 256-particle system indicate that these transitions are an artifact of system size, The ‘nonaqueous suspensions were liquidlike and weakly flocculated at rest, and the more concentrated suspensions (4 = 0.4) underwent an order-disorder transition with increasing shear rate. Shear thinning was accompanied by the disruption of the flocs and by the reorganization of the particles into a triangular lattice of strings similar, but not identical, to that found by Heyes.” In contrast to Heyes’ proposalI4 for the limited role of the attractive potential, Wilemski also found that the attractive van der Waals potential had a significant effect on the calculated viscosity at lower shear rates. In this paper, we explore in more detail the order-disorder transition found earlierI7 for the nonaqueous system with attractive forces using 108 and 256 particles. Unlike the behavior of the aqueous s y ~ t e m , ~the ~ , order-disorder ~* transition is found to occur in both 108- and 256-particle systems. In the larger system, the transition to the string-ordered phase is sharp, and the strings are packed in a triangular lattice resembling that found by Xue and Grest.I6 The viscosity and stress versus shear rate curves exhibit (28)Cook, R.; Wilemski, G . J . Phys. Chem., this issue.
The Journal of Physical Chemistry, Vol. 96, No. 10, 1992 3983
Order-Disorder Transition in Colloidal Suspensions
TABLE I: Summary of NEBD Simulations for N = 108 and 256 Particles (f Is the Shear Rate, Step, and Aa/qa the Average Excess Viscosity) N = 108 9 , s-1 10 20 50 100 150 200 36 250 200 330 n/ 1000 40 160 5 5 5 5 80115 5 At X lo6,s A1110 4.14 2.77 1.48 0.986' 0.751,"0.703b 0.448'
9 , s-'
10 n/1000 60 At X 106,s 5 A1190 5.24
20 230 5 3.45d
50 40
100 240
5
5
1.82 1.11'
N = 256 130 135 140 100 100 160 5 80113 80/13.5 4017 0.9309 0.555 0.528 0.5159 125 160
150 200 80115 0.500,"J 0.509b"
R
the Number of Steps per Run, A t the Time
250
300 120 1613 0.211
88 5 0.2729
200 160 5 0.329g
250 128 5 0.274"
300 70 1613 0.243c
400 100 5 0.135 400 60 5 0.193f
500 100
5 0.101
500 96 5 0.1748
750 90 200175 0.1468
"Starting from 9 = 100 s-l. *Starting from = 200 s-l. 'Averaged over last 200000 time steps. "Adveraged over last 80000 time steps, 'Averaged over last 30000 time steps. /Averaged over last 50000 time steps. 8Averaged over last 40000 time steps.
discontinuous shear thinning similar to that found experimentally by Chen and Zukoski.' We believe that this is the first report of an NEBD simulation of discontinuous shear thinning in steady shear flow, although the phenomenon has been found in previous NEMD simulation^.^^^^^^'^^'^ In the transition region, we also observe thixotropic behavi~?~ (timedependent stress and viscosity) with a time constant roughly equal to the free particle diffusion time. Thixotropic behavior has not been reported for previous NEBD simulations, although it has been noted in NEMD simulations.2
Simulation Methodology The simulations were performed using the Brownian dynamics computer code of Wilemskil' which does not include hydrodynamic interactions. Systems of 32, 108, and 256 particles of diameter d = 20 = 1.2 pm were simulated at a volume fraction 4 = 0.4. The interparticle potential was a sum of attractive van der Waals force3O and polymeric steric repulsion forces3' as a function of r, the distance between particle centers
kTS(2Ro + r)(Ro- r)* (1) Ro = d
+ 26
J) is the Hamaker constant, k is the where A (=5.0 X Boltzmann constant, T (=293 K) is the absolute temperature, and 6 (=0.02 pm) is the thickness of the adsorbed polymer layer. S (=5.82 X 10l6~ m - is ~ a) composite parameter characterizing the steric repulsion. The potential is shown in Figure 1. The particles' positions are updated according to the equation Ari = (ui + DE;/kT)At
+ Ri
(2)
where ri is the position vector of the ith particle, D is the selfdiffusion coefficient of the particle (evaluated with the StokesEinstein expression), Fi is the total force on the particle, ui is the undisturbed fluid velocity at the location of the particle, and Ri is the random displacement of the particle due to Brownian motion in the time interval At. The simulations were run using simple shear flow in the x direction with the velocity gradient in t h e y direction, that is, ui = +y,e,, where is the shear rate (s-I), yi is t h e y coordinate of the ith particle, and e, is a Cartesian unit vector. In order to simulate shearing in an infinite system, periodic boundary conditions and coordinate space homogeneous shear conditions24were used. The contribution of the interparticle forces to the stress uxyis given by
+
(3) (29) Bauer, W. H.; Collins, E. A. In Rheology, Theory and Applications; Eirich, F. R.,Ed.; Academic: New York, 1967; Vol. 4,p 423. (30)Hamaker, H.C. Physica 1937, 4, 1058. (31) Ottewill, R. H.; Walker, T. Kolloid 2.2.Polym. 1968, 227, 108.
4
I-
1
0.1
I 10
a 0
I
I
I
l
l
100
1
1 - 1
I
I 1000
Shear Rate, 7 T
Figure 2. Shear rate dependence of the excess viscosity. The leastsquares fits to the equation AT = 4-p give p for the indicated ranges of shear rates. N = 108: 0.64(10-150), 1.68 (150-500),1.45 (250-500). N = 256: 0.69 (10-125),0.62 (200-750),0.41 (400-750).
where yij = yi - y,, FXiiis the x component of the force on particle i due to its interaction with particle j , V is the system volume, and the angle brackets indicate a time average. The corresponding exviscosity of the suspension is Aq = u,,/+. The systems were studied over shear rates of 10-750 s-I, numerically equal to the dimensionless shear rate +T (or Peclet number) = a2+/D = 6?rq,,a3+/kTby the specific choice of a, T, and vo Pa s), the viscosity of the pure suspending fluid. The time scale is set by particle diffusion, T = a2/D, where T conveniently equals 1 s for our system parameters. Using time steps between 2 X 10% and 7 X lo+%, runs of 40000-200000 time steps were executed to ensure that the system was not in a metastable state such as was mentioned in ref 17. Generally, the final particle coordinates from calculations a t one shear rate were used as initial particle coordinates for a slightly different shear rate. Results and Discussion The simulation results for the excess viscosity, Aq, in Figure 2 indicate that this nonaqueous suspension undergoes shear thinning: its viscosity decreases with increasing shear rate. The power law dependence of Av on .i/ (Aq 0: +-p) is similar for the two system sizes at low + but differs substantially at high 9 . Values of p are noted in the caption to Figure 2. Conditions for each run are tabulated in Table I. The most interesting feature of Figure 2 is the shear thinning transition that occurs discontinuously at a critical shear rate of + c ~= 125 for N = 256 particles and continuously for N = 108 = 150. The system size differences in particles beginning at .i/c~ the transition are emphasized in the stressshear rate plots of Figure 3. Below the transition, the suspensions are disordered and liquidlike, and the stress grows linearly with In ( + T ) . Above
3984
Rigos and Wilemski
The Journal of Physical Chemistry, Vol. 96, No. 10, 1992 125,
50
__
1
1
I
1
0
O 0 100
10
150
(b)fT-
V
1000
Shear Rate, i.r
z
Figure 3. Shear rate dependence of the stress.
(d); . . . T = 500
... N-256
1
0
m
0
B
0
a Y
0
0
0
e
W
I
I
0.4 0.6 Time (117)
I
0.8
B B
B
a
rr
1
I
0
B
m
0.2
S
B
B
O 0
80 I 0.0
0
m
Q
0
I
I
1.0
(a) f~ = 20
Figure 4. Time dependence of the stress for N = 256 at shear rates of 100, 125, and 150. The fluctuation plot is generated by dividing the run into intervals of l/(TAt) steps and connecting the average values for successive intervals with straight-line segments.
the critical shear rate, the stress drops abruptly for N = 256 due to the reorganization of the particles into an ordered lamellar structure with lines or "strings" of particles moving easily past each other. In the smaller, 108-particle system, the transition is postponed to a slightly higher shear rate, and the reorganization occurs progressively, not abruptly, with increasing shear rate. This leads to the continuous decrease in the stress seen in Figure 3 as the shear rate is increased past the critical value. What is not evident in Figures 2 or 3 is the time-dependent character of the suspension's response to changes in shear rate. Figures 2 and 3 show only the long-time, steady-state values. However, the suspension is t h i x o t r ~ p i c ,and ~ ~ considerable time must be allowed for the suspension's structure to fully relax after each change in the shear rate. This is particularly important for large changes in shear rate and for changes that traverse the critical shear rate in either direction. This behavior is shown clearly in the stress fluctuation plots in Figure 4 for TT = 150 and 125. These curves reflect the time evolution of the suspension from a disordered to an ordered state (47 = 150), and vice versa (TT = 125). They contrast sharply with the stationary behavior at j 1 =~ 100 which involves no major structural reorganization of the suspension. The details of the time dependence will be easier to explain after we discuss the structural changes that underlie the transition. The nature of the transition is evident from snapshots of particle configurations above and below the critical shear rate. These are projections of the centers of all the particles onto the yz plane, which is perpendicular to the direction of the shear flow ( x direction) and which contains the velocity gradient 0, direction). These nonaqueous systems have liquidlike order at equilibrium, and this disordered state (Figure 5a (N = 108) and Figure 6a (N = 256)) persists up to the critical shear rate (Figure 5b and Figure 6b). At TT = 130,just above the critical value, the larger system is well-ordered with the particles organized into an array of strings (Figure 6c) packed on a triangular lattice. In contrast, at 4s = 200, well above the critical shear rate, the smaller system is only partially reorganized into the string array (Figure 5c). This organization improves with shear rate and reaches nearly perfect
(b)fT=125
Y
(d)fr=750 "
m
0
E
r B
d
r a
P E
0
"
o
a
s Y
0
O
o
m a
O
m
~ B
I @
~
W ~
O
m
m %
O
o
o
I
O
W
m
n
o o I B @ a n P o o w o J o
z Z Figure 6. Snapshots of particle configurations for N = 256 projected onto the plane (yz)normal to the flow ( x direction) at the indicated shear rates. The velocity gradient is in the y direction.
ordering at j~ = 500 (Figure 5d). This gradual increase in ordering correlates with the continuous decline in the stress seen in Figure 3 above the critical shear rate. In the larger system, the degree of ordering is already high just above the critical shear rate and the improvement at high shear rates, e.g. +T = 750 (Figure 6d):is less dramatic. This accounts for the differences in the behavior of the stress in two systems that were noted earlier. The quantitative difference in the stress at shear rates >300 probably occurs because the two-dimensional b z ) packing density of the strings is about 5% higher for N = 256 than for N = 108, so particles in the larger system may be interacting a bit more strongly on average than in the smaller system. The structural change into strings is also evident in yx snapshots at +T = 125, 130,and 750 (Figure 7, a, b, and c) for N = 256.
The Journal of Physical Chemistry, Vol. 96, No. 10, 1992 3985
Order-Disorder Transition in Colloidal Suspensions (a)?T = 125
( b ) f r = 130
I
I
$
5
I 10
,
,
/
,
,
100
,
,
,,I 1000
Shear Rate , f r
Figure 9. Shear rate dependence of the number of nearest neighbors for N = 108 and N = 256.
X
X
Figure 7. Snapshots of particle configurations for N = 256. (a), (b), and (c) are projections onto the yx plane at the indicated shear rates; (d) is
one layer of particles with -2.3 < y < -2.1 projected onto the xz plane for a shear rate of 750. The flow is along the x axis, and the velocity gradient is in the y direction. 14
12 10
-c m
N = 256
-
y t = 100
r
8 -
6 4 -
2 0.
-
I
.
+
fr=750
12
6 64 -
02 0.9
1 1.3
1.7
2.1 r/d
= 1.79, 2.06, and 2.68 match almost perfectly the second, third, and fourth nearest-neighbor distances for the planar hexagonal packing of spheres with effective diameters equal to the first nearest-neighbor distance, r / d = 1.03. This packing typifies the layers found in the ordered suspensions (Figure 7d). It is also instructive to examine the behavior of the number of nearest neighbors with shear rate, since only nearest-neighbor interactions contribute significantly to the average stress. The average number of nearest neighbors was determined by counting only pairs of particles separated by r / d C 1.2; beyond this distance, the attractive potential has effectively decayed to zero. Figure 9 shows that the two systems behave similarly below +c. At these low shear rates, the suspension has only short-range order arising from weakly flocculated groups of particles held together by the shallow secondary minimum of the interparticle potential. With increasing shear rate, the number of nearest neighbors decreases steadily, indicating breakup of the flocs. Despite this, the stress increases with shear rate (Figure 3) because particles are driven slightly closer together by the irrotational component of the shear flow, and they interact more strongly. Above +c, the behavior of the two systems differs. In the larger system, in spite of the discontinuous increase in the number of nearest neighbors, the excess viscosity and stress both drop abruptly at qC. This happens because most of the nearest neighbors are now nearly coplanar, and coplanar interactions do not contribute significantly to the stress in the suspension. In the smaller system, the number of coplanar nearest neighbors rises gradually, not abruptly, for > +c. This behavior correlates with the gradual decline in the stress seen earlier in Figure 3. The stress does not reach its minimum until the degree of string ordering is very high, and in the smaller system this requires higher shear rates. This tendency is further exaggerated with decreasing system size. In a series of runs with 32 particles, we found the string phase to be poorly defined even at +T = 500, the highest shear rate used.32 Thus, in small systems the transition is diffuse, and behavior seen in the larger system at lower shear rates is postponed to higher shear rates. Now we return to consider the basis for the thixotropic behavior shown in Figure 4. This figure contains plots of subaverages of uxyversus time that are valuable for assessing when a true steady state has been achieved.33 As noted earlier, the run at ?T = 100, which was started with particle coordinates from the end of a run at +r = 50, shows stationary behavior since no structural reorganization of the suspension occurs in passing between these shear rates. However, the run at +T = 150 required a long relaxation period before a steady state was achieved. The explanation for this behavior is that the initial disordered state (the starting coordinates were from the end of the +r = 100 run) persists for roughly 0.421. During this induction period, the viscosity and stress remain essentially constant while the system searches through configuration space for a pathway to the ordered string phase. The reorganization of the particles then proceeds for another time interval of 0.47 until the stable ordered phase is
2.5
2.9
3.3
Figure 8. Spherically averaged pair correlation functions for N = 256 at different shear rates: (a) 100 and (b) 750.
Above the critical shear rate, the layering in the flow direction ( x ) contrasts sharply with the random particle arrangement below the transition. The particles in these layers, defined by a narrow spread of y coordinates, are found in an imperfect hexagonal packing, as seen in Figure 7d, which is the xz projection of such a layer for a shear rate of 750. The particles in one of these ‘xz planes” all move at nearly the same x velocity, the local convective velocity, so we can imagine these layers of particles moving relatively rigidly past each other. Due to the short range of the interparticle force, all of the stress in the ordered suspensions arises from the interaction of particles in adjacent layers. The order in the suspension is also exhibited by the spherically averaged pair correlation function g(r) in Figure 8 for N = 256 at two different shear rates: +T = 100 and 750. For low shear rates (Figure 8a), the randomness of the Brownian particles results 1 and a small nextin a single nearest-neighbor peak at r l d nearest-neighbor ‘bump” at r l d 2. Above a shear rate of 125 (Figure 8b), the suspension has become ordered, and the character of g(r) changes. New “bumps” appear at several r / d values and grow more defined at high shear rates. The three peaks at r l d
- -
(32) Rigos, A. A,; Wilemski, G. Unpublished results. (33) The discrepancies between results presented here and the previous results” arise from the fact that, in the earlier work, runs of 40 OOO time steps were deemed sufficiently long without evaluating the stability of the final state.
Rigos and Wilemski
3986 The Journal of Physical Chemistry, Vol. 96, No. IO, 1992
realized. For the last time period of 0.257, the suspension is ordered, and the amplitude of the fluctuations is greatly reduced. The opposite behavior is demonstrated by the run at a shear rate of 125, begun with particle coordinates from the end of the run with j 7 = 150. The suspension is initially ordered, and the stress is relatively low. However, the shear rate of 125 is not large enough to maintain the ordered state. The Brownian motion of the particles results in the progressive disordering of the suspension, and the stress reaches a steady value after a time of about 0.457. From these two cases, we see that when the shear rate is abruptly changed from a value above j c to a value below, or vice versa, the stress and excess viscosity change slowly, on a macroscopic time scale as the structural transition proceeds. This thixotropic behavior is due here to the simple physical mechanism of particle diffusion. Consider the disordering transition which occurs via the yz diffusion of particles into the 'free volume" between strings. With the particle strings initially separated by about 1.03d, each particle must diffuse a distance p of about 1.03d/(3)Il2 for full disordering to occur. For two-dimensional diffusion, the time tD required for a particle with diffusion coefficient D to incur a mean-squared displacement of p2 is given by the familiar relation pz = 4DtD. Evaluating D with the Stokes-Einstein expression, we find that t D 0.35s. This diffusion time compares favorably with the times found in the simulation and indicates that many-body interactions do not play a decisive role in the disordering dynamics, most likely because of the short diffusion distances and short-range forces involved. Finally, we note a curious quirk in the behavior of the 108particle system: the appearance of two distinct string packing patterns. One type of packing can be seen in Figure 5c,d and Figure 7 of ref 17. The other type resembles that found for the 256-particle system as seen in Figure 6c,d. The patterns differ by a 30' rotation as seen by comparing Figure 5d with Figure 6d. Once formed, each packing pattern was indefinitely stable within the high shear rate region; we presume this was due to the role of the periodic boundary conditions in preserving the pattern that formed first. Comparison runs made at shear rates of 400 and 500 showed that the two packings consistently gave slightly different average stresses and, therefore, excess viscosities. This behavior is clearly an artifact since the steady-state viscosity or suspension structure should not be multiply valued. For the 256-particle system, only one type of string packing was observed.
-
Conclusions We have explored the system size and time dependence of shear thinning in nonequilibrium Brownian dynamics simulations of a thixotropic, sterically stabilized colloidal suspension. The order-disarder shear thinning transition, identified previ~usly,'~ was found to be sharp with 256 particles and diffuse with 108. We believe that Figures 2 and 3 contain the first discontinuous shear thinning curves generated with NEBD for steady shear flow. Thixotropy (time-dependent viscosity) was associated with the time needed to organize or break down the ordered string structure depending on the direction of the transition. This time dependence correlated well with the time scale for free particle diffusion. The 108-particle system had to be more strongly driven convectively before achieving a highly ordered string packing. This is re-
sponsible for the diffuseness of the transition in the small system. We also found two different ordered states with slightly different stresses and excess viscosities for the 108-particle system. This artifact may be a consequence of small system size and/or the simulation algorithm. One obvious implication of our results is that the small system sizes (32 and 108 particles) used previously14J5J7must be avoided in NEBD simulations of shear ordering in colloidal suspensions. Hanley et al.23have noted previously that nonunique, ordered steady states occurred in NEMD simulations when the systems were improperly thermostated. The artifacts we observe cannot be due to thermostating, which is not needed in Brownian simulations, but they may nonetheless be a consequence of the simulation algorithm which is strictly valid only for dilute systems. The simulation algorithm may also be responsible for the appearance of the string phase. As seen in eq 2, with the neglect of hydrodynamic interactions, the convective contribution to the particle displacement varies linearly with they coordinate of the particle, and fluctuations in the local fluid velocity field are neglected. The particles are thus driven to assume a linear velocity profile as if they were hydrodynamically independent, mimicking somewhat the effect of biased themostating found in the NEMD simulation^.^^*^^ In concentrated systems, the particles will hydrodynamically interact strongly at close range, and the resulting velocity fluctuations may be large enough to prevent string formation even though they contribute only modestly to the suspension viscosity. Tentative, indirect evidence for this proposition may be found in the monolayer simulations of Bossis and Brad ~ which ? include ~ hydrodynamic ~ ~ interactions. They find no evidence of string phases over a very wide range of shear rates unless the monolayer is highly concentrated or very strongly interacting.24 However, their systems are small, and we have found string phases to be poorly defined in three-dimensional simulations of 32-particle systems.32 In large two-dimensional NEMD simulations, string phases are found when a linear velocity profile is imposed via thermostating, and they are absent when unbiased thermostating is used.22923It would indeed be disappointing if the interesting ordered string phase found in NEBD simulations were an artifact due to neglect of hydrodynamic interactions, but a decisive answer to this question may be difficult to obtain for large three-dimensional systems because of the formidable computational effort needed to include hydrodynamic interaction^^^-^^,^^ in the simulations. Acknowledgment. We dedicate this paper to Marshall Fixman whose many contributions to science constantly inspire us. Computer time was provided by the Memmack College Computer Facility (MCCF) and the Supercomputer Facility at the Institute for Molecular Science (IMS), Okazaki, Japan. A.A.R. thanks Rand Hall from MCCF, Dr. Umpei Nagashima and Kazuhiko Honda from IMS for their computational assistance, and Mark Gelinas for his assistance in the plotting of the data. (34) Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 1352. Fixman, M. J. Chem. Phys. 1978,69, 1527, 1538; Macromolecules 1981,14, 1710. Ansell, G. C.; Dickinson, E.; Ludvigsen, M. J. Chem. SOC.,Faraduy Trans. 2 1985,81,1269. van Megen, W.; Snook, I. J . Chem. Phys. 1988,88, 1185.
