## Abstract

We present data for the rheology of suspensions of monodisperse particles of varying aspect ratio, from oblate to prolate, and covering particle volume fractions *ϕ* from dilute to highly concentrated. Rheology is characterized by fitting the experimental data to the model of Herschel & Bulkley (Herschel & Bulkley 1926 *Kolloid Z.* **39**, 291–300 (doi:10.1007/BF01432034)) yielding three rheometric parameters: consistency *K* (cognate with viscosity); flow index *n* (a measure of shear-thinning); yield stress *τ*_{0}. The consistency *K* of suspensions of particles of arbitrary aspect ratio can be accurately predicted by the model of Maron & Pierce (Maron & Pierce 1956 *J. Colloid Sci.* **11**, 80–95 (doi:10.1016/0095-8522(56)90023-X)) with the maximum packing fraction *ϕ*_{m} as the only fitted parameter. We derive empirical relationships for *ϕ*_{m} and *n* as a function of average particle aspect ratio *r*_{p} and for *τ*_{0} as a function of *ϕ*_{m} and a fitting parameter *τ**. These relationships can be used to predict the rheology of suspensions of prolate particles from measured *ϕ* and *r*_{p}. By recasting our data in terms of the Einstein coefficient, we relate our rheological observations to the underlying particle motions via Jeffery’s (Jeffery 1922 *Proc. R. Soc. Lond. A* **102**, 161–179 (doi:10.1098/rspa.1922.0078)) theory. We extend Jeffery’s work to calculate, numerically, the Einstein coefficient for a suspension of many, initially randomly oriented particles. This provides a physical, microstructural explanation of our observations, including transient oscillations seen during run start-up and changes of rheological regime as *ϕ* increases.

## 1. Introduction

Suspensions of particles in liquids are ubiquitous in nature and industry. The food, cosmetic, plastics, pharmaceutical, oil and mineral separation industries involve a plethora of multiphase substances, many of which are processed as particle suspensions, regardless of how they are finally presented to the consumer. Natural hazards are frequently associated with fluid flows, nearly all of which carry an appreciable fraction of suspended solid particles (e.g. landslides, mudflows, avalanches, magma flows). Civil defence authorities, charged with managing hazard scenarios, require input from scientists regarding the likely behaviour of the expected flow, which is controlled by its rheology. Our ability to control industrial flow processes accurately, and to predict natural ones, depends critically on our knowledge of the rheology of particle suspensions.

The rheology of a material is commonly characterized by a constitutive equation which relates imposed shear stress *τ* to the resultant material strain rate via the material’s physical properties; for a Newtonian liquid, , where *μ* is the Newtonian viscosity. The addition of a dispersed phase to a Newtonian liquid, whether liquid to form an emulsion (e.g. Frankel & Acrivos 1970), gaseous to form a bubble suspension (e.g. Llewellin *et al.*2002*a*,*b*), or solid as considered here, can lead to the introduction of all kinds of non-Newtonian behaviour. It is common in experimental rheometry to describe the rheology in terms of an apparent viscosity measured at a particular stress or strain rate. For a Newtonian fluid, *η*≡*μ*. However for a non-Newtonian fluid, apparent viscosity is a function of strain rate, hence a single measurement of *η* gives a rather poor indication of rheology. For suspensions, the apparent viscosity of the suspension *η*_{s} is often normalized by the viscosity of the suspending liquid *μ*_{0} (for a Newtonian liquid) and reported as the relative apparent viscosity *η*_{r}=*η*_{s}/*μ*_{0}.

The rheology of a particle suspension is a complex function of its physical properties and of processes that occur at the scale of the suspended particles. The most important factors are particle volume fraction *ϕ*, particle shape, interactions between particles, the spatial arrangement of particles and the nature of the bulk flow field. Other factors which are important in certain particle suspensions, but have received less attention, are the size- and shape-distribution of the particles and inter-particle forces (e.g. electroviscous effects). The latter are more relevant in colloidal and aqueous systems (see review by Jeffrey & Acrivos 1976).

Much research has focused on suspensions of spheres, however particle anisometry introduces several additional effects, most notably: (i) the local flow around a non-spherical particle is different from that around a spherical particle, hence the contribution of the particle to suspension viscosity is also different; (ii) non-spherical particles are orientable and their contribution to suspension viscosity depends on their orientation; (iii) particle interactions are strongly influenced by particle shape; in general, at the same particle volume fraction, the degree of interaction among non-spherical particles will be greater than among spherical particles.

Published empirical research on particle suspensions has focused on determining a relationship between (relative) apparent viscosity and particle volume fraction *η*_{r}(*ϕ*). This work suffers from two problems that have been commented upon in previous reviews (Rutgers 1962*a*,*b*; Thomas 1965). First, the focus on apparent viscosity means that non-Newtonian behaviour typically receives only a qualitative description. Second, experimental conditions are not always well constrained: experiments may span a wide range of Peclet numbers and particle Reynolds numbers, even in a single viscosity determination; particle size distributions may be highly polydisperse; suspending fluids may have non-Newtonian rheology; inter-particle forces may be present. As a consequence of these twin problems, results for *η*_{r}(*ϕ*) obtained by different workers may differ by orders of magnitude. Moreover, the link between these empirical observations of viscosity and rheology and the hydrodynamics on the particle-scale, as deduced in theoretical and experimental studies, has not been fully explored, hence a physical explanation for the observations remains elusive. These are the problems that we address in this study.

We first review past research and theoretically explore the link between particle-scale motions and rheology by considering the contribution that an individual particle makes to the viscosity of the suspension. We then present rheological data for suspensions of particles of varying degrees of anisometry and varying particle volume fraction, in a strictly Newtonian suspending liquid. Combining these two strands of investigation allows us to develop constitutive equations that quantify the observed non-Newtonian behaviour of suspensions of spherical and non-spherical particles across a wide range of particle volume fractions, and to explain bulk rheological behaviour in terms of the underlying hydrodynamics and particle motions.

## 2. Suspension viscosity and particle motions

### (a) Suspensions of spherical particles

Rutgers (1962*a*,*b*) and Thomas (1965) present critical summaries of early research which concentrate on the relationship between relative viscosity and particle volume fraction. Both find marked disagreement among the various experimental studies. Rutgers estimates an ‘average’ curve for *η*_{r}(*ϕ*) while Thomas, recognizing that the data scatter is in part due to systematic variations in experimental conditions, processes the data to minimize the scatter before attempting a functional fit to the data. Both curves are presented later in §4*a* where we compare them with our experimental data.

Both Rutgers and Thomas find three regimes of *η*_{r}(*ϕ*): (i) A dilute regime, restricted to (Thomas) or (Rutgers), where *η*_{r}(*ϕ*) is approximately linear and rheology is Newtonian. (ii) A semi-dilute regime, , where *η*_{r} shows a higher order dependence on *ϕ* but behaviour is still approximately Newtonian. (iii) A concentrated regime, starting near *ϕ*=0.25, characterized by rapid growth of apparent viscosity and increasingly non-Newtonian behaviour (particularly shear-thinning) with particle volume fraction.

Suspension behaviour in the dilute limit was first addressed theoretically by Einstein (1906, corrected 1911), who derived an analytical solution for the hydrodynamics around an isolated sphere which yields:
2.1
where the constant *B* is variously referred to as the ‘Einstein coefficient’ or the ‘intrinsic viscosity’, and takes the value *B*=2.5. Despite extensive experimental effort over the last century, this value for the coefficient *B* has not been incontrovertibly validated, with different researchers favouring values covering the range ; (Happel 1957; Brenner 1970; Jeffrey & Acrivos 1976; Pabst *et al.* 2006). Notwithstanding uncertainty over Einstein’s value for *B*, the restriction to very dilute suspensions demonstrated by Rutgers (1962*a*,*b*) and Thomas (1965) severely limits the practical utility of Einstein’s result.

Work in the semi-dilute regime has focused on finding coefficients for the higher order terms in *ϕ* neglected by Einstein, producing relationships of the form:
2.2
with *B*=2.5 and derived from consideration of particle–particle interactions (Guth & Gold 1938; Vand 1948; Manley & Mason 1955). Lower values of *B*_{1} have been found when Brownian motion (Saito 1950; Batchelor 1977) and inertia (de Bruijn 1942) are important.

The polynomial relationships discussed above typically describe experimental data poorly when *ϕ*>0.25 and get worse as particle volume fraction increases. One reason for this is that they predict a finite value of viscosity as . This is clearly unphysical since the densest possible packing for monodisperse spherical particles is *ϕ*_{m}≈0.74, at which particle concentration the viscosity must be infinite. Note that the maximum packing that is achieved in disordered suspensions is typically rather lower: numerical simulations yield a value for random close-packing of monodisperse spheres of *ϕ*_{m}≈0.64 (Rintoul & Torquato 1996). Reported values for sheared suspensions, e.g. *ϕ*_{m}≈0.67 (Rutgers 1962*a*) and *ϕ*_{m}≈0.68 (Kitano *et al.* 1981), indicate that shear can impose additional structure within a suspension compared with random close-packing.

In view of the above discussion, the *η*_{r}(*ϕ*) relationships that most successfully model experimental data in the concentrated regime include *ϕ*_{m} as a parameter. Krieger & Dougherty (1959) considered the contribution of successive packets of suspension to the total particle volume fraction and to the suspension viscosity (termed the ‘functional equation’ approach by Pabst 2004) deriving the following relationship:
2.3
where *B* is the Einstein coefficient. This relationship is a popular choice for fitting to experimental data, in which case *B* and *ϕ*_{m} are considered to be fitting parameters (Jeffrey & Acrivos 1976; Pabst 2004; Pabst *et al.* 2006).

Maron & Pierce (1956) derived an equation with the same functional form as equation (2.3) from consideration of the Ree–Eyring flow equations:
2.4
This relationship can be used for data fitting when only one fitting parameter is desired and has found widespread application. Comparison with equation (2.3) implies *B**ϕ*_{m}=2. Note that binomial expansion of both equations (2.3) and (2.4) recovers the polynomial given in equation (2.2).

Other approaches to dealing with concentrated suspensions have been proposed. Frankel & Acrivos (1967) assume that the viscosity of concentrated suspensions is controlled by energy dissipation in the narrow gaps between the particles and derive a theoretical relationship on those grounds. Hoffman (1974) assumes that, at high concentrations, particles organize themselves into sheets that flow over one another and derives a model for suspension rheology that includes non-Newtonian effects.

Observations of non-Newtonian rheology in suspensions of spheres, evidenced by shear-thinning or non-zero yield stress, have been reported by several workers; for example, Rutgers (1962*a*) reports non-Newtonian behaviour in 24 of 34 experimental datasets reviewed, with typical onset at *ϕ*>0.25 to 0.45. Stickel & Powell (2005) give a comprehensive summary of the types of non-Newtonian behaviour observed in concentrated suspensions of force-free spheres and discuss their origins in terms of suspension microstructure (note that a broader, though less detailed, discussion of non-Newtonian effects that includes anisometric particles and inter-particle forces is given by Jeffrey & Acrivos 1976). Stickel & Powell (2005) focus primarily on the effects of Brownian motion and of inertia that have been recognized as important at least since de Bruijn (1951).

The impact of Brownian motion is governed by the Peclet number (Frankel & Acrivos 1967; Brenner 1974; Stickel & Powell 2005):
2.5
where *k*=1.38×10^{−23} J K^{−1} is the Boltzmann constant, *T* the absolute temperature and *a* the particle radius. Experiments show that viscosity becomes independent of Peclet number for (Stickel & Powell 2005), as Brownian interactions are disrupted by hydrodynamic forces (Krieger 1972; Brenner 1974).

The impact of inertia is governed by the particle Reynolds number (Stickel & Powell 2005):
2.6
where *ρ*_{0} is the density of the suspending fluid. Stickel & Powell (2005) link inertial effects, which become important at , to the shear-thickening rheology which has been observed experimentally by several workers (e.g. Hoffman 1972; Barnes 1989; and references therein).

The Stokes number measures the coupling between the solid and fluid phases, variations in which may lead to non-Newtonian behaviour (Coussot & Ancey 1999):
2.7
where *ρ*_{p} is the particle density and λ the characteristic length of the particles. When *St*≪1, fluid–solid coupling is strong.

Another non-Newtonian phenomenon commonly observed in concentrated suspensions is yield stress (reviewed by Barnes 1999). Heymann *et al.* (2002) present evidence for an elastic component of deformation at high particle concentrations and conclude that, at low shear stress, particles interact as an elastic network, which breaks up when yield stress is reached. They find that yield stress is a function of particle volume fraction and is larger for smaller spheres because of the greater density of particle–particle interactions, which gives rise to a more robust particle network.

### (b) Suspensions of non-spherical particles

The behaviour of non-spherical particles in shearing flow has been well studied for particles that have an axis of rotational symmetry and have fore–aft symmetry along that axis (e.g. prolate and oblate spheroids, cylinders, rods, discs, dumb-bells, etc.). In this work we only consider particles which have approximately this symmetry. Such particles may be characterized by the particle aspect ratio, *r*_{p}=*l*_{a}/*l*_{b}, where *l*_{a} is the length of the particle’s axis of rotational symmetry (the *a*-axis) and *l*_{b} is its maximum diameter perpendicular to that axis. Brenner (1974) showed that any particle with such rotational and fore–aft symmetry can be treated as an equivalent spheroid with aspect ratio *r*_{e}=*a*/*b* where *a* is the length of the semi-axis of rotational symmetry and *b* is the length of the orthogonal semi-axis; the relationship between *r*_{p} and *r*_{e} for a given shape must usually be determined by the experiment. The motion of a spheroidal particle in shearing flow, and of the fluid surrounding it, was solved analytically by Jeffery (1922). In the following sections, we review this work and deduce first the contribution to viscosity made by an individual particle as a function of its orientation in the flow and then the resultant viscosity of a suspension of randomly oriented particles. The equivalence noted by Brenner (1974) means that the analysis is applicable to all the particles considered in this study, and to any particle with both rotational and fore–aft symmetry.

#### (i) Motion of a non-spherical particle in a shearing fluid

Particles suspended in a viscous fluid undergoing simple shear rotate about the flow’s vorticity vector. The particles’ rotation accommodates some of the bulk strain experienced by the suspension, hence a particle that is rotating increases the suspension viscosity by a smaller amount than a particle that cannot rotate: Brenner (1970) showed that the Einstein coefficient (equation (2.1)) for a dilute suspension of spherical particles is *B*=4 for non-rotating particles compared with Einstein’s value of *B*=2.5 for particles that are free to rotate.

Figure 1 shows a prolate spheroid in simple shearing flow and defines the coordinate system and angles which describe the particle’s motion. It is instructive to consider two special cases of particle motion: (i) If the particle is aligned so that its *a*-axis is parallel to the vorticity vector (*θ*=0), the particle rotates about that axis (measured by *ψ*), the motion is steady with a constant rotation rate , and angle *φ* is undefined. (ii) If the particle is aligned so that its *a*-axis is perpendicular to the vorticity vector (*θ*=*π*/2) the particle tumbles ‘end over end’. The orientation of the particle is given by *φ* which is equal to zero when the *a*-axis is parallel to the velocity gradient (points along the *y*-axis in figure 1). The rotation rate of the particle about the vorticity vector is periodic and varies with particle orientation; for a prolate particle, the rotation rate is maximal when *φ*=0,*π* and minimal when *φ*=*π*/2,3*π*/2. In any other orientation (0<*θ*<*π*/2), the particle’s *a*-axis describes an elliptical cone about the vorticity vector; the motion of the particle’s *a*-axis is periodic and is sometimes called a ‘Jeffery’s orbit’.

Jeffery (1922) determined the motion of a spheroidal particle in a dilute, shearing suspension. For simple shearing flow, the tumbling rate is given by:
2.8
Integration of this equation gives the period of rotation of the particle about the vorticity vector:
2.9
Note that tumbling rate and period are independent of *θ*, which is given by:
2.10
where *C* is the ‘orbit constant’, which can be uniquely determined from equation (2.10) if *θ* and *φ* are known simultaneously at any point (e.g. initial orientation). The likelihood of finding a particle in a particular orientation is inversely proportional to the particle’s angular velocity in that orientation. The *φ*-orientation probability density function is therefore given by , where *T* enters as the normalizing constant to ensure that the integral over all possible orientations is 1.

Figure 2 shows the angular velocity and *φ*-orientation probability density function of oblate (*r*_{e}<1) and prolate (*r*_{e}>1) spheroidal particles in a simple shearing flow as a function of particle orientation for particles with aspect ratios in the range 0.1≤*r*_{e}≤10. In general, non-spherical particles rotate rapidly when they project a large profile across the flow (for a prolate particle, this is maximal when *φ*=0; for an oblate particle, this is maximal when *φ*=*π*/2) and rotate slowly when they present a small profile across the flow (*φ*=*π*/2 for prolate and *φ*=0 for oblate particles). Consequently, the particles spend proportionally more of their rotational period ‘aligned’ with the flow than they do ‘across’ the flow. This effect increases as aspect ratio deviates further from unity. Jeffery’s equations for the motion of particles have been validated numerically (Ježek *et al.* 1999; Jiang 2007) and experimentally (Trevelyan & Mason 1951; Goldsmith & Mason 1962; Anczurowski & Mason 1968). However, Yarin *et al.* (2000) show that, in the general case of ellipsoidal particles not having rotational symmetry, chaotic motions are predicted for high aspect ratio particles.

#### (ii) Particle orientation and viscosity

The rate of work that must be done to shear a volume of material *V* is related to the material’s apparent viscosity *η* by (e.g. Happel 1957):
2.11
If that material is a packet of fluid suspending a particle, then the total rate of doing work is the sum of the rate of work that must be done to shear the equivalent volume of particle-free fluid and the extra rate of work owing to the presence of the particle . Equation (2.11) can be rewritten to give the suspension’s apparent viscosity as a sum of these two contributions:
2.12
For a Newtonian suspending fluid, the first term on the right-hand side is identically the fluid’s Newtonian viscosity *μ*_{0}, hence we can rewrite equation (2.12) to give the relative viscosity of the suspension:
2.13
where *V*_{p} is the volume of the particle and *ϕ*=*V*_{p}/*V* is the particle volume fraction. Note that equation (2.13) has the form of the Einstein equation (2.1) with .

Jeffery (1922) derives an equation for (his eqn 61) from which, via our equation (2.13), we can deduce a particle’s instantaneous contribution to the viscosity of a dilute suspension for any particle orbit. For orbit *C*=0, the particle’s axis of symmetry is aligned with the vorticity vector; the particle ‘rolls’ along in the flow with steady motion hence the particle’s contribution to viscosity is constant with strain. For any other orbit, a particle’s contribution to viscosity varies with angle *φ*, hence also with strain. The variation is maximal when the particle is in orbit , i.e. when the particle is tumbling in the flow with its axis of rotational symmetry perpendicular to the vorticity vector; this case is plotted in figure 3. When the particle’s longest dimension is parallel to the flow (*φ*=*π*/2 for a prolate particle and *φ*=0 for an oblate particle), flowline distortion is minimal and the fluid must do little extra work to flow around the particle, hence *B* is minimal. When the particle is oriented with its longest dimension perpendicular to the flow, the particle is rotating most rapidly—moving at almost the same speed as the adjacent fluid—consequently flowline distortion is, again, minimal and *B* is minimal. When the particle is at *π*/4 to the flow, it presents a large profile to the flow but is moving more slowly than the flow, consequently flowline distortion is maximal, hence *B* is maximal. Figure 3*b*,*d* show the contribution to viscosity with strain (equivalently, time, for constant strain rate) over a complete revolution. Note that, since particles spend most time roughly aligned with flow, the time-averaged contribution to viscosity is nearer the minimal value than the maximal value.

#### (iii) Viscosity of suspensions of non-spherical particles

Jeffery (1922) determines the time-averaged Einstein coefficient *B* by integrating the work that must be done to rotate a particle over a complete orbit (his eqn 62). Jeffery tabulates *B* for a range of aspect ratios for two limiting orbits: *C*=0 and (*a*-axis, respectively, parallel to and perpendicular to the vorticity vector; see figure 1 and equation (2.10)). In fact, Jeffery’s eqn (62) can be used (via our equation (2.13)) to determine the exact, time-averaged value of *B* for any combination of aspect ratio *r*_{e} and orbit constant *C*, noting that his constant of integration *k*=*b*/*C*. Figure 4*a* presents a suite of *B*(*r*_{e}) curves for various orbit constants.

We have extended this approach to determine the viscosity of a suspension of particles with a distribution of orientations. We calculate the time-averaged value of *B* by numerical integration of Jeffery’s eqn (61), over a complete rotational period, for one million identical particles with random initial orientations, in the high *Pe* limit (no Brownian diffusion of orientations). In the dilute limit particles do not interact, hence, their individual contributions to viscosity can be averaged to determine the suspension viscosity. We repeated this calculation for suspensions of particles with aspect ratios in the range 0.01≤*r*_{e}≤100; results are plotted in figure 4*b*. The results are in good agreement with values of *B* calculated from the equations and tables of goniometric factors (high *Pe* limit) presented in Brenner (1974) except at extreme values of *r*_{e}, where the results are sensitive to small errors in the tabulated values of goniometric factors.

The viscosity of a suspension of non-spherical particles depends on hydrodynamic and particle–particle interactions as well as particle orientation. We follow Doi & Edwards (1978) who define three regimes, given *N* particles of length *L* and diameter *D* in unit volume: dilute, where interactions are rare, *N*≪1/*L*^{3}; semi-dilute, where hydrodynamic interactions are important, 1/*L*^{3}<*N*<1/*L*^{2}*D*; concentrated, where particle–particle contact is important, *N*>1/*L*^{2}*D*. For a spheroidal particle with aspect ratio *r*_{e}=10, the regime boundaries are *ϕ*=0.005 and *ϕ*=0.05.

In the dilute regime, Jeffery theorized that particle orientations would drift over time to the orientation that minimizes the work done to rotate the particle: *C*=0 for prolate and for oblate particles. This configuration corresponds to the minimum line in figure 4*b*. This hypothesis has been confirmed experimentally for the dilute regime by Anczurowski & Mason (1967). Two datapoints determined from their table 6 are plotted in figure 4*b* and show that the viscosity of dilute suspensions at large strains is lower than would be expected for randomly oriented particles. Behaviour in the semi-dilute regime has been considered experimentally by Stover *et al.* (1992). They find that the distribution of particle orbits remains close to random over time indicating that hydrodynamic interactions between particles lead to rotary diffusion of particle orientations that counters the drift towards low *B* predicted by Jeffery. These results are supported by the numerical study of Rahnama *et al.* (1995). Ferec *et al.* (2009) investigated behaviour in the concentrated regime numerically, and found that particle–particle contact interactions push the viscosity higher than predicted for randomly oriented particles.

Experiments involving particle volume fractions in the concentrated regime (e.g. Powell 1991; Djalili-Moghaddam & Toll 2006; Pabst *et al.* 2006) show that the relationship between apparent viscosity and particle volume fraction is linear only for the lowest concentrations measured, hence, in general, the Einstein relationship (equation (2.1)) is not valid. The *η*_{r}(*ϕ*) relationship is also strongly dependent on *r*_{e}; the dependence of *η*_{r} on *ϕ* is stronger as *r*_{e} increases for prolate particles (and, presumably, as *r*_{e} decreases for oblate particles). Pabst *et al.* (2006) use the Krieger–Dougherty relationship (equation (2.3)) to fit their data. This provides a link to the ‘dilute’ behaviour through the Einstein coefficient that appears in the exponent; they find values of *B*=4.99 and *B*=10.2 for *r*_{e}=5 and *r*_{e}=16, respectively. These datapoints plot significantly higher than the maximum value allowable by Jeffery’s analysis (figure 4*b*). This is to be expected since the values of the Einstein coefficient are derived from a dataset that is entirely in the concentrated regime where hydrodynamic and collision interactions between particles—assumed negligible in Jeffery’s analysis—are important.

## 3. Experiments

### (a) Materials and techniques

Measurements were performed on suspensions of particles with five different geometries (figure 5 and table 1) in silicone oil (Cannon Viscosity Standard N15000: 41.32 Pa s at 25^{°}C). This non-aqueous, viscous liquid was chosen to prevent electroviscous effects and minimize particle settling during the measurements. Particle dimensions were determined optically from images taken with a Zeiss SteREO V.8 stereomicroscope and using image analysis software JMICROVISION. Densities were determined via He-pycnometry. Oblate particles (figure 5*a*) are ultra-fine, polyacrylic art glitter by Efco. Spherical particles (figure 5*b*) are soda-lime silica glass beads from Potters-Ballotini ltd. Prolate particles of three different aspect ratios were used: silicon carbide grit (prolate A, figure 5*c*) and two types of acicular wollastonite crystals (prolate B, figure 5*d*; and prolate C, figure 5*e*) obtained from Osthoff-Petrasch, Germany.

Flow curves, , were determined using a ThermoHaake MARS II rotational rheometer with a parallel plate sensor system (diameter: 2*R*=35 mm, gap size: *h*=1.5 mm, temperature: 25^{°}C) in controlled-stress mode (figure 6*a*). A 20-step ‘up ramp’ of incrementally increasing shear stress was followed by a 20-step ‘down ramp’. At each step the rheometer records the strain rate when equilibrium stress conditions are reached. In all parallel plate systems, the strain rate depends on the radial position and it is conventional to use the strain rate at the outer margin of the upper sensor. The maximum stress value was usually set at 500 Pa, except for some highly viscous samples, where a lower stress was used in order to prevent wall-slippage. Maximum strain rates achieved were approximately 0.03–0.05 s^{−1} for concentrated samples, and up to 12 s^{−1} for very dilute suspensions.

### (b) Transient effects at flow initiation and pre-shear treatment

Suspensions of non-spherical particles must be sheared to large total strain before repeatable rheometric results can be obtained (Hinch & Leal 1973; Okagawa *et al.* 1973; Ivanov *et al.* 1982; Powell 1991). This ‘pre-shear treatment’ is required because particles are, initially, randomly oriented and equilibrium orientation distributions are established only after a period of shearing. Before equilibrium is achieved, transient oscillations in apparent viscosity are observed (figure 7*a*).

The oscillations arise because the motion of a non-spherical particle in simple-shearing flow is periodic and its contribution to suspension viscosity is also periodic (see §2*b*(i), (ii) and figure 3). The contribution to viscosity of a randomly oriented particle will decay as it rotates to become aligned with flow, then spike periodically as it flips over. The rotational period of a particle depends only on its aspect ratio and the local strain rate (equation (2.9)) so, for identical particles in a uniformly shearing flow and in the absence of particle interactions, all of the particles rotate back to their initial orientation with the same period, leading to viscosity oscillations. In our experiments, the suspended particles have a range of aspect ratios, hence they also have a range of periods of rotation. As strain increases, particle flipping becomes asynchronous causing the viscosity oscillations to decay. Variations in strain rate throughout the flow enhance asynchrony of rotation and the decay of oscillations. We find that oscillations are damped for all suspensions for , defining a minimum required strain for the ‘pre-shear treatment’.

Figure 7*b* shows the result of numerical determination of the expected viscosity of each suspension using our equation (2.13) and equation 61 of Jeffery (1922) following the approach we outline in §2*b*(iii). The motion of 100 000 particles with random initial orientations was calculated as a function of strain using Jeffery’s equations for the motion of ellipsoidal particles. Particle aspect ratios were distributed according to the means and standard deviations presented in table 1 and the strain-rate distribution owing to the geometry of the sensor system was taken into account. After each small increment of strain the instantaneous contribution to viscosity of each particle was calculated. In the dilute limit, particles do not interact, hence their individual contributions to viscosity can be averaged to determine the viscosity of the suspension. The experimental and numerical data in figure 7 show some qualitative agreement. In both datasets, the oscillations in the viscosity of suspensions of prolate particles—where aspect ratio has a large standard deviation leading to a commensurately wide range of particle rotational periods—show the greatest amplitude, but decay quickly. For oblate particles, which have a uniform aspect ratio, oscillations have small amplitude but decay slowly. The experimentally determined viscosity of the suspension of spheres is roughly constant with strain, consistent with theory.

The differences between the experimental and theoretical datasets are also informative. Firstly, the measured viscosities are higher than the numerical values; this point is returned to later in §5*a*. Secondly, the oscillations decay more rapidly in the experimental data. Okagawa *et al.* (1973) predict that particle interactions should lead to a more rapid approach to equilibrium by perturbing particles’ rotational periods. This suggests that particle interactions are important in our experiments even at *ϕ*≈0.05.

### (c) Rheological characterization

The shape of the flow curves (§3*a*) shows a systematic dependence on particle volume fraction. At low concentrations ( for spherical particles) the flow curve is a straight line through the origin, indicating Newtonian rheology (figure 6*b*). For higher particle concentrations the flow curve becomes convex upwards, indicating shear-thinning rheology (figure 6*c*). At the highest particle concentrations ( for spherical particles) the flow curve has a non-zero intercept on the stress axis, indicating a yield stress (figure 6*d*).

It has been common throughout previous experimental studies of suspension rheology to present only the apparent viscosity of the suspension . This one-dimensional measure is insufficient to capture the non-Newtonian rheology shown in figure 6*c*,*d*. As our first data-processing step, we fit the three parameter Herschel–Bulkley model (Herschel & Bulkley 1926) to each raw flow curve:
3.1
where *τ*_{0} is the yield stress below which there is no flow, *K* is the consistency (which takes the same value as *η* evaluated at ) and *n* is the flow index which defines the degree of non-Newtonian behaviour (shear thickening for *n*>1, and shear-thinning for *n*<1). Characterizing a liquid in terms of these three parameters gives a much richer description of its rheology than is possible with a single value of apparent viscosity.

The precision of the measurements is estimated at ±2 per cent, and includes instrumental resolution, uncertainty in filling the sensor, and the fact that, for very prolate particles, the ratio of gap width to particle axis length can be less than 10. In addition, the fit of the Herschel–Bulkley equation (3.1) to the experimental data leads to a standard error in the mean value of each of the rheological parameters *τ*_{0}, *K* and *n*. The combination of these uncertainties are given as *y*-error bars in figures 8 and 9. Uncertainty in the value of the particle concentration *ϕ* determined for a sample may arise both from particle inhomogeneities (polydisperse sizes and aspect ratios) and uncertainties during the sample preparation procedure (i.e. weighing of the phases, and mixing of the suspension). These are estimated to be ±3 per cent and are displayed as *x*-error bars in figures 8 and 9.

## 4. Rheology of suspensions of spherical particles

The yield stress *τ*_{0}, consistency *K* and flow index *n* (equation (3.1)) were determined for 24 suspensions of monodisperse spherical particles with particle volume fractions in the range 0.005≤*ϕ*≤0.59. Results are plotted in figure 8.

### (a) Consistency of suspensions of spheres

The consistency *K* is the parameter most closely related to the apparent viscosity that is typically quoted in experimental studies of suspension rheology; for a Newtonian suspension, where *n*=1 and *τ*_{0}=0, *K*≡*η* and has the dimensions of dynamic viscosity: Pa s. For the case of shear-thinning, where *n*≠1, *K* has non-integer dimensions Pa s^{n}. By analogy with the definition of *η*_{r}, we define the relative consistency *K*_{r}=*K*/*μ*_{0}.

Figure 8*a* demonstrates that suspension consistency shows a systematic variation with particle volume fraction and compares our *K*_{r}(*ϕ*) data with the *η*_{r}(*ϕ*) results given by Rutgers (1962*a*) and Thomas (1965). In both cases, the data are in excellent agreement for ; the agreement with Rutgers’ data is excellent for . For our data fall between the two curves. The comparison of *K*_{r}(*ϕ*) with *η*_{r}(*ϕ*) breaks down for higher *ϕ* as the flow index *n* deviates from 1. Figure 8*a* also shows the best fit of the Maron–Pierce model (equation (2.4)) to our data, with *ϕ*_{m} as an adjustable parameter. The fitting process, which was done on to avoid biasing the fit to large values of *K*_{r}, was performed using the freely available statistical analysis package SimFit (http://www.simfit.man.ac.uk) and yielded a value of *ϕ*_{m}=0.633. The fit is excellent across the dataset (*R*^{2}=0.997) and is within the error bars of all but three of the 24 datapoints. The two-parameter Krieger–Dougherty model (equation (2.3)) was also fitted to the data, yielding values of *B*=3.27 and *ϕ*_{m}=0.641 (*R*^{2}=0.998). With these parameter values, the exponent of equation (2.3) has the value *B**ϕ*_{m}=2.10, which agrees well with the exponent of 2 in the Maron–Pierce model. The Maron–Pierce and Krieger–Dougherty best-fit curves are almost indistinguishable. We note that the values of *ϕ*_{m} we determine are very close to the value of *ϕ*_{m}≈0.64 determined for the random close-packing of spheres by Rintoul & Torquato (1996) and only slightly lower than the value of *ϕ*_{m}≈0.67 identified by Rutgers for a sheared suspension.

Figure 8*b* shows the same data, but focuses on the semi-dilute and dilute (inset) regimes. In the semi-dilute regime, the polynomial presented in equation (2.2) (to second order in *ϕ*) shows excellent agreement with our data for when the parameter values *B*=2.5,*B*_{1}=14.1 determined by Guth & Gold (1938) are used. In the dilute regime both the Maron–Pierce and the Krieger–Dougherty relationship slightly overpredict the consistency, but are within experimental errors. The agreement between the Einstein relationship and the data is reasonable for ; a rather higher concentration than the values of 0.01 and 0.02 suggested by Thomas (1965) and Rutgers (1962*a*), respectively.

### (b) Shear-thinning in suspensions of spheres

We observe shear-thinning for moderate particle volume fractions (figure 8*c*); it becomes progressively more important with increasing concentration for . In our experiments, shear-thinning must result from a decrease in viscosity as strain rate increases rather than from a dependence of viscosity on total strain (thixotropy) since the data from the ‘up-ramp’ (increasing ) and ‘down-ramp’ (decreasing ) approximately coincide, while total strain increases throughout the experiment (figure 6*c*,*d*). The dependence of shear-thinning on *ϕ* suggests that its origin is in the micro-mechanics of interactions between particles.

Previous theoretical and experimental studies have attributed shear-thinning to the dependence of viscosity on the suspension’s Peclet number, particle Reynolds number or Stokes number (see §2*a*). Stickel & Powell (2005) indicate Newtonian behaviour for values of and . Using equation (2.5) with *μ*_{0}=42 Pa s, *a*=90 μm, s^{−1} and *T*=300 K, we obtain *Pe*=1.4×10^{8} as a minimum for our experiments; consequently all our experiments are conducted in the hydrodynamic regime, where Brownian motion is negligible. Similarly, the maximum values of the particle Reynolds and Stokes numbers attained can be calculated from equations (2.6) and (2.7). Assuming *ρ*_{0}=894 kg m^{−3}, *ρ*_{p}=2500 kg m^{−3}, *a*=135 μm, s^{−1}, *μ*_{0}=42 Pa s and λ≈*a* we obtain *Re*=7.7×10^{−6} and *St*=2.2×10^{−5}, respectively; consequently, all of our experiments are conducted in the viscous regime where particle inertia is negligible and where particle and fluid motions are fully coupled.

At least one other experimental study (Chong *et al.* 1971) has also observed shear-thinning in suspensions of spheres in the regime of *Pe*≫10^{3} and *Re*≪10^{−3}. In both that study and this, Brownian and inertial effects are negligible, furthermore, the particles are not charged and are too large for van der Waals forces to be significant. Consequently, it is difficult to imagine a mechanism through which the micro-structure of the suspension could be controlled by strain rate, hence we conclude that the dependence of viscosity on strain rate has another origin.

We hypothesize that highly localized viscous heating of the suspending fluid could be responsible. In shearing flow, the rate of heat dissipation per unit volume (power density) in a packet of shearing fluid depends strongly on the local strain rate (equation (2.11)). In general the strain rate in the small gap between two particles in a shearing suspension will be much higher than the bulk strain rate, leading to a higher power density. As the gap between two particles decreases, this will reach a point when diffusion of heat out of the packet of fluid will be insufficient to prevent a rise in the temperature of the packet. This will lead to a decrease in the viscosity of the liquid in the gap, effectively lubricating the passage of the particles past one another. Since the viscosity of concentrated particle suspensions is dominated by the viscous forces in these small gaps (Frankel & Acrivos 1967), this mechanism could exert a strong control on suspension viscosity. Heating is enhanced as bulk strain rate increases leading to shear-thinning rheology. Note that, although this model requires that the energy dissipated per unit volume of fluid in the packet is high, the volume of the packet is very small so that the rise in the bulk temperature of the suspension may be negligible.

### (c) Yield stress of suspensions of spheres

At the highest particle volume fractions (), the suspensions of spheres show a yield stress (figure 8*d*). This phenomenon was investigated in detail by Heymann *et al.* (2002) who attributed the yield stress to the formation of networks of particles. They found that the network deforms elastically in response to an applied stress, which is transmitted through the network via direct particle–particle contact. The yield stress is reached when the applied stress is sufficient to cause the network to break up, beyond which point the suspension flows viscously. Heymann *et al.* find that the yield stress increases with increasing particle volume fraction and that their data for *τ*_{0}(*ϕ*) are well-described by a modified form of the Maron–Pierce relationship:
4.1
where *τ** is a fitting parameter; physically it is the value of the yield stress at and is related to the size of the spheres. Figure 8*d* shows the best fit of equation (4.1) to our data.

A reasonable fit (*R*^{2}=0.977) is obtained when we use the value of *ϕ*_{m}=0.633 determined in §4 *a* by fitting the Maron–Pierce relationship (equation (2.4)) to the cognate data for relative consistency *K*_{r}(*ϕ*). This gives *τ**=0.153 Pa. If, instead, *ϕ*_{m} is treated as a fitting parameter, then an excellent fit (*R*^{2}=0.999) is obtained for *ϕ*_{m}=0.611 and *τ**=4.83×10^{−2} Pa. Whichever value of *τ** is used, it fits with the trend of decreasing *τ** with increasing particle size noted by Heymann *et al.*; we present our data and theirs in table 2. Further experiments with particles of different sizes should allow a predictive *τ**(*a*) relationship to be determined.

## 5. Rheology of suspensions of particles of arbitrary shape

Experimental data for relative consistency *K*_{r}, yield stress *τ*_{0} and flow index *n* for all particles (see §3*a* and table 1) are presented in figure 9. Raw flow curves were generated and processed following the methodology presented in §3. The data for suspensions of particles of all shapes follow the same trends in *K*_{r}, *τ*_{0} and *n* that are identified and discussed for suspensions of spherical particles in §4. The chief difference between the data for suspensions of spherical and non-spherical particles is that *K*_{r}, *τ*_{0} and *n* all show a stronger dependence on *ϕ* as particle anisometry increases, i.e. *K*_{r} increases more-rapidly with *ϕ*, and yield stress and shear-thinning develop at lower *ϕ*, as anisometry increases.

Figure 9*a* shows that the Maron–Pierce relationship (equation (2.4)) provides an excellent fit to the *K*_{r}(*ϕ*) data for all particles. The values of *ϕ*_{m} which give the best fit for each type of particle are presented in table 3, which also presents the best fit of the two-parameter Krieger–Dougherty relationship (equation (2.3)). There is a clear trend in the maximum packing fraction with the aspect ratio of the suspended particles, which is explored further in §5 *a* below. Figure 9*b*,*d*,*f* show the same data presented in figure 9*a*,*c*,*e* but with the particle volume fraction normalized by the value for *ϕ*_{m} determined from the Maron–Pierce fit for the appropriate particle aspect ratio. When normalized in this way, the *K*_{r}(*ϕ*/*ϕ*_{m}) data collapse to a single curve. The *τ*_{0}(*ϕ*/*ϕ*_{m}) curves also collapse fairly well, but the *n*(*ϕ*/*ϕ*_{m}) curves do not. These three datasets are discussed in more detail below.

### (a) Consistency of suspensions of particles of arbitrary aspect ratio

The aspect ratio of the particles in a suspension has a strong influence on the consistency of the suspension in the case of prolate particles. The curves in figure 9*a* show that, for any given particle volume fraction, the relative consistency is lowest for spherical particles and increases with increasing *r*_{p} for prolate particles. The single dataset for oblate particles indicates that their influence on consistency is similar to that of spherical particles.

The trend in each of the datasets in figure 9*a* is strikingly similar and, consistent with the findings of Kitano *et al.* (1981) and Pabst *et al.* (2006), we find that the curves are well-described by the Maron–Pierce relationship. Table 3 gives the best-fit values of *ϕ*_{m} alongside *r*_{p} values for each suspension. Physically, the quantity *ϕ*_{m} is the particle volume fraction at which there is no longer sufficient ‘free space’ in the suspension for the particles to be able to move past one another, hence, the suspension becomes ‘jammed’ and its viscosity becomes infinite. The data indicate that suspensions of particles with high aspect ratio become jammed at lower particle volume fractions than suspensions of more equant particles. This is physically intuitive since the degree of interaction between elongate particles is higher than for equant particles. The collapse of the *K*_{r}(*ϕ*/*ϕ*_{m}) data to a single curve (figure 9*b*) indicates that the dependence of *K*_{r} on particle aspect ratio results from the dependence of *ϕ*_{m} on particle aspect ratio.

The dependence of the maximum packing fraction, determined from the Maron–Pierce fits, on particle aspect ratio is plotted in figure 10*a*. The high quality of fit for the Maron–Pierce relationship in figure 9*b* indicates that the product *B**ϕ*_{m}≈2 is valid across a wide range of particle aspect ratios, providing a link with the exact analyses for the dilute limit of Einstein (1906, 1911) and Jeffery (1922) (both discussed in §2) and allowing us to recast the *ϕ*_{m}(*r*_{p}) data as *B*(*r*_{p}) (figure 10*b*).

Previously, Kitano *et al.* (1981) and Pabst *et al.* (2006) have proposed models for *ϕ*_{m}(*r*_{p}) based on their experimental data; both models and datasets are included in figure 10; note that their models are linear in *r*_{p} and predict negative values for *ϕ*_{m} for large *r*_{p}. For comparison, the numerically determined value of *B*(*r*_{p}) is also plotted (assuming equilibrium conditions have been reached for a monodisperse suspension of initially randomly oriented, ellipsoidal particles with *r*_{e}=*r*_{p}; see §2*b*(iii) and figure 4*b*). In the light of the discussion at the end of §2*b*(iii), it is not surprising that the experimental data all plot higher than the numerical curve in figure 10*b* since the suspensions of high-aspect ratio particles used in our experiments cannot be considered ‘dilute’ (*sensu*Doi & Edwards 1978) even for the lowest particle volume fractions that we investigate.

Transforming our *ϕ*_{m}(*r*_{p}) data to *B*(*r*_{p}) reveals a relationship that is linear to a good approximation and is valid for :
5.1
This relationship is based on a linear regression (*R*^{2}=0.986) of the data excluding the datapoint for oblate particles. The data and best-fit line plot between those of Kitano *et al.* (1981) and Pabst *et al.* (2006). Note that Kitano *et al.* used a non-Newtonian suspending liquid and the *r*_{p}=1 experiments of Pabst *et al.* were conducted using suspensions of starch globules which we consider to be unlikely to behave as hard spheres. The simple empirical model that we propose in equation (5.1) is perhaps more useful than those proposed by Kitano *et al.* and Pabst *et al.* since *ϕ*_{m} is nonlinear in *r*_{p} and is non-negative for all aspect ratios.

### (b) Yield stress of suspensions of particles of arbitrary aspect ratio

The plots of *τ*_{0}(*ϕ*) presented in figure 9*c* show that suspensions of particles of all aspect ratios exhibit a negligible yield stress over low and intermediate particle volume fractions, rising steeply at high particle volume fractions. The particle volume fraction at which appreciable yield stress is developed varies systematically with particle aspect ratio. Suspensions of high aspect ratio, prolate particles develop a yield stress at the lowest concentration, and suspensions of spherical particles at the highest concentration. Presumably this trend exists because elongate particles interact with one another to a greater extent at lower particle volume fractions than spherical particles do. Consequently, particle networks capable of accommodating stress elastically also form at lower particle volume fractions for elongate particles.

Figure 9*d* plots yield stress against normalized particle volume fraction *ϕ*/*ϕ*_{m} where *ϕ*_{m} is taken from the Maron–Pierce fits to the *K*_{r}(*ϕ*) data (table 3). The data for the suspensions of particles of different aspect ratio collapse reasonably well to a single curve. In all cases, appreciable yield stress is only developed for (shaded in figure 9*d*). The collapse of the data when normalized to *ϕ*_{m} indicates that, as for the consistency discussed in the previous section, the dependence of yield stress on particle aspect ratio, particularly the onset of appreciable yield stress, is accounted for by the relationship between particle aspect ratio and *ϕ*_{m}.

Following the analytical approach outlined in §4*c*, we fit the *τ*_{0}(*ϕ*/*ϕ*_{m}) relationship (equation (4.1)) proposed by Heymann *et al.* (2002) to the normalized data with *τ** as the only adjustable parameter and using the *ϕ*_{m} values from the Maron–Pierce fits to the *K*_{r}(*ϕ*) data (table 3). The fits obtained, presented in table 4, are good for suspensions of prolate and spherical particles, but rather poor for the suspension of oblate particles. There is no clear systematic trend in the values of *τ** determined from the fit with particle aspect ratio. Heymann *et al.* (2002) suggest that, for suspensions of spherical particles, there is an inverse relationship between *τ** and particle size; i.e. yield stress is higher for suspensions of smaller particles owing to the greater number of particle–particle contacts per unit volume. This was supported by our data for spherical particles presented in table 2. Cross-referencing tables 1 and 4 shows that no such simple trend is evidenced for suspensions of non-spherical particles. This may be owing to variations in inter-particle friction for different particle materials.

### (c) Shear-thinning of suspensions of particles of arbitrary aspect ratio

The onset of shear-thinning behaviour can be observed for suspensions of particles of all aspect ratios at intermediate *ϕ* and becomes more pronounced as *ϕ* increases (figure 9*e*). This non-Newtonian behaviour also shows a systematic dependence on particle aspect ratio: it is weakest for suspensions of spherical particles and strongest for high aspect ratio, prolate particles. This echoes the trends in consistency and yield stress with aspect ratio discussed in the previous sections. As before, for each *n*(*ϕ*) dataset we normalize the particle volume fraction by the maximum packing fraction determined from the best fit of the Maron–Pierce relationship (equation (2.4)) to the *K*_{r}(*ϕ*) data (table 3). The results are plotted in figure 9*f*. In this case it is clear that the normalization is insufficient to collapse the data to a single curve. This indicates that shear-thinning has a stronger dependence on particle aspect ratio than consistency and yield stress do, since the dependence of shear-thinning on aspect ratio is not fully accounted for by the dependence of aspect ratio on *ϕ*_{m}.

The increase in the importance of shear-thinning with increasing particle volume fraction is evidenced by a decrease in *n* with increasing *ϕ*/*ϕ*_{m} in figure 9*f*. For suspensions of all particle aspect ratios, this trend reverses for the most concentrated suspensions. Comparison of figure 9*f*,*d* shows that this reversal corresponds to the onset of yield stress for each suspension. This suggests that there is a link between the microstructural origins of yield stress and shear-thinning but it is not clear what that link is. The region where yield stress appears to influence shear-thinning is shaded in figure 9*f*. In the region where *n*(*ϕ*/*ϕ*_{m}) shows a systematic dependence on particle aspect ratio, corresponding to , we find that the following, purely empirical relationship provides a good fit to the data (*R*^{2}=0.944):
5.2
The values of the parameters in this relationship (the factor of 0.2 and exponent of 4) were determined by simultaneous fit of equation (5.2) to the *n*(*ϕ*/*ϕ*_{m}) datasets for the suspensions of spherical and prolate particles (but not oblate particles). Each dataset was first filtered to remove the datapoints where appreciable yield stress was measured. The filtering retained all data for *ϕ*/*ϕ*_{m}≤0.7 for the suspensions with the most prolate particles (*r*_{p}=9.17) and all data for *ϕ*/*ϕ*_{m}≤0.8 for all other suspensions. Curves of equation (5.2) are plotted for each dataset, up to the appropriate *ϕ*/*ϕ*_{m} limit, in figure 9*f*.

In §4*b* we proposed that shear-thinning may result from localized viscous heating and the consequent decrease in viscosity of the fluid in the small gap between approaching particles. Particles of high aspect ratio rotate more rapidly as they flip over in shearing flow than spherical particles (§2*b*(i)); hence, their approach velocity is likely to be higher, leading to greater viscous heating. Furthermore, while spherical particles will slide past one another fairly rapidly, elongate particles may be in contact with one another for much longer. These microstructural processes might explain why shear-thinning has a stronger dependence on particle aspect ratio than consistency and yield stress do.

## 6. Conclusions

The results presented above provide a suite of equations that can be used to predict the rheology of suspensions of monodisperse particles of varying aspect ratio over particle volume fractions from dilute to highly concentrated. We find that the model of Herschel & Bulkley (1926; our eqn 3.1) provides an excellent basis for a rich rheological characterization via its three parameters—consistency *K* (cognate with viscosity), flow index *n* (a measure of shear-thinning) and yield stress *τ*_{0}. The consistency *K* of particles of arbitrary aspect ratio can be accurately predicted by the model of Maron & Pierce (1956; our equation (2.4)) with a single fitting parameter *ϕ*_{m}. We propose empirical relationships for *ϕ*_{m} (equation (5.1)) and the flow index *n* (equation (5.2)) as a function of average particle aspect ratio *r*_{p}. Appreciable yield stress is seen in all suspensions for . The model developed by Heymann *et al.* (2002; our equation (4.1)) delivers an excellent fit to the data for suspensions of spherical and prolate particles with *τ** taking a value in the range Pa, but showing no apparent systematic dependence on particle aspect ratio. The value of *τ** is likely to be related to the density of particle–particle contacts and the coefficient of friction between the particles. The fit for oblate particle is less good and involves a higher *τ**.

Recasting our data in terms of the Einstein coefficient *B*=2/*ϕ*_{m} enables us to relate our rheological observations to the underlying particle motions during shearing. We extend Jeffery’s (1922) theory to calculate the Einstein coefficient for a suspension of many, initially randomly oriented particles. This provides a physical explanation for both the transient effects observed at the start of each run (figure 7) and the steady-state results (figure 10). We attribute discrepancies between the data and theory to the fact that our suspensions violate the dilute assumption of the theory.

We can use the observations and theory developed here to produce a microstructural explanation for the observed rheological changes as *ϕ* increases (figure 11). The most important microstructural change that occurs with increasing *ϕ* is a decrease in the typical minimum separation between nearest neighbour particles. Maximum packing fraction *ϕ*_{m} is reached when neighbouring particles are in permanent contact, forming a jammed network with a typical minimum separation of zero. Although figure 11 uses spherical particles to illustrate the changes in suspension rheology, the same microstructural processes also occur in suspensions of non-spherical particles. The chief difference is that non-spherical particles occupy a greater volume as they rotate in the flow than spheres at the same concentration, hence the typical minimum separation at a given *ϕ* is lower, and hence *ϕ*_{m} is also lower. A striking feature of the analysis of the data for consistency *K* and yield stress *τ*_{0} shown in figure 9 is the collapse of the curves to a single curve when *ϕ* is normalized by *ϕ*_{m}. The implication of this is that the ratio *ϕ*/*ϕ*_{m} represents a measure of the typical minimum separation that is independent of particle aspect ratio and so each part of figure 11 can be identified with a range of *ϕ*/*ϕ*_{m} independent of particle aspect ratio.

At very low particle volume fractions (figure 11*a*), the suspended particles are sufficiently well separated that interactions between them are negligible. In this regime, the particles increase the viscosity by causing the fluid to do extra work in flowing around them. Particles rotate in the flow according to the equations presented in Jeffery (1922). Rheology is Newtonian (*n*=1, *τ*_{0}=0, hence *K*=*η*) and viscosity (consistency) increases linearly with particle volume fraction (figure 8*b* inset). At slightly higher particle volume fractions (figure 11*b*), we find that rheology is still Newtonian but the viscosity now increases nonlinearly with *ϕ* (figure 8*b*). This is because particles are now close enough that they can interact hydrodynamically. Particles form doublets and higher multiplets that rotate together in the flow enhancing the distortion of the flow lines and increasing viscosity.

At intermediate particle volume fractions (figure 11*c*) shear-thinning is observed (*n*<1). The typical minimum separation between particles is sufficiently small that squeezing-flow in the small gaps between approaching particles causes local viscous heating, lubricating the passage of the particles leading to a decrease in viscosity at high strain rates (§4*b*). Shear-thinning occurs at lower values of *ϕ*/*ϕ*_{m} for particles of high aspect ratio than for spheres (figure 9 and §5*c*).

At still higher particle volume fractions (figure 11*d*), a yield stress is developed (*τ*_{0}>0) (figure 8*d*). This is because of the formation of chains and networks of touching particles. The magnitude of the yield stress is likely to be a function of particle size and the coefficient of friction between particles (§5*b*). The networks are capable of accommodating stress elastically until the yield stress is reached. At this point the networks break down and shear-thinning behaviour is observed as discussed above.

At the highest particle volume fractions (figure 11*e*), once the suspension has yielded, shear-thinning is less important and the stress–strain rate relationship is, once again, approximately linear; note the rise in *n* for large *ϕ*/*ϕ*_{m} in figure 9*f*. A plausible explanation for this is that the networks do not break down fully, but particle-free shear planes develop, which accommodate the bulk of the strain. Since the fluid in the shear planes is Newtonian, shear-thinning is less pronounced.

## Acknowledgements

We thank Antonio Costa for his helpful contributions. The comments of two anonymous reviewers helped to improve the paper. S.M. is supported by NERC grant NE/D003091/1, E.W.L. is supported by NERC Research Fellowship NE/D009758/2.

## Footnotes

- Received August 22, 2009.
- Accepted November 11, 2009.

- © 2009 The Royal Society