Spatial-Temporal Correlations at the Onset of Flow in 

Concentrated Suspensions 

Nicos S. Martysa, Didier Lootensb, William L. Georgec, Steven G. Satterfieldc, Pascal Hébraudd

a

National Institute of Standards and Technology,100 Bureau Drive, Stop 8615, Gaithersburg, MD 20899-8615 b

SIKA TechnologyA.G., Tuffenwies 16, CH-8048 Zurich, Switzerland c

National Institute of Standards and Technology,100 Bureau Drive, Stop 8615, Gaithersburg, MD 20899-8615 

 dIPCMS,UMR 7504,23 rue du Loess,67034 Strasbourg Cedex 2, France 

Abstract. Spatial-temporal correlations in the startup-up flow of colloidal suspensions composed of attractive spherical 

particles under static and shear conditions are studied by computer simulation. The evolution of particle organization with 

time is followed as equilibrium is approached. The time dependence of the radial distribution and organization of nearest 

neighbors is tracked as equilibrium is approached and then as a constant shear rate is applied. 

Keywords: rheology, colloidal suspensions, Van der Waals forces, coordination number, structure, thixotropy. 

PACS: 82.70.Dd, 8360.Rs, 83.80.Hj 

Colloidal suspensions [1] are of fundamental importance for many applications, such as cementitous materials. 

They exhibit a rich variety of behavior such as shear thinning, viscoelasticity, shear thickening, and jamming. Here 

we study how a colloidal system, starting from rest, responds to an applied strain. The initial organization of the 

microstructure and its subsequent change under applied stress, play an important role in determining the rheological 

properties of the suspension. Experimental study of the microstructure during the flow can be made with neutron 

and light scattering techniques, but these techniques require long time averages making it difficult to study the early 

stages of flow. In this paper, we present results of a computational study of concentrated colloidal systems under 

flow making the link between microstrucure and rheological properties. 

The simulation approach used in this study is based on a modified dissipative particle dynamics method [2]. This 

model accounts for lubrication forces when spherical aggregates are very close and includes a Van der Waals 

attractive force with an effective hard sphere repulsive force. The interaction potential chosen for this study has a 

depth of 25KBT with spherical particles of size of 100 nm. Three volume fractions of spheres, 20 %, 40 % and 50 %, 

are studied using mono-size spheres numbering 3760 to 9616. The initial microstructure is established without 

interaction: particles are randomly packed and then allowed to equilibrate for a short period. During this 

equilibration stage, the Van der Waals force is gradually turned on in order to avoid the introduction of a large 

kinetic energy as the spheres could be very close initially. Once the forces are fully “turned on” we define time zero. 

The suspension is then allowed to undergo a secondary equilibrium phase that can be related to a thixotropy. We 

then study how the suspension responds to the introduction of a constant applied strain rate at different stages of 

secondary equilibration. 

 The temporal evolution of the coordination number, Z, is first considered. Z is obtained from the radial 

distribution function, which is integrated from zero out to a small distance associated with the first minimum of the 

radial distribution function. Figure 1 shows the growth of Z as a function of time. For all volume fractions, Z tends 

to a plateau value of approximately 6.6. A constant strain rate is then applied at different equilibrium time. We find 

that under lower shear rates Z increases slowly as if the shearing is enhancing the aging process. On the other hand, 

at higher shear rates, the coordination number decreases approximately following an exponential decay. 
10-4 10-3 10-2 10-1 100

4.0

4.2

4.4

4.6

4.8

5.0

5.2

5.4

5.6

5.8

6.0

6.2

6.4

6.6

Z

TIME

2.2

10-2 10-1 100

4.0

4.2

4.4

4.6

4.8

5.0

5.2

5.4

5.6

5.8

6.0

6.2

6.4

6.6

 Z

DEFORMATION

10-4 10-3 10-2 10-1 100

4.0

4.2

4.4

4.6

4.8

5.0

5.2

5.4

5.6

5.8

6.0

6.2

6.4

6.6

Z

TIME

22

10-3 10-2 10-1 100

4.0

4.2

4.4

4.6

4.8

5.0

5.2

5.4

5.6

5.8

6.0

6.2

6.4

6.6

 Z

DEFORMATION

10-4 10-3 10-2 10-1 100

4.0

4.2

4.4

4.6

4.8

5.0

5.2

5.4

5.6

5.8

6.0

6.2

6.4

6.6

Z

TIME

220

10-3 10-2 10-1 100

4.0

4.2

4.4

4.6

4.8

5.0

5.2

5.4

5.6

5.8

6.0

6.2

6.4

6.6

 Z

DEFORMATION

FIGURE 1. Coordination number as a function of the time (logarithmic scale) and deformation (strain). Time is in units

corresponding to the sphere diffusing a distance equal to its radius. Shear rate is equal 2.2, 22.0 and 220.0 inverse time units 

(from the left to the right).

Stress growth curves are calculated as a function of the equilibration times and are represented in Figure 2 for a 20

% volume fraction suspension. For lower equilibration times, the stress grows slowly with strain until a critical

strain is reached where a more rapid increase ensues. This is indicative of the fact that the suspension has not 

equilibrated long enough to form a spanning cluster and that the small strain facilitates the formation of clusters. 

However, for longer equilibration times the stress rises more rapidly with strain. Even though the slope of the stress 

growth curve increases, the maximum reached at the equilibrium does not change significantly indicating that the

cluster formation is still tenuous. It is possible to collapse all the curves onto a universal master curve by a linearly 

translating the stress growth curve along the strain axis. Hence the shearing produces an “aging” effect for the 

suspension as the shearing facilitates the formation of larger clusters. 

0.0 0.2 0.4 0.6 0.8 1.0

5

10

15

20

25

30

 0.015

 0.15

 0.61

 0.91

 1.83

σ/μγ

STRAIN

equilibrium time

0.0 0.2 0.4 0.6 0.8 1.0 1.2

5

10

15

20

25

30

 0.015

 0.15

 0.61

 0.91

 1.83

σ/μγ

STRAIN (SHIFT)

equilibrium time

FIGURE 2. Stress growth curves for a 20 % volume fraction suspension where the strain is applied at different points of

equilibration. By shifting the strain, the stress growth curves can be made or overlap, forming a universal stress growth curve. 

On the other hand, at higher volume fractions the stress growth curves take a different form as the system ages 

(Figure 3). At early times, the stress growth curve increases monotonically but at later times it is found, initially, to 

grow much more rapidly. This rapid increase of stress with strain indicates that the suspension has to overcome a 

barrier associated with the rearrangement of its microstructure. Indeed, at higher volume fractions, the stress could

increase rapidly and then drop as quickly. The maximum of the stress can be interpreted as a static yield stress. 

Let us now consider the angular distribution of radial contacts between spheres (Figure 4). In the low volume

fraction range, we find that the contacts are largely oriented in the extension quadrant. As the spheres are pulled 

apart by the applied strain, they have a tendency to form transient string like objects that break up as it interact with 

nearby clusters. In general the first normal stress is negative for this simulation case. At higher volume fractions and

shear rates, the orientation of the angular distribution of contacts becomes more populated in the compression 

quadrant. Here one can see dramatic fluctuations of stress chains that form as the shear rate is increasing. This 

behavior is akin to giant stress fluctuations seen in suspension at the onset of a jamming transition. In this regime,

the first normal stress was found to be positive. 
0.0 0.2 0.4 0.6 0.8 1.0

0

20

40

60

80

100

 0.015

 0.15

 0.61

 0.91

 1.83

σ/μγ

 STRAIN

EQUILIBRIUM TIME

0.0 0.2 0.4 0.6 0.8 1.0 1.2

0

20

40

60

80

100

 0.015

 0.15

 0.61

 0.91

 1.83

σ/μγ

 STAIN (SHIFT)

EQUILIBRIUM TIME

FIGURE 3. Stress growth curves for a 40 % volume fraction suspension where the strain is applied at different points of

equilibration. The stress growth curves could not be made to overlap by shifting the strain. 

FIGURE 4. Angular distribution of contacts for volume fraction 20 % (left) and 50 % (right). The shear flow is in the x

direction. The population of angular contacts are projected in the vorticity, or xy plane (blue), yz plane (green) and zx plane

(violet). For the 20 % volume fraction simulation, the contact populations projected in the vorticity plane were largely found in 

the extension quadrant (between 0o

 and 90o). In the 50 % volume fraction, the contacts largely populate the compression

quadrant(between -90o

 and 0o). 

In summary, at low volume fractions, a slowly applied strain had a tendency to facilitate aging of the suspension. 

Stress growth curves could collapse onto a universal curve. The maximum stress was not dependent on the aging

and could be reached sooner with the applied strain. At volume fractions close to 50 %, the stress growth curves did 

not exhibit a universal behavior. As the suspension aged, a barrier had to be overcome to begin flow. This barrier

can be thought of as a yield stress. For this case, contacts largely oriented in the compression quadrant and 

fluctuating stress chains were observed as the shear rate was increased. 

ACKNOWLEDGMENTS 

We would like to gratefully acknowledge support from the Virtual Cement and Concrete Testing Laboratory

consortium (VCCTL). The flow simulations were performed under award SMD-05-A-0129, “Modeling the 

Rheological Properties of Suspensions: Application to Cement Based Materials,” for NASA’s National Leadership

Computing System Initiative on the “Columbia” supercomputer at the NASA Ames Research Center. 

REFERENCES 

1. R. G. Larson, The Structure and Rheology of Complex Fluids, New York: Oxford University Press. 

2. N. S. Martys, J. Rheol., 49(2), 401-424, (2005).