Contact and stress anisotropies in start-up flow of colloidal suspensions

Nicos S. Martys1, Didier Lootens2, William George1, and Pascal H ́ebraud3

1 National Institute of Standards and Technology,

100 Bureau Drive, Stop 8615, Gaithersburg, MD 20899-8615, USA

2 SIKA Technology A.G., Tuffenwies 16, CH-8048 Z ̈urich, Switzerland

3

IPCMS UMR 7504 23 rue du Loess 67034 Strasbourg Cedex 02, France

Spatio-temporal correlations in start-up flows of attractive colloids are explored by numerical

simulations as a function of their volume fraction and shear rate. The suspension is first allowed

to flocculate during a time tw, then the stress necessary to induce its flow is computed. We find

that, at low volume fractions, the stress is a universal function of the strain. On the contrary,

at high volume fractions, this scaling behavior is no longer observed and a supplementary stress

becomes necessary to induce flow. To better understand the physical origin of the supplementary

stress, we examine the creation, disruption and orientation of contacts between the particles, and

the corresponding contribution to stress as a function of strain. Our simulations show that the onset

of flow is dominated by the creation of contacts between the particles at low shear rates, and by

their disruption at high shear rates. However, neither the evolution of the number of contacts with

strain, nor their orientation can fully account for the non-scaling behavior of the stress at high

volume fractions. At small strains, the relative importance of forcing in the compression quadrant

increases with volume fraction, and with flocculation time. This mechanism of stress transmission

through the compression quadrant is not accounted for in the usual description of yield stress, which

considers the breaking of bonds oriented in the extension quadrant.

PACS numbers: 83.80.Hj,87.10.Tf,83.60.Df

1. INTRODUCTION

Many complex fluids such as dense suspensions, emulsions or colloidal glasses are composed of interacting particles

that self-organize into an isotropic microstructure under quiescent conditions. The application of a rate of deformation

disrupts this isotropy such that a collective behavior emerges. As a consequence, a rich variety of structures emerge [1–

3], that determine mechanical macroscopic flow properties characterized by non-linear relationships between the stress

and the shear rate : among them, plastic behavior [4] (a power-law relationship between the stress and the shear rate),

jamming [5–7] (flow stops under the application of a high enough stress) and yielding [8] (flow does not start although

a non-zero small stress is applied). It is of crucial importance to understand how these materials start to flow and

to tune their physicochemical properties in order to reduce the amount of energy necessary to make them flow. An

outstanding question about such systems thus emerges : what physical mechanisms are responsible for the observed

macroscopic behavior as the initial isotropy is broken ? For example, mechanical measurements reveal a very nonlinear relationship between the applied stress and the induced shear strain or shear rate : it is necessary to apply a

non-zero value of the stress, the yield stress, in order to induce the flow. The microscopic origin of the yield stress

is not generally understood, and its value, its measurement, and even its existence are not well established [9, 10].
2

The yield stress is generally thought as a manifestation of the maximum force needed to break contacts between

particles [4]. But, when a suspension flows, particles are not only drawn apart, but also pushed together. That is,

under simple shear, two colloids whose center-to-center segment is oriented in the extension quadrant are drawn apart,

whereas when they approach each other due to the flow, their center-to-center segment is oriented in the compression

quadrant (Fig. 1). In this paper, we show that this phenomenon leads to a new physical mechanism that produces an

excess of stress at small strains. Indeed, as flow develops, aggregates of particles form [11, 12]. These aggregates are

not isotropic, and their orientation relative to the flow direction depends on the shear rate and the volume fraction

of the suspension, as may be seen from direct observation of the particles inside the shear cell (fig. 3). We show that,

at high volume fraction, a supplementary stress develops at the onset of flow. This supplementary stress is due to the

organization of contacts between particles in the compression quadrant of the shear flow, and it is not fully accounted

for by the simple breaking of colloidal contacts.

The paper is organized as follows. The computational model is briefly introduced and the preparation of the suspension,

before the application of shear is described. A constant shear is then applied and it is shown that, at low concentrations,

the growth of stress with strain exhibits a universal behavior. This universal behaviour is not observed at small

strains and the higher concentration studied. Indeed, an additional stress, supplementary to that indicated by simple

scaling arguments, is needed to make the suspension flow. The microstructure of the suspension responsible of this

supplementary stress is then investigated in the last part of the paper, where the orientations of pairs of particles in

contact, and the stresses between them are studied.

2. THE DPD MODEL

We use a dissipative particle dynamics (DPD) based model [16] to study the onset of flow of a colloidal suspension

of attractive particles. This approach can be thought of as a Lagrangian formulation of Navier Stokes with thermal

fluctuations [24]. As this computational model is largely described in a previous work [16], we will only highlight its

main features and focus, more so, on the relevant aspects of the model to this work. The DPD simulation is similar in

structure to molecular dynamics simulation, but, instead of modeling all the molecular properties of the system, the

motion of mesoscopic DPD particles that represent a coarse grained fluid, are considered [15, 16]. The DPD particles

are subjected to dissipative, random and conservative forces. The dissipative and random forces control the viscosity

of the fluid and, to maintain a well defined temperature, are related by fluctuation-dissipation theorem [23]. The

conservative force is utilized to control the compressibility of the system, so that it matches that of water [14]. In

contrast to the individual DPD particles, a colloid is then defined as an assembly of constrained DPD particles so

that they form a rigid body (Fig. 2). As described previously [16, 21], when modeling a dense suspension of hard

spheres (or colloids in our case), the DPD particle interactions are not sufficiently strong to prevent overlaps of the

spheres so this model includes lubrication forces [16, 17]. The computational model, employed in this study, has been

shown to recover the Einstein intrinsic viscosity for a dilute suspension of spheres [15, 16], reasonable estimates of the

Huggins coefficient for moderately dense suspension [16, 21] and good agreement with Stokesian dynamics predictions

of relative viscosity as a function of shear rate in dense suspensions [16, 20].

For this study, an additional interaction potential is incorporated into the model. Here, we use a Derjaguin type

approximation of the Lennard-Jones interactions [18] between colloid spheres A and B at short distances. This
3

potential scales as the sum of a hard sphere repulsion term, AHS/s7

AB and a Van der Waals potential, H/sAB,

where AHS and H are constants related to the hard-sphere potential and to the Hamaker constant, and sAB is

the distance between the two particles surfaces. The colloid’s radius is the only length scale explicitly included in

the description of the suspension. The parameters chosen for the simulation could correspond to spherical alumina

particles of approximately 100 nm diameter. The interaction potential between two particles exhibits a minimum of

25 kBT at a distance of 6 nm between the surfaces of the particles. The time integration scheme for the DPD particles

and the colloids is based on a velocity-Verlet algorithm described in [16, 22].

The stress tensor for this model suspension has several contributions. There are contributions from the propagation

of momentum and interparticle forces of the DPD particles and is given by :

σαβ =

1

V m

X

i

p ̃iαp ̃iβ +

1

2V

X

i,j

f

α

ij (ri − rj )β (1)

where i and j refer to two different DPD particles, α and β to cartesian coordinate axes, fij is the force between

particles i and j, and  ̃piis the momentum of particle i relative to the macroscopic velocity field. V is the total volume

of the system and m the mass of the DPD particles. In addition there are contributions from the colloidal interactions.

The stress contribution from the Van der Waals forces is given by :

1

2V

X

A,B

F

α

AB(rA − rB)β (2)

where A and B refer to two different colloids and FAB the Van der Waals forces between colloids A and B described

above. Stress contributions from the lubrication forces are also included. They imply several modes. The most important one is the squeeze mode, which accounts for forces that develop as two spheres directly approach each other.

This force is proportional to the velocity difference between the spheres and is inversely proportional to the nearest

surface to surface distance. For the case of monosize spheres the lubrication force

Flub =

3

2

πηa2(vA − vB)/sAB (3)

where η is the continuous fluid viscosity, a is the sphere radius, vA is the velocity of colloid A and sAB is the nearest

surface to surface distance between spheres A and B. Additional contributions to the squeeze mode as well are other

modes all scale logarithmically and are largely dominated by the squeeze mode [17]. In our computations, all the

modes of the lubrication force are taking into account up to the first order, including terms that scale as 1/sAB,

ln sAB and sAB ln sAB.

There are additional corrections to the stress tensor due to the constraint forces on the DPD particles that make up

the colloid. The constraint forces are determined by accounting for the rigid body motion in the individual particles

displacements and change in velocity at each time step of the velocity Verlet algorithm. Details of this approach are

described in [22].

We study three volume concentrations, φ = 20 %, 40 % and 50 % where the number of colloidal particles ranges

from 3760 to 9616. As we are interested in the onset of flow we study flow in a simple shear geometry with shear

strain γ smaller than 1. In each experiment, we control the shear rate  ̇γ, and measure the macroscopic stress σ..
4

3. PREPARATION OF THE SAMPLE

The colloidal particles are attractive, and in the absence of flow, the suspension evolves to form a colloidal gel. We

follow a well-defined preparation protocol of the suspension. It is prepared by randomly placing the spheres, in a cubic

cell, with small overlaps being possible. Then, a repulsive force is applied between spheres if they overlap. Once all the

spheres no longer overlap, the system is allowed to equilibrate under the action of Brownian and hydrodynamic forces.

As the system equilibrates, the pair distribution function is evaluated and found to become reasonable consistent

with the Percus Yevick [19] approximation. When a stable radial distribution function is obtained, the Van der

Waals and lubrication forces are introduced by turning them on slowly enough so that the additional kinetic energy

produced in the system, because two spheres are too close each other, is allowed to dissipate. We choose the end of the

introduction of the Van der Waals and the lubrication forces as the initial time of the system. Then, for a given time

tw (tw = 9.3, 93, or 1120), the suspension is allowed to evolve towards a flocculated colloidal gel. As the suspension

ages during a time tw, it is still out of equilibrium, and during this period, the pair distribution function evolves

(fig. 4(a)) : it becomes increasingly peaked at contact and develops secondary peaks. The increase of g(r) at contact

with time (fig. 4(a) inset), is indicative of the increase in number of colloids in near contact [25]. For φ = 20 %, a

marked peak at a distance between colloids’center equal to √

3r is obtained, showing that a local hexagonal structure

develops (Fig. 4(b)).

4. RESCALING OF THE FLOW CURVES

At several tw values, we then begin to make the suspension flow. A shear deformation is imposed by translating the

upper and lower boundaries of the cube according to the Lee Edwards boundary condition [27]. The direction of the

flow is denoted ux and its gradient uy (Fig. 1). In order to compare our results with measured stress and shear rate

values, we need to identify the characteristic stress  ̃σ and time (or shear rate, γ ̃ ̇) units of the flow. As the interparticle

interaction is strong, the relevant energy scale in our system is the interaction potential between two particles, rather

than thermal energy, so that in the studied flow regime the interaction forces between two particles compete with

hydrodynamic forces. Given that Fmax is the maximal interaction force between two colloids, a their radius, η the

viscosity of the continuous medium, the natural stress, shear rate (and thus also time) units of our flow are given by :

σ ̃ = φ

2 Fmax

a

2

(4)

γ ̃ ̇ =

Fmax

6πηa2

(5)

where  ̃σ is the product of Fmax with the number of particle bonds that cross the unit area, φ

2/a2

[26], γ ̃ ̇ is the shear

rate, that equilibrates the intercolloidal stress, Fmax/a2, with the hydrodynamic stress, 6πηγ ̇ . Our simulation results

may thus be compared to real systems in which the energy of interactions between colloids is larger than thermal

energy. Considering the 100 nm alumina particles, the maximal interaction force, given by the maximum of the slope

of the interaction potential is of the order of 4.10−13 N, so that the stress scale  ̃σ is 40 Pa. The corresponding shear

rate γ ̃ ̇ is 210 s−1.

The response of the suspension to the application of shear is given in Fig. 5. At extremely low deformations, we

do observe an elastic regime, in which the stress is proportional to strain, and goes to zero at null shear values.
5

Nevertheless, thermal activation plays a role in this elastic regime, as the apparent elastic modulus measured

decreases when the shear rate decreases. For this paper, we are not interested in this elastic regime, but will describe

in detail the onset of flow regime. The value of the stress at the onset of flow depends on the equilibrium time : the

longer tw, the higher the initial stress value. Then, the stress progressively increases up to a maximum value that

does not depend on the equilibration time.

In addition, the stress growth curves of low concentration suspensions exhibit a remarkable feature : if the origin of

strain is chosen to take a non-zero value, γ0 , then all the stress, as a function of strain, superimpose (fig. 5, right

column). γ0 is an increasing function of the waiting time tw. As a consequence, the application of flow plays the same

qualitative role as the aging of the suspension. Such self-similarity is characteristic of the aging dynamics of out of

equilibrium systems, and has been observed in the case of concentrated soft colloidal particles [28].

In contrast to the low concentration suspensions, for the higher concentrations studied and when the waiting time

is long enough, the evolution of the stress as a function of time exhibits a different shape. In this case, the stress

growth curves may no longer be superimposed by shifting the strain. At long waiting times, the value of the stress,

(ignoring the elastic regime), at null strain now tends towards a non-zero finite value. The stress, needed to make

the suspension flow, exceeds that which would be necessary according to the aging behavior. This departure from the

scaling behavior occurs once the flow of the suspension has been started, but before the establishment of a steady flow

regime at larger shear values. The excess of stress, relative to the scaling behavior, needed to make the suspension

flow thus differs from the idealized concept of yield stress, which is a non-zero value of the stress required to induce

the flow, as  ̇γ −→ 0, and we call this excess stress a ”supplementary stress”. The higher the suspension concentration

and the higher the shear rate, the greater the supplementary stress needed. As yield stress is generally associated with

the breakage of bonds between attractive particles [4], we look for the microscopic origin of the supplementary stress

that we measure macroscopically, and consider :

– the evolution of the number of contacts between particles when shear is applied,

– the orientation of contacts under shear,

– the angular distribution of interparticle stress between particles in contact.

5. MICROSTRUCTURAL ORIGIN OF THE SUPPLEMENTARY STRESS

5.1. Evolution of the number of contacts

The microstructure of the suspension may be described with the help of the pair distribution function, g(r). We

have observed that it exhibits a first peak at contact, near r = 2a, and a series of 2 correlation peaks. As the

stress is localized in the thinnest gaps between particles [29], we focus our analysis on the particles in near contact

and consider the pair distribution function g(r) of the suspensions, at different times as the suspension equilibrates

and as the suspension is later strained. First note, g(r) exhibits a well-defined maximum near r = 2a, where the

particles are at contact. From the value of this maximum, one can compute the average number of contacts Z of the
6

particles [30, 31] :

Z =

N

V

Z 2a+δa

2a

4πr2g(r)dr (6)

a being the radius of the particles, δa the width of the first peak of the pair distribution function, N their number

and V the volume of the shearing cell. In the absence of flow, whatever the volume fraction, the number of contacts

increases with time. It tends towards a plateau value close to 7. The higher the concentration and the longer the

equilibrium waiting time, the faster the plateau value of Z is reached. When we apply a low shear rate deformation,

the number of contacts generally increases with time. On the contrary, at high shear rates, the total number of

contacts decreases as time elapses, before reaching an equilibrium. Whatever the applied shear rate value, we found

that the evolution of the number of contacts as a function of the shear follows a simple exponential behavior (Fig. 6) :

Z(γ) = Z∞ + (Z0 − Z∞)e

−γ/γc

, where Z0 is the number of contacts before the application of shear, Z∞ the number

of contacts once equilibrium is established, and γc is a characteristic strain of contact formation (at low shear rate),

or contact disruption (at high shear rates). At low shear rate, γc is of the order of 100 %, corresponding to the

characteristic shear necessary to induce collision between two colloids. γc thus decreases when the volume fraction

of the suspension increases (Fig. 6d). In contrast, at higher shear rates, the number of contacts decreases, and

disruption of contacts dominates the dynamics. Here, γc no longer depends on the volume fraction of the suspension

(Fig. 6d). Whatever the volume fraction and the shear rate, the evolution of the number of contacts as a function

of shear is always a simple exponential function and we do not observe a significant increase of the contact number that could be responsible for the non-scaling behavior of the flow curves when the suspension has aged long enough.

5.2. Orientation of contacts

The evolution of the number of contacts alone cannot explain the appearance of a supplementary stress. We

next consider their orientation relative to the shear direction. We do not observe any orientation of the direction of

contacts either in the plane normal to the shear direction, or in the plane normal to the shear. In the shear/shear

gradient plane, the contacts become anisotropic under flow (fig. 7). A very different behavior is observed at low

and at high volume fractions. Indeed, for φ = 20 and 40 %, the ratio of the number of contacts in the compression

quadrant over those in the extension quadrant, Ncomp/Next, decreases when the shear increases. On the contrary,

for the highest volume fraction, φ = 50 %, it increases. This observation agrees with experimental observation of the

formation of aggregates in the compression quadrant, close to the direction of the flow [2].

More remarkably, the contacts anisotropy curves for tw = 9.3 and tw = 1120 superimpose (fig. 7), whereas the

measured macroscopic stresses exhibit very different behavior, the longest aged suspension no longer obeying to the

scaling behavior. Once again, this observation rules out the hypothesis that the contact anisotropy is responsible for

the excess of stress.
7

5.3. Angular distribution of interparticle stress between particles in contact

We thus consider the anisotropy of the stress stored between contacts and define σ

comp

12 /σext

12 =

R 0

−π/2

σ12(θ)dθ/ R π/2

0

σ12(θ)dθ, the ratio of stress between contacts in the compression over the extension quadrant.

The stress between contacts for the highest concentrated and longest aged suspension (φ = 50 %, tw = 1120), accumulates in the compression quadrant of the flow, whereas stress is only localized in contacts oriented in the extension

quadrant for the less aged suspension (fig. 8a). More generally, for φ = 40 % and φ = 50 %, the stress along contacts in

the compression quadrant no longer tend towards 0 when the strain γ goes to 0 (fig. 8b). Thus, the macroscopic stress

necessary to induce the flow is stored along contacts oriented in the compression quadrant. This result is important

as it complements the general observation that the yield stress is due to the breaking of contacts in the extension

direction [4]. As shown in fig. 8a, for concentrated flocculated suspensions, the creation of contacts and their setting

under tension in the compression quadrant of the flow is quantitatively as important as their breaking in the extension

direction.

A macroscopic consequence of the the creation and disruption of contacts between colloids may also be found in the

evolution of the stress with strain of the suspension when the depth of the Van der Waals potential minimum is

different. We have performed simulations for φ = 20 % and φ = 50 %, for two different interaction potential depths,

25 kBT and 125 kBT. Although one may expect that the maximum stress would increase by the same factor of the

depth of the interaction potential between colloids, we observe that while it increases by a factor of 4.6 for φ = 20 %,

there is only a factor of 3.7 increase for φ = 50 % (fig. 9). This indicates that some stress is not only due to the Van

der Waals interaction between the colloids, but also due to the creation and disruption of contacts as the colloids

rearrange themselves under shear. Such forces play a larger role when the suspension is more concentrated.

6. CONCLUSION

In summary, the beginning of the application of shear strain to a colloidal suspension induces an anisotropy of the

orientation of contacts between particles. This anisotropy is qualitatively different at low and high volume fractions. In

the first case, an excess of contacts develop, whose relative orientation lies in the extension quadrant of the shear flow,

in the latter, there is an increase of particles whose relative orientation is in the compression quadrant. At low volume

fractions, the stress necessary to induce flow is the stress needed to break bonds, and microscopic stress is stored in

the extension quadrant, which corresponds to the standard explanation of the origin of yield stress. Nevertheless, at

high volume fraction, a new physical mechanism plays a role in the appearance of yield stress : the accumulation of

stress in the compression quadrant of the flow.

Acknowledgements 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 :

Applications to Cement Based Materials” for NASA’s National Leadership Computing System Initiative on the ”Columbia” supercomputer

at the NASA Ames Research Center. This research used resources of the Argonne Leadership Computing Facility at Argonne National
8

Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract DE-AC02-06CH11357.

[1] P. Varadan and M.J. Solomon Direct visualization of flow-induced microstructure in dense colloidal gels by confocal laser

scanning microscopy. Journal of Rheology, 47, 4, 943-968, 2003.

[2] B.J. Maranzano and N.J. Wagner Flow-small angle neutron scattering measurements of colloidal dispersion microstructure

evolution through the shear thickening Journal of Chemical Physics, 117, 22, 10291-10302, 2002.

[3] P. Schall, D.A. Weitz, F. Spaepen Structural rearrangements that govern flow in colloidal glasses. Science, 318, 1895-1900,

2007.

[4] R.G. Larson The structure and rheology of complex fluids. Oxford University Press, New York, 1999.

[5] H.A. Barnes Shear-thickening (”Dilatancy”) in suspensions of non-aggregating solid particles dispersed in newtonian

liquids. Journal of Rheology, 33, 2, 329-366, 1989.

[6] J. Bender and N.J. Wagner Reversible shear thickening in monodisperse and bidisperse colloidal dispersions. Journal of

Rheology, 40, 5, 899-916, 1996.

[7] D. Lootens, H.Van Damme, and P. H ́ebraud. Giant stress fluctuations at the jamming transition. Physical Review Letters,

90, 178301, 2003.

[8] H.A. Barnes, K. Walters The Yield Stress Myth. Rheologica Acta, 24, 323-326, 1985.

[9] P.C.F. Møller, J. Mewis and D. Bonn Yield stress and thixotropy : on the difficulty of measuring yield stresses in practice.

Soft Matter, 2, 274-283, 2006.

[10] J.R. Stokes and J.H. Telford Journal of Non-Newtonian Fluid Mechanics, 124, 137-146, 2004.

[11] W. Wolthers, D. Van den Ende, V. Breedveld, M.H.G. Duits, A.A. Potanin, R.H.W. Wientjes and J. Mellema Linear

viscoelastic behavior of aggregated colloidal dispersions. Physical Review E, 56, 5, 5726-2733, 1997.

[12] T. Serra, J. Colomer and X. Casamitjana Aggregation and breakup of particles in a shear flow. Journal of Colloids and

Interface Science, 187, 466-473, 1997.

[13] R.J. Flatt, N. Martys and L. Bergstrom The rheology of cementitious materials. MRS Bulletin, 29, 5, 314-318, 2004.

[14] R.D. Groot and P.B. Warren Dissipative particle dynamics : Bridging the gap between atomistic and mesoscopic simulation.

Journal of Chemical Physics, 107, 4423-4435, 1997.

[15] P.J. Hoogerbrugge and J.M.V.A. Koelman Simulating microscopic hydrodynamic phenomena with dissipative particle

dynamics. Europhysics Letters, 19, 155160, 1992.

[16] N.S. Martys Study of a dissipative particle dynamics based approach for modeling suspensions. Journal of Rheology, 49,

2, 401-424, 2005.

[17] S. Kim and S.J. Karrila Microhydrodynamics. Butterwort-Heinemann, 1991.

[18] J.N. Israelachvili Intermolecular and Surface Forces Academic Press, 1992.

[19] T.M. Reed and K.E. Gubbins Applied Statistical Mechanics :Thermodynamics and Transport Properties of Fluids

Butterwort-Heinemann, 1991.

[20] D.R. Foss and J.F. Brady ”Structure, diffusion, and rheology of Brownian suspensions by Stokesian dynamics simulations”

J. Fluid Mechanics, 407,167-200, 2000.

[21] E. Boek, P.V. Coveney,H.N. Lekkerkerker, and P.Van der Schoot ”Simulating the rheology of dense colloidal suspensions

using dissipative particle dynamics simulations” Phys Rev E, 55,3124-3133,1997.

[22] N.S. Martys, R.D. Mountain Velocity Verlet algorithm for dissipative-particle-dynamics-based models of suspension Physical Review E, 59, 3733-3736, 1999.

[23] P. Espa ̃nol and P. Warren Thermodynamically admissible form of discrete hydrodynamics. Europhysics Letters, 30,
9

191-196, 1995.

[24] P. Espa ̃nol and M. Revenga Smoothed dissipative particle dynamics Physical Review , 67, 026705, 2003.

[25] R.J.M. d’Arjuzon, W. Frith and J.R. Melrose Brownian dynamics simulations of aging colloidal gels. Physical Review E,

67, 061404, 2003.

[26] W.B. Russel, D.A. Saville and W.R. Schowalter Colloidal Dispersions. Cambridge University Press, 1989.

[27] M.P Allen D.J. Tildesley Computer Simulation of Liquids Claredon, Oxford, 1987.

[28] M. Cloitre, R. Borrega and L. Leibler Rheological aging and rejuvenation in microgel pastes. Physical Review Letters, 85,

22, 4819-4923, 2000.

[29] L.E. Silbert, R.S. Farr, J.R. Melrose and R.C. Ball Stress distributions in flowing aggregated colloidal suspensions. Journal

of Chemical Physics, 111(10), 4780-4789, 1999.

[30] J.-P. Hansen and I.R. MacDonald Theory of simple liquids. Academic Press, 1976.

[31] T. Aste, M. Saadatfar and T.J. Senden Geometrical structure of disordered sphere packings. Physical review E, 71, 061302,

2006.

[32] D. Lootens, H. van Damme, Y. H ́emar and P. H ́ebraud Dilatant flow of concentrated suspensions of rough particles.

Physical Review Letters, 95, 268302, 2005.
10

A B

D C

g

q>0

q<0

y

x

Fig. 1: Strain of a square under a simple shear. Particles A and C, whose center-to-center segment lies in the compression

quadrant (hatching along the first diagonal), are drawn together due to the flow, whereas particles B and D, in the extension

quadrant (hatching along the second diagonal, are pulled apart. θ is the angle of the center-to-center segment and the direction

of the flow.
11

sAB

-20

0

20

U

/

k

B

T

0 10 20 30 40 50 60

sAB(nm)

A

B

i j

Fig. 2: Definition of the DPD particles and of the colloids. The suspension fluid is represented by DPD particles (blue disks).

A and B are two colloids that consist in an assembly of DPD particles, subjected to constraints so that they form a rigid body.

In addition to the the DPD based forces between all (DPD) particles, the colloids A and B are subjected to Van der Waals

repulsive force and lubrication forces. Inset : additional interaction potential between the colloids (Derjaguin approximation of

a Lennard-Jones potential).
12

(a) (b)

Fig. 3: Images from a simulation of a sheared colloidal suspension for volume fraction φ = 0.20 (a) and φ = 0.50 (b) and at

a strain γ of about 0.5. All views are in the shear gradient plane with the top of the suspension moving to the right and the

bottom moving to the left. The figures, offset to the lower right, represent a typical slice through the respective suspension.

Included in the lower images is an indication of stress chains where the largest stresses are found. Note that for the φ = 0.20,

the chains form in the extension quadrant. In constrast, for the φ = 0.50 case, many chains are formed both in extension

and compression. Typically, at the lowest shear rates there are fewer stress chains in compression however there are more by

proportion as the volume fraction and shear rate is increased.
13

4

3

2

1

0

g(r)

1.0 2.0 3.0 4.0

r

4

3

2

1

0

g(r)

1.0 2.0 3.0 4.0

r

10

8

6

4

2

0

1.00 1.10 1.20 1.30

(a) (b)

Fig. 4: (a) Evolution of the pair distribution function as a function of aging time, for φ = 50 %, tw = 0.93 (short dashed line),

tw = 9.3 (long dashed line), tw = 560 (continuous line) and tw = 1120 (thick continuous line). Inset : the pair distribution

functions for the three studied aging times, for particles close from contact. (b)Pair distribution function of the suspensions

after tw = 1120, for φ = 20 % (continuous line), φ = 40 % (short dashed line), φ = 50 % (long dashed line).
14

0.01

2

4

6

0.1

0.01

2 4

0.1

2 4

1 γ-γ0

σsupp

0.01

2

4

6

0.1

0.01

2 4

0.1

2 4

1 γ-γ0

0.01

2

4

6

0.1

0.01

2 4

0.1

2 4

1 γ-γ0

σsupp

0.12

0.08

0.04

0.00

σ

0.0 0.2 0.4 0.6 0.8 1.0

γ

0.12

0.08

0.04

0.00

σ

0.2 0.4 0.6 0.8 1.0

γ

0.12

0.08

0.04

0.00

σ

0.2 0.4 0.6 0.8 1.0

γ

(a)

(c)

(b)

Fig. 5: (a),(b),(c) Dimensionless stress as a function of shear (left column) and time (right column) when a constant shear

rate is applied (  ̇γ = 3.5.10−3). (a) φ = 20 %, (b) φ = 40 %, (c) φ = 50 %. Before the application of shear, the suspension was

allowed to equilibrate for times tw = 9.3 (continuous line), tw = 93 (dashed line), or tw = 1120 (bold continuous line). When

the shear scale is translated by γ0 (right column), all the curves superimpose at low concentration, and exhibit a supplementary

stress σsupp at the two largest concentrations.
15

6.5

6.0

5.5

5.0

4.5

Z

0.1 100

tw

4 6

0.01

2 4 6

0.1

2 4 6

1

γ

6.0

5.5

5.0

4.5

6

0.01

2 4 6

0.1

2 4 6

γ

6.0

5.5

5.0

4.5

Z

100 103

tw

6.4

6.0

5.6

5.2

Z

0.1 100

tw

6

0.01

2 4 6

0.1

2 4 6

γ

(a)

(c)

(b)

0.1

2

4

6

1

2

4

6

10

γ

c

10

2 3 4 5 6 7 8 9

100

φ

(d)

Fig. 6: Evolution of the number of contacts as a function of tw, before the application of shear rate (•), and under flow as a

function of the shear γ, for three volume fractions, φ = 20 % ((a)), φ = 40 % ((b)) and φ = 50 % ((c)). Two shear rate values

are represented :  ̇γ = 3.5.10−3(thin dashed lines) and  ̇γ = 3.5.10−1(thick dashed lines). Continuous lines are exponential fits

of the data, Z(γ) = Z∞ + (Z0 − Z∞)e

−γ/γc

. (d) Critical shear values deduced from the monoexponential fits, at  ̇γ = 3.5.10−3

(•) and  ̇γ = 3.5.10−1(), as a function of the volume fraction, φ.
16

1.00

0.95

0.90

0.85

Ncomp/Next

0.2 0.4 0.6 0.8 1.0

γ

Fig. 7: Evolution of the ratio of contacts in the compression quadrant over the contacts in the extension quadrant, as a function

of shear, for φ = 20 % (circles), φ = 40 % (squares) and φ = 50 % (triangles). Solid lines and filled symbols (resp. dashed lines

and empty symbols) correspond to suspensions that have aged for tw = 9.3 (resp. tw = 1120).
17

2.0

1.5

1.0

0.5

0.0

-0.5

σ12comp

/

σ12

ext

0.0 0.2 0.4 0.6 0.8 1.0

γ

500

400

300

200

100

0

σ12

/∆θ

-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5

θ

(a)

(b)

Fig. 8: (a) Angular distribution of the stress between contacts (stress σ12 localized along contacts that form an angle [θ, θ+∆θ]

with the flow direction (see fig. 7(a)), with ∆θ = π/180, for a suspension whose volume fraction is φ = 50 % and after a shear of

0.05. Dashed lines (resp. continuous lines) correspond to suspensions that have aged for tw = 9.3 (resp. tw = 1120). (b)Evolution

of the ratio of stress stored along contacts in the compression quadrant over the stress stored along contacts in the extension

quadrant, as a function of shear, for φ = 20 % (circles), φ = 40 % (squares) and φ = 50 % (triangles). Continuous lines and

filled symbols (resp. dashed lines and empty symbols) correspond to suspensions that have aged for tw = 9.3 (resp. tw = 1120).

Arrows point at the non-zero limits of the stress ratio, at φ = 40 % and φ = 50 %, and tw = 1120
18

0.01

2

3

4

5

6

0.1

s

0.2 0.4 0.6

g

0.001

2

4

0.01

2

4

0.1

2

s

0.2 0.3 0.4 0.5 0.6 0.7

g

(a) (b)

Fig. 9: Stress as a function of deformation for two different volume fractions, φ = 20 % (a) and φ = 50 % (b). Continuous lines

correspond to an interaction potential between colloids whose minimum value Umin is 25 kBT (continuous lines), and 125 kBT

(dashed lines). The stress curves for Umin = 125kBT have been divided by 4.6 (φ = 20 %) and 3.2 (φ = 50 %) in order to

match the stress curves for the less deep minimum of interaction potential (Umin = 25kBT).