Rheol Acta (2010) 49:1059–1069

DOI 10.1007/s00397-010-0480-7

ORIGINAL CONTRIBUTION

A smoothed particle hydrodynamics-based fluid model

with a spatially dependent viscosity: application to flow

of a suspension with a non-Newtonian fluid matrix

Nicos S. Martys · William L. George ·

Byong-Wa Chun · Didier Lootens

Received: 30 October 2009 / Revised: 16 August 2010 / Accepted: 18 August 2010 / Published online: 9 September 2010

© US Government 2010

Abstract A smoothed particle hydrodynamics approach is utilized to model a non-Newtonian fluid with

a spatially varying viscosity. In the limit of constant

viscosity, this approach recovers an earlier model for

Newtonian fluids of Español and Revenga (Phys Rev E

67:026705, 2003). Results are compared with numerical

solutions of the general Navier–Strokes equation using

the “regularized” Bingham model of Papanastasiou

(J Rheol 31:385–404, 1987) that has a shear-ratedependent viscosity. As an application of this model,

the effect of having a non-Newtonian fluid matrix, with

a shear-rate-dependent viscosity in a moderately dense

suspension, is examined. Simulation results are then

compared with experiments on mono-size silica spheres

in a shear-thinning fluid and for sand in a calcium

carbonate paste. Excellent agreement is found between

simulation and experiment. These results indicate that

measurements of the shear viscosity of simple shearrate-dependent non-Newtonian fluids may be used in

simulation to predict the viscosity of concentrated suspensions having the same matrix fluid.

N. S. Martys (B) · W. L. George

National Institute of Standards and Technology,

100 Bureau Drive, Stop 8615,

Gaithersburg, MD 20899-8615, USA

e-mail: nicos.martys@nist.gov

B.-W. Chun

Grace Construction Products,

Cambridge, MA 02140, USA

D. Lootens

SIKA Technology A.G.,

Tuffenwies 16, 8048 Zürich, Switzerland

Keywords Smoothed particle hydrodynamics ·

Suspension · Shear-thinning fluids ·

Non-Newtonian fluids

Introduction

The flow- and time-dependent properties of nonNewtonian fluids are important in many biological,

environmental, and technological processes including

blood flow, food processing, and the placement of

concrete. In many cases, the non-Newtonian fluid of

interest is a suspension of rigid bodies immersed in

a matrix fluid. Such systems exhibit a wide variety of

phenomena such as shear thinning (Larson 1999), shear

thickening (Bender and Wagner 1996; Barnes 1989),

shear banding (Moller et al. 2006; Fielding and Olmsted

2004), dilatancy (Barnes 1989), and jamming (Lootens

et al. 2003). Most theoretical studies assume that the

matrix fluid of the suspension is Newtonian, where the

constitutive relation describes a linear stress evolution

as a function of shear rate with zero yield stress. Except

for the case of dilute (i.e., low-volume fraction) suspensions, the inherent complexity of suspensions makes

it difficult to predict the rheological properties from

purely theoretical means. Instead, there is an effort to

develop computationally based models to help improve

our understanding of complex fluids. The detailed modeling of suspension flow is challenging because the time

evolution of a rigid body’s motion entails keeping track

of a moving boundary at the fluid-solid interface. Still,

many advances have been made in understanding and

predicting the flow properties of suspensions where the

matrix fluid is Newtonian as an array of approaches
1060 Rheol Acta (2010) 49:1059–1069

have been developed including traditional computational fluid dynamics, Stokesian dynamics (Foss and

Brady 2000) and Lagrangian-based approaches like

dissipative particle dynamics (DPD) (Hoogerbrugge

and Koelman 1992; Boek et al. 1997; Martys 2005).

However, for many suspensions of interest, the matrix

fluid, for instance a cement paste (Flatt et al. 2004)

or polymeric-based material, is non-Newtonian. Often

such fluids can be described using Bingham, Herschel

Bulkley, or power-law models that have a shear-ratedependent stress or viscosity with a yield stress. The

computational modeling of such suspensions is greatly

complicated by the fact that local shear rates and hence

the viscosity can widely vary throughout the suspension. This feature is not accounted for in current computational models of suspensions.

In this paper, a numerical model is presented for

modeling non-Newtonian fluids with application to

flow of a suspension with a non-Newtonian fluid matrix

using a smoothed particle hydrodynamics (SPH)-based

approach. This model builds on and, in some respects,

is a synthesis of results found in Monaghan (2005) and

Español and Revenga (2003). For this paper, the fluid

modeled may be described as having a viscosity that is

allowed to have a spatial dependence. That is, the viscosity may, for example, depend on the local shear rate

or temperature which can, in turn, vary spatially. It is

shown that the model agrees well with direct numerical

solution of the general Navier–Strokes equation for an

idealized Bingham model of Papanastasiou, in the cases

of Couette and Poiseuille flow. As an application of

this model, the effect of having a non-Newtonian fluid

matrix, with a shear rate dependent viscosity in a moderately dense suspension, is then examined. Simulation

results are compared with experimental data of monosized silica particles in a shear-thinning fluid matrix. We

also consider the case of non-spherical particles where

simulation results are compared with experiments on

sand particles in a calcium carbonate paste. This paper

is organized as follows. Firstly, the numerical equations for the simulation are derived and the computational model is then described. The simulation model

is then validated for some simple flows like Couette

and Poiseuille flow using the “regularized” Bingham

model of Papanastasiou (1987) which has a shear-ratedependent viscosity. Using the same model fluid, a

sphere suspension is then simulated to understand the

consequences of having a non-Newtonian fluid matrix.

Next, comparison between experiment and simulations

are given for the cases of monosize spherical silica

particles in a shear-thinning fluid and finally for sand

in a calcium carbonate/water paste.

Numerical model

The model presented in this paper is based on the

SPH approach. SPH is a Lagrangian formulation of the

Navier–Stokes equations that was originally developed

by Gingold and Monaghan (1977) and Lucy (1977)

to study astrophysical problems. Recent advances include the inclusion of thermal fluctuations (Español

and Revenga 2003) in hydrodynamic flows and models describing viscoelastic fluids. While SPH has many

guises, the approach used in this paper largely relies on

the discretization of an integral representation of the

Navier–Stokes equations. The basic idea is to construct

an integral equation, composed of the fluid variables,

and a specified weight function having properties such

that when suitable dynamical variables are Taylor expanded in the integral, key parts of the Navier–Stokes

equations are recovered. Once a satisfactory form of

integrand is determined, the integral is represented as

a discrete summation over a weighted set of points.

Each point corresponds to a SPH “particle” which

carries information related to local fluid properties

(e.g., density, velocity...). For example, consider the

Taylor expansion, about r of the following integral:



dr





f(r



) − f(r)

r

− r



F



|r

− r|



≈



dr





r

− r



· ∇ f(r)(r

− r)F



|r

− r|



= ∇ f(r) + O(3). (1)

Here, f(r) could be a function that describes the fluid

and F(r) has the following property,  drF(r)rr = 1.

Additional relationships and properties of F(r) are described in the Appendix. This integral is an alternate

representation of the gradient of the function f(r) up to

O(3) in gradients. The next step is to approximate this

integral as a summation over neighboring points noting

that to each point, rq, there is associated a fluid property, fq. Consider, first, a “smoothed” representation

of f(r), given as

f(r) = 

q

fq S



|rq − r|

 (2)

with

S



rq − r



= W

|rq − r|





p W

|rp − r|

 . (3)
Rheol Acta (2010) 49:1059–1069 1061

Here, W(|rq − r|) is a weight function such that number

density of particles, dp, can be written as

dp = 

q

W|rq − rp|



, (4)

thus insuring that the representation of f(r) is properly

normalized. The fluid density, ρp, at point, rp, is then

given by ρp = mdp, where, m, is the mass. It can be

shown that the gradient of f at a point labeled p can

be written

∇ fp=∇ f(rp)= 1

dp



q



fq − fp

rq − rp

F



|rq − rp|



,

(5)

where ∇W(r) = −rF(r). On inspection, it is clear that

Eq. 5 is consistent with Eq. 1.

Note, in this paper i and j will be reserved to indicate components of vectors, whereas p and q will

correspond to actual SPH particles. Also, there is no

distinction between raised or lowered indices.

Using this approach, the time evolution of a fluid

may be represented as a set of particles, located at rp,

carrying local fluid properties, that undergo motion

in response to effective “interparticle” forces that are

defined by a descretized version of an integral representation of the Navier–Stokes equations. To achieve this,

consider the Lagrangian formulation of the continuity equation and the general Navier–Stokes equations

(Landau and Lifshitz 1987):

∂ρ

∂t = −ρ∇ · v, (6)

and

ρ

∂vi

∂t = −∂ P∂xi

+

∂

∂xk



μ

	 ∂vi

∂xk

+

∂vk

∂xi

− 2

3

δik∇ · v




+

∂

∂xi

(ζ∇ · v). (7)

Here, ρ, is the fluid density, P is pressure, v is velocity,

μ and ζ are the shear and bulk viscosities, respectively.

In these equations, the bulk and shear viscosities cannot

be taken outside the gradient operator because they can

be spatially dependent. In the limit that the viscosity

is a constant (Eq. 7) reduces to the usual Navier–

Stokes equations. Note, the Lagrangian formulation

is preferred because this approach can give us more

flexibility in handling moving boundaries as shown in

a previous work for the case of a DPD-based model.

Given the result shown in Eq. 5, it is easy to see that

a discrete representation of the continuity equation is

Monaghan (2005):

	∂ρ

∂t




p

= m

q

F



|rp − rq|

rp − rq



·



vp − vq



. (8)

To construct an integral representation of the general Navier–Stokes equation, a good starting point

(Zhu et al. 2010)1 is to consider the following integral

equation:2

Ai = 5



dr





μ(r



) + μ(r)


v j

(r



) − v j(r)




×



(r− r)i(r− r)j

(r− r)2



F(|r

− r|). (9)

Expanding μ(r) and v j(r) about r and taking advantage of the properties of F(r) one obtains

Ai =∇· μ∂iv



+∇· μ∇vi



+ ∂i



μ∇ · v



. (10)

Incorporating this approximation into Eq. 7 the following result is obtained

ρ

∂vi

∂t = −∂ P∂xi

+ Ai +

∂

∂xi

		ζ − 5

3

μ




∇ · v




= −∂ P

∂xi

+ Ai + Bi (11)

with Bi = ∂

∂xi



(ζ − 5

3μ)∇ · v



. Hence to O(3) in gradients, the Taylor expansion of Eq. 9, alone, i.e., Bi = 0,

recovers the velocity gradient portion of the Navier–

Stokes equations with the assumption that ζ = 5

3μ. The

SPH representation of Ai is given by

Ai

p =5



q

F



|rp − rq|



ρq

×



μp + μq

rp−rqi

rp−rq



·



vp−vq





rp−rq

2



. (12)

1This earlier paper provides an alternate SPH-based construction

of the general Navier–Stokes equations by utilizing the representation of gradients, as found in Eq. 1. As a result, obtaining

second order derivatives requires repeated application of Eq. 1.

The formulation of the general Navier–Stokes equations in the

current paper avoids the direct construction of many of the

second-order derivatives through the use of Eq. 9.

2Note, a similar integral form is used for the heat equation with a

spatially dependent thermal conductivity (Monaghan 2005) and

for the Navier–Stokes equations with constant viscosity (Español

and Revenga 2003).
1062 Rheol Acta (2010) 49:1059–1069

To incorporate the third term on the right of Eq. 11 so

that the ratio of shear to bulk viscosity is not fixed to 5

3 ,

one may use the derivative of a product property and

the following approximations. Taking a = ζ − 5

3μ then

Bi = ∂

∂xi

(a∇ · v) = (∂ia)∇ · v + a∂i∇ · v. (13)

Using standard SPH-based representations of derivatives and gradients (Monaghan 2005) in Eq. 11, one

obtains:



Ba

i



p = (∂ia)p

= m

ρp



q

F



|rp − rq|

rp − rq

i



ap − aq



, (14)



Bb

i



p = (∇ · v)p

= m

ρp



q

F



|rp − rq|

rp − rq

·



vp − vq



, (15)

and from Español and Revenga (2003), we deduce



Bc

i



p =

a∂i∇·v



p = map



q

F



|rp−rq|



ρq

×



5



rp−rq

i

rp−rq



·



vp−vq





rp−rq

2



−vp−vq

i



,

(16)

so that

Bi

p = 

Ba

i



p

∗



Bb

i



p + 

Bc

i



p . (17)

An alternative representation would be to determine

a∇ · v for each SPH particle and then apply the gradient

representation in Eq. 5. However, the Eq. 17 version

will be used for this work. Note that the term, ∇ · v, can

be directly obtained from evaluation of the continuity

equation so that no additional calculation is needed.

To account for the pressure gradient term, we follow

the previous work of Monaghan (2005):

(∇ P)p = −mρp



q



Pp

ρ2

p

+

Pq

ρ2

q



F



|rp − rq|

rp − rq



.

(18)

A pressure term, commonly used to model incompressible fluids, is given as P = c2(ρ − ρeq) where c is

related to the speed of sound and ρeq is an equilibrium

density.

When determining the local value of the viscosity it

will be first necessary to evaluate the local shear rate

tensor. A representation of the discretized shear rate is

given by:



γ ̇ij

p = 

q

F



|rp − rq|



ρq/m



rp − rq



i



vp − vq



j

. (19)

The local magnitude of the shear rate is then

γ ̇p =



ij 

γ ̇ij2

p

2 . (20)

Note that in the limit of constant viscosity, the momentum equations are the same as those derived by

Español and Revenga (2003), so this approach is a slight

generalization of theirs. In another work (Monaghan

2005), the Bi term is missing, which, again, is the same

as taking ζ = 5

3μ or a = 0. In this work, we will take the

bulk viscosity as zero, that is, we set a = −5

3μ. Including

the bulk viscosity, as in previous models (Monaghan

2005), by effectively taking the parameter a = 0, can

introduce oscillations in the flow field which is not

studied this paper.

Model for rigid body inclusions

The model for colloidal particles largely follows that

used for DPD-based simulations in Hoogerbrugge and

Koelman (1992) and Martys (2005). As in the case

of the DPD-based models, a colloid is defined as an

assembly of constrained SPH particles so that they form

a rigid body. The rigid body motion is then determined by summing the forces due the neighboring SPH

particles and other auxillary forces (e.g., lubrication,

colloidal and body forces). As described previously in

Boek et al. (1997) and Martys (2005), when modeling a

dense suspension of hard spheres, DPD, or in this case

SPH, particle interactions are not sufficiently strong

enough to prevent overlaps of the colloidal spheres. To

model realistic flows between neighboring spheres in

very close proximity would require a very fine resolution, or a number density of particles that is too high

to make simulations tractable over reasonable times.

To avoid this problem, lubrication forces (Kim and

Karrila 1991; Martys 2005) are explicitly included in the

simulation to account for hydrodynamic interactions

between neighboring spheres. There are several forms

or “modes” of lubrication force interactions between

hard spheres. The most well known and important is

called the squeeze mode, that accounts for forces that

develop as two spheres directly approach each other.
Rheol Acta (2010) 49:1059–1069 1063

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, is equal

to 3

2πμa2

r (VA − VB)/sAB, where ar is the sphere radius, VA and VB are the velocities of spheres “A”

and “B”, and sAB is the nearest surface to surface

distance between spheres labeled A and B. There are

additional contributions to the squeeze mode as well as

other modes like the twist mode, which, as it sounds,

describes the effect of one sphere rotating relative to a

neighboring sphere. These additional contributions all

scale logarithmically and are largely dominated by the

squeeze mode.

For the case of hard spheres, all the modes of the

lubrication force, up to the first order, including terms

that scale as 1/sAB, ln sAB and sAB ln sAB, will be used in

simulations. When the colloidal particle is not spherical,

the lubrication force must be modified to account for

local curvature. In this case the lubrication force is

limited to the squeeze mode for this paper.

An additional modification of the lubrication force is

made by adjusting the viscosity of the lubrication force

term to account for the local shear rate between rigid

bodies near the sphere surface. This modification is

admittedly heuristic, however, as a first approximation

of shear-thinning effects of two rigid bodies in near

contact, may be considered reasonable. Further study

of this effect is needed.

Although lubrication forces will prevent overlaps

of spheres, it is well known that in the course of a

simulation such rigid bodies are not restricted from extremely close approaches. When the spheres are nearly

in contact, very refined time steps are needed such that

the simulation becomes too time consuming. To avoid

this problem, yet still account for the lubrication force

effects, a short range repulsive force, FAB ∼ C

s8

AB

, directed radially between colloids A and B, was included

in the simulation. The functional form of repulsion is

akin to a DLVO approximation of a Lennard–Jones

hard sphere repulsion force (Israelachvili 1992).

All simulations will be carried out in the limit of

low Reynolds number. Re = ρLv/μ = ρa2

r γ/μ <<  ̇ 1

where, a, is the sphere radius. As there are no thermal

fluctuations present in this simulation, the results will

correspond to the infinite Peclet number limit.

Numerical integration approach

A modified velocity Verlet time integration algorithm

is used for both the “fluid” particles and the rigid

body inclusions (Groot and Warren 1997; Martys 2005;

Martys and Mountain 1999). The basic algorithm is

given in the following equations for an incremental time

step δt:

x(δt) = x(0) + v(0)δt +

δt

2

2

a(0), (21)

v(δt) = v(0) +

δt

2



a(0) + a(δt)




, (22)

where a(0) = f(x(0), v( ̃ 0))/m is the acceleration of the

SPH particle due to an effective force f described in

the next section and v ̃ is taken to be (Groot and Warren

1997)

v(δ  ̃ t) = v(0) +

δt

2

a(0). (23)

The translation of the rigid body inclusion was also

updated using a similar modified velocity-Verlet approach while a quaternion based algorithm accounted

for the rotation of the object. Details of this algorithm

are beyond the scope of this paper and are given in

Martys and Mountain (1999). It should be mentioned

that this approach has been used in previous studies

of suspension flow based on the DPD method. The

SPH-based computational model described in this paper is similar in structure to that of DPD so that it

was relatively easy to adapt the SPH model into the

DPD-based code.

Stress tensor

To evaluate the rheological properties of the composite

fluid (i.e., fluid matrix and rigid bodies combined),

it is convenient to construct a stress tensor based on

Kirkwood’s approach and similar to that used in molecular dynamics and DPD (Hoogerbrugge and Koelman

1992; Martys 2005; Allen and Tildesley 1987). We first

rewrite:

m

	∂vi

∂t




p

= m

ρp



−

	∂ P

∂xi




p

+ (Ai)p + (Bi)p



, (24)

as

m

	∂vi

∂t




p

= 

q=p

fi

pq, (25)

where m is a mass that will be taken to equal 1 in

the simulation. We can think of each neighboring SPH

particle, labeled q, contributing a force fpq on SPH

particle, p, as all terms on the right side of Eq. 19 result
1064 Rheol Acta (2010) 49:1059–1069

from a summation over the neighboring particles. The

stress tensor for this model suspension has several contributions. There are contributions from the propagation of momentum and inter-particle forces of the SPH

particles that are given by

σij = 1

Vtm



q

p ̃

q

i p ̃

q

j +

1

2Vt



p,q

fi

pq(r p − rq)j, (26)

where p ̃q is the momentum of particle q relative to the

macroscopic velocity field and Vt is the total volume of

the system. There are corrections to the stress tensor

due to the constraint forces on the SPH 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.

Determination of the correction for constraint forces is

described in Martys and Mountain (1999).

In addition, there are contributions to the stress tensor from the colloidal particle interactions. The stress

contribution from a central force between spherical

bodies is given by

σc

ij = 1

2V



A,B

Fi

AB(rA − rB)j, (27)

where A and B refer to two different colloids and FAB

is the force between colloids particles A and B described above. Stress calculations from lubrication

forces are also included. The forms of these equations

are described in Kim and Karrila (1991). The viscosity

of the entire system is then obtained from evaluation of

the total stress tensor μ = σxy/γ ̇.

For the simulations in the study, the number density

of SPH particles ranged from 3 to approximately 4.6

(where length is defined in units of h, see Appendix).

For a simple Couette flow the lower value was adequate

to reproduce the correct viscosity as a function of shear

rate. To improve resolution for the spheres simulation,

the higher value of density was used.

For the case of simulating Couette flow, a Lees–

Edwards Allen and Tildesley (1987) boundary condition was utilized. This boundary condition is equivalent

to application of a constant rate of strain at the boundaries. Other boundary conditions will be described as

utilized.

Regularized Bingham fluid

A simple model non-Newtonian fluid of Papanastasiou

(1987) was used for preliminary testing of this model.

For this model the shear rate dependent viscosity is

given by

μp = μpl + μ0

M



1 − exp −Mγ ̇p




γ ̇p

, (28)

where μpl is the viscosity in the limit of infinite shear

rate, usually called plastic viscosity, μ0 + μpl is the

viscosity in the limit of zero shear rate. The factor M,

can be thought of as a time scale where the suspension recovers its low shear rate behavior. Although

approximating a Bingham fluid, this model actually has

a zero yield stress. This model can be thought of as

a regularization of the Bingham model, which avoids

the singularity in viscosity in the limit of zero shear

rate. Unlike viscoelastic fluids, there is no inherent

memory effect in this fluid as the viscosity is known

immediately from the local shear rate. Such effects have

been considered elsewhere for model SPH fluids and

could be incorporated into this approach.

A dimensionless number called the Bingham number is given as Bn = μ0L

MμplV where L is the simulation

cell width and V is a characteristic velocity. For the case

of Couette flow V

L =  ̇γ .

Results

Couette flow

As a first test of this simulation approach, the case of

Couette flow is studied for the model of Papanastasiou with a Bn = 0.1. Figure 1 shows early and latter

stage flow fields for this model fluid. Note, due to the

application of the Lees–Edwards boundary condition,

there is a higher shear rate near the region of the

initial application of the shear, thus causing the fluid to

temporarily have a lower viscosity there. On the other

hand, because the fluid is more viscous at low shear

rates, there is much less flow in the middle portion

of the simulation cell making the flow profile appear

flat there. However as time progresses the entire fluid

evolves to the same shear rate and a linear velocity

profile finally develops. The evaluation of the stress

tensor for this latter stage flow produces a viscosity

consistent with that given by the constitutive relation.

Similar results were obtained for other values of Bn.

Poiseuille flow

A second test was to study Poiseuille flow, where

a body force-driven flow between parallel plates is
Rheol Acta (2010) 49:1059–1069 1065

-0,4 -0,2 0,0 0,2 0,4

-0,4

-0,2

0,0

0,2

0,4

V/Δ

Vmax

X/L

Bn=0,1

Fig. 1 Velocity profile, normalized to the maximum velocity

difference across the simulation cell, V/	Vmax, of width, L,

for Couette flow at early (f illed circles) and late stage (open

squares) flow. X is the distance from the zero velocity plane in the

simulation cell. Note for the early stage flow, near the application

of shear on the left and right hand sides of figure, the shear rate

is higher so that the viscosity is lower here. In the middle region,

at early stages, the velocity profile is more so horizontal as the

viscosity is higher here. At later stages of flow, the velocity profile

becomes linear

modeled. This requires the introduction of a no-slip

boundary condition at the fluid plate interface. To do

this, the following method described in Martys (2005)

is utilized. Here, the simulation cell is divided in half

and a body force is applied in one direction, in one half

of the simulation cell, and in the opposite direction in

the other half. Figure 2 shows the velocity profile in a

half cell for a set of Bn. Also shown is the numerical

solution of the general Navier–Stokes equation using

the same constitutive equation. Again agreement between simulation and numerical solution is reasonably

good. Note that there is a small variation on either side

of the velocity profile near the “walls”. This example

illustrates the sensitivity of the simulation to resolution

and is worth indicating. This artifact is a consequence

of not using a perfectly symmetric initial configuration

of points in the simulation cell. While slightly different

velocity profiles are obtained the agreement is still very

good with the independent numerical solution, thus

demonstrating the robustness of algorithm.

Suspension with an idealized Bingham fluid

Next, some general features of suspension flow for the

regularized Bingham model of Papanastasiou (1987)

will be discussed. Parameters for the Papanastasiou

model were chosen such that the viscosity limit for the

-0,5 -0,4 -0,3 -0,2 -0,1 0,0 0,1 0,2 0,3 0,4 0,5

0,0

0,2

0,4

0,6

0,8

1,0

Bn=0.3

Bn=0.1

Bn=0

V/Vmax

X/W

Bn=0.2

Fig. 2 Equilibrium Poiseuille flow fields for different Bn. The

open symbols represent results from SPH simulations. The solid

lines are numerical solutions of the generalized Navier–Stokes

equations. W is the width of the cell and X is the distance

from the center line of symmetry of the cell. The velocity, V, is

normalized to the maximum velocity of the case of Bn = 0. The

Bn = 0 case corresponds to the Newtonian fluid solution which

is parabolic. As Bn, increases the central portion of the velocity

profiles becomes flattened

low shear rate regime was about 100 times greater than

the high shear rate regime. For many cases this model

parameter choice would be a good approximation of

a Bingham fluid. Figure 3 shows the viscosity of the

fluid normalized to the infinite shear-rate value of the

viscosity. A set of monosize spheres, making a solid

volume fraction of 0.40, are added to the simulation

cell. Figure 3 also shows the relative viscosity of this

suspension. Again, the suspension viscosity is given

relative to the infinite shear rate value of the matrix

fluid. While the general trend is that the viscosity of the

suspension decreases with shear rate, there are three

regimes of interest for this study. Consider comparison

of the suspension viscosity with that of the matrix

fluid for the same shear rate. If the matrix fluid was

Newtonian, the suspension viscosity should be about a

factor of ten higher than that of the matrix fluid (Foss

and Brady 2000; Martys 2005). Note that at the very

low or very high shear rates studied, this limit is being

approached in the simulation. The reason is that, at the

highest shear rates, the relevant viscosity values sampled through out the simulation cell correspond to the

high shear rate limit of the matrix fluid. Other regions

in the suspension do not contribute as significantly to

the overall stress in the system. A similar argument

can be made at the lowest shear rate regime where

the viscosity will not vary as much through out the
1066 Rheol Acta (2010) 49:1059–1069

10-3 10-2 10-1 100101

100

101

102

103

μ/

μp

Mγ

.

Fig. 3 Comparison of matrix fluid viscosity (solid line) and relative viscosity of a suspension composed of monosize spheres with

volume fraction 0.40 (open squares). Values of viscosity, μ, for

both the matrix fluid and suspension are normalized to the plastic

viscosity μp of the matrix fluid. The matrix fluids properties

are based on the regularized Bingham model of Papanastasiou

(1987). Parameters for the model were chosen such that viscosity

increased by a factor of 100 as the shear rate decreases. The

open squares represent simulation data for the sphere suspension.

At the lowest and highest shear rates, the relative viscosity

approaches that of a similar sphere suspension in a Newtonian

fluid (approximately a factor of 10 higher than the matrix fluid).

At the intermediate shear rates shown, the relative viscosity of

the suspension is only about a factor of 2 to 3 higher than the

corresponding matrix fluid viscosity

suspension. An interesting shear rate regime is away

from where the viscosity plateaus. Here, the suspension

viscosity is only about a factor of two to three higher

than the matrix viscosity at each shear rate. In this case,

it is found that, as the suspension is sheared, the local

shear rate between spheres is much higher than that of

the bulk or global shear rate. As a result, the viscosity

between the spheres takes on the high shear rate limit

so that the overall viscosity is lower than that expected

for the relative viscosity.

Comparison with experimental measurements

of suspensions with spherical silica particles

Now, comparison is made between simulation and

actual experiments on idealized sphere suspensions,

Firstly, a shear- thinning fluid was prepared using a

solution of 5% by volume Methylhydroxypropylcellulose in water. This fluid has some similarity with the

regularized Bingham fluid with a viscosity plateau at

low shear rates. However, for the data obtained, there

was no indication of a plateau at the highest shear rates

tested. Experimental data were obtained with an Anton

10-3 10-2 10-1 100101102103104 10-1

100

101

102

103

γ (s-1)

49%

39%

30%

19%

0%

μ (Pa.s)

.

Fig. 4 Comparison of experimental measurements (solid lines)

of viscosity, μ (Pa-s) as a function of the shear rate γ ̇, (1/s) for

a shear-thinning fluid and suspensions composed of silica sphere

at the designated volume fractions φ with simulation data. The

dashed line represents a fit to the experimental measurement of

the matrix fluid. The open symbols represent simulation data.

Based on sample preparation, experimental results varied less

10% from that shown in the figure

Paar rheometer3 either with a cone/plate (50 μm gap)

or with a plate/plate (1 mm gap) geometry of 25 mm

diameter. Rheological measurements were then made

on systems composed of the same matrix fluid but

with silica particles added for a set of volume fractions

ranging from 0.19 to 0.49. The silica particles are type

P0060-Beads with diameter ranging between 100 to

140 μm. Figure 4 shows the viscosity as a function

of shear rate of the matrix fluid and the suspension

made by adding the silica spheres to the matrix fluid.

Tests were performed to check for the influence of

sedimentation and migration on the measured viscosity.

We expected that effects due to sedimentation would

be small due to the relatively high viscosity of the

matrix fluid. As a check, measurements were made at

a fixed shear rate for two minutes and did not show

a significant change in viscosity Repeating the measurements with the same sample gave similar results.

We also found there was no hysteresis of the viscosity

as a function of shear rate when increasing and then

decreasing the shear rate. This indicated there was no

3Certain commercial equipment, instruments, or materials are

identified in this paper in order to specify the experimental

procedure adequately. Such identification is not intended to

imply recommendation or endorsement by the National Institute

of Standards and Technology, nor is it intended to imply that

the materials or equipment identified are necessarily the best

available for the purpose.
Rheol Acta (2010) 49:1059–1069 1067

significant sedimentation or migration of the particles

over the time scales of the measurement.

The viscosity versus shear rate for the experimental

matrix fluid was then fit to a function that was used to

define a constitutive law for the fluid. This function was

used in the simulation code to define the local viscosity

for a given shear rate. Figure 4 shows that the correct

viscosity is self consistently recovered by simulation.

Next spheres were added to the simulation cell at volume fractions close to that of the experiment. As can be

seen in Fig. 4, the simulation results follow the experimental trends closely. Note, at the lowest shear rates,

the viscosity of the sphere suspensions, for different

volume fraction, are consistent with that of that of a suspension with a Newtonian fluid (Foss and Brady 2000)

as the viscosity being sampled is approaching its plateau

value. On the other hand, as the shear rates increases,

the viscosity of the higher volume fraction suspensions

had a tendency to decrease at lower values of shear rate

than the lower volume fraction systems. This is a result,

for the higher volume fraction suspension, of the typical

spacing between spheres being smaller—hence, larger

relative shear rates are present between spheres that

result in lower localized viscosities.

Comparison with experimental measurements

of a model mortar suspension

A second fluid system was considered that is representative of many suspensions (i.e., cement based materials). Here the matrix fluid was calcium carbonate/water

slurry with a sand (F-95 (US silica)) added to make a

mortar like suspension. Of course the calcium carbonate slurry is also a suspension but the constituent particle sizes are much smaller (less than 10 μm) than that of

sand so that it can be thought of as a continuum fluid.

The mean sand size was about 0.2 mm in diameter and

ranged about a factor of 5 about the mean. A Bohlin

CVO coaxial cylinder rheometer with a gap of about

2.5 mm was used to measure the suspension viscosity.

In contrast to the previous data, the matrix fluid did

not show a plateau in viscosity at lowest shear rates

tested, but, did begin to indicate a viscosity plateau

at the higher shear rates. As in the previous case the

simulation reproduced the correct constitutive relation

for the matrix fluid. A model for the sand particles was

used in the simulation that was based on reconstructed

tomographic images (Garboczi et al. 2001) of a similar

sand (Fig. 5). The tomographic images were analyzed

so that the local surface properties (surface normals,

principle curvature) of the particles could be incorporated into the simulation in order to better estimate

Fig. 5 Visualization of aggregates used in simulation of a suspension composed of sand embedded in a calcium carbonate paste.

The volume fraction φ = 0.50. The sand shapes are based on

X-ray microtomography images of a related sand (Garboczi et al.

2001)

lubrication forces. For the two volume fractions studied

(0.40 and 0.50) there was good agreement (Fig. 6) between simulation and experiment, especially given that

10-3 10-2 10-1 100 101

100

101

102

103

Sand/CaCO3 paste=40%

CaCO3 paste

γ (s-1)

μ/μp

.

Sand/CaCO3 paste=50%

-

Fig. 6 Comparison of simulation (open symbols) results to experimental measurements (solid lines) of relative viscosity, vs shear

rate γ ̇, (1/s) for a suspension composed of sand embedded in

a calcium carbonate paste. The matrix fluid (bottom solid line)

and an approximate fit (dashed line) to the data are also shown.

Data for volume fractions φ = 0.40 and 0.50 are shown. Based on

sample preparation, experimental results could vary about 10%

from that shown in the figure
1068 Rheol Acta (2010) 49:1059–1069

the model inclusions were not necessarily identical to

the sand particles. Note the higher values obtained for

the simulation for the 0.50 volume fraction could in part

be attributable to a sensitivity of particle shape variation in the simulation. Clearly more study is needed

here but the results are encouraging.

Conclusion

In conclusion, a set of numerical equations for modeling suspensions with a non-Newtonian fluid matrix

has been presented. In a series of tests of this model,

very good agreement was found between simulation

and independent numerical solutions, and, with a set

of experiments performed on two different suspensions. This work demonstrates that one may measure

the rheological properties of a fluid and then use this

measurement as input into a simulation code to determine the properties of a suspension having the same

matrix fluid. Extending this model to other fluid types

(e.g., viscoelastic, temperature dependent viscosity...)

is straight forward. Further research is needed to better understand how to properly modify lubrication

forces to account for the presence of a non-Newtonian

fluid matrix. In addition, correct modeling of thermal

fluctuation effects to account for low Peclet number

flow would also be desirable as this is not currently

understood for the case of a non-Newtonian fluid.

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 also used resources

of the Argonne Leadership Computing Facility at Argonne

National Laboratory, which is supported by the Office of Science

of the US Department of Energy under contract DE-AC02-

06CH11357.

Appendix

Proper construction of the weight function W(r) is important for the physically correct transmission of matter

or forces between neighboring SPH particles. Here

some key properties of W(r) and its derivative are given

for convenience. The weight function and alternate

formulations are discussed more fully in Español and

Revenga (2003), Monaghan (2005).



drW(r) = 1. (29)

∇W(r) = −rF(r) (30)

In this work, the SPH Lucy function is utilized for a

weight function.

W(r) = 105

16πh3



1 + 3

r

h

 1 − r

h

3

(31)

and

F(r) =

	 315

4πh5




1 − r

h

2

. (32)

For for the remainder of this paper, we set h = 1. Some

of the more important properties of F(r) are given

below.



drF(r)rr = 1. (33)



drF(r)

xxxx

r2 = 35

. (34)



drF(r)

xxyy

r2 = 15

. (35)

References

Allen MP, Tildesley DJ (1987) Computer simulation of liquids.

Claredon, Oxford

Barnes HA (1989) Shear-thickening (“dilatancy”) in suspensions

of non-aggregating solid particles dispersed in Newtonian

liquids. J Rheol 33(2):329–366

Bender J, Wagner NJ (1996) Reversible shear thickening in

monodisperse and bidisperse colloidal dispersions. J Rheol

40(5):899–916

Boek E, Coveney PV, Lekkerkerker HN, van der Schoot P (1997)

Simulating the rheology of dense colloidal suspensions

using dissipative particle dynamics simulations. Phys Rev E

55:3124–3133

Español P, Revenga M (2003) Smoothed dissipative particle dynamics. Phys Rev E 67:026705

Fielding SM, Olmsted PD (2004) Spatiotemporal oscillations and

rheochaos in a simple model of shear banding. Phys Rev Lett

92:084502

Flatt RJ, Martys N, Bergstrom L (2004) The rheology of cementitious materials. MRS Bull 29(5):314–318

Foss DR, Brady JF (2000) Structure, diffusion, and rheology of

Brownian suspensions by Stokesian dynamics simulations.

J Fluid Mech 407:167–200
Rheol Acta (2010) 49:1059–1069 1069

Garboczi EJ, Martys NS, Saleh HH, Livingston RH (2001)

Acquiring, analyzing, and using complete three dimensional aggregate shape information. In: Aggregates, concretes, bases and fines, ninth annual symposium proceedings

(CD-ROM). International center for aggregate research,

April 22–25, 2001. Austin Texas, pp 13

Gingold RA, Monaghan JJ (1977) Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon

Not R Astron Soc 181:375–389

Groot RD, Warren PB (1997) Dissipative particle dynamics:

bridging the gap between atomistic and mesoscopic simulation. J Chem Phys 107:4423–4435

Hoogerbrugge PJ, Koelman JMVA (1992) Simulating microscopic hydrodynamic phenomena with dissipative particle

dynamics. Europhys Lett 19:155:96–160

Israelachvili JN (1992) Intermolecular and surface forces.

Academic Press

Kim S, Karrila SJ (1991) Microhydrodynamics. Butterwort–

Heinemann

Landau LD, Lifshitz EM (1987) Fluid mechanic. Pergamon,

Oxford

Larson RG (1999) The structure and rheology of complex fluids.

Oxford University Press, New York

Lootens D, Van Damme H, Hébraud P (2003) Giant stress

fluctuations at the jamming transition. Phys Rev Lett 90:

178301

Lucy LB (1977) A numerical approach to the testing of the fission

hypotheses. Astron J 82:1013–1024

Martys NS (2005) Study of a dissipative particle dynamics

based approach for modeling suspensions. J Rheol 49:401–

424

Martys NS, Mountain RD (1999) Velocity Verlet algorithm for

dissipative-particle-dynamics-based models of suspension.

Phys Rev E 59:3733–3736

Moller PCF, Mewis J, Bonn D (2006) Yield Stress and thixotropy:

on the difficulty of measuring yield stresses in practice. Soft

Matter 2006 2:274–283

Monaghan JJ (2005) Smoothed particle hydrodynamics. Rep

Prog Phys 68:1703–1759

Papanastasiou TC (1987) Flows of materials with yield. J Rheol

31:385–404

Zhu H, Martys NS, Ferraris C, DeKee D (2010) A numerical

study of the flow of Bingham-like fluids in two-dimensional

vane and cylinder rheometers using a smoothed particle hydrodynamics (SPH) based approach. J Non-Newton Fluid

Mech 165:362–375