Paper No 7t" International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Hindered settling in thixotropic liquids
Jos Derksen, Orest Shardt
Chemical & Materials Engineering, University of Alberta, Edmonton, AB, Canada T6G 2V4
jos@ualberta.ca, shardt@ualberta.ca
Keywords: nonNewtonian liquids, direct numerical simulation, latticeBoltzmann, liquidsolid suspensions
Abstract
The flow of suspensions consisting of monodisperse, solid spheres and nonNewtonian liquids has been studied via direct
numerical simulations. The liquids are purely viscous (i.e. nonelastic) with shear thinning and/or thixotropic (timedependent)
behavior. The interstitial liquid flow is solved by means of the latticeBoltzmann method. Thixotropy enters via a
networkintegrity parameter that relates to the local, apparent viscosity and for which a transport equation has been solved. In
the main portion of this paper the sphere assemblies are static. These simulations show that in terms of the drag force the
shearthinning character of the liquid manifests itself more pronounced at higher solids volume fractions. Thixotropy tends to
increase the drag force due to the decoupling of locations of high deformation rates and low viscosity. In addition, assemblies
of moving (settling) spheres have been studied and we show the profound role of the lubrication force (and the modelling
choices it involves) on the onset of instabilities.
Introduction
The waste streams in oil sands processing (called tailings,
see e.g. Masliyah, et al 2004) consist of clay suspensions
(fine, surface active particles in water), residual bitumen and
sand. Separation of these components is an important
environmental issue as it relates to land reclamation at the
end of the oil sands mining cycle. The separation of coarser
particles (flocks, sand) from the clay suspension is a
nontrivial and potentially troublesome process since clay
suspensions exhibit nonNewtonian (yield stress,
shearthinning, thixotropic) behavior. This behavior is due
to a structural network being formed as a result of
longrange interactions of clay particles, and the breakdown
of the network due to liquid deformation. Thixotropy
(timedependent rheology) enters the physical picture
because buildup and breakdown of the network is not an
instantaneous process; it is subject to ionictransport
limitations.
In this paper we explore a computational approach to the
above problem: We perform direct numerical simulations
(DNS's) with full resolution of the solidliquid interfaces of
dense suspensions consisting of solid spheres (representing
the coarse particles) in thixotropic liquids. In the first set of
simulations thixotropic liquids are forced through a static,
random assembly of monosized spheres. Goal of these
simulations is to measure the drag force as a function of
solids volume fraction and liquid properties.
In a subsequent set of simulations the monosized spherical
particles are allowed to move and we study their hindered
settling through thixotropic liquids in a narrow, fully
periodic domain. In such systems with Newtonian liquids,
wave instabilities are know to occur with voidage waves
traveling opposite to the direction of gravity (Duru et al.
2002). We expect an impact of thixotropy on the onset and
strength of the instabilities: the voidages may favor the
buildup of the structural network in the liquid, thus making
it more difficult for solids to penetrate the void.
This paper is organized in the following manner: we first
discuss our (relatively simple) thixotropy model. Its
parameters, together with the flow conditions allow us to
define a nondimensional parameter space. Then we briefly
go into the numerical procedure we are using that is based
on a latticeBoltzmann procedure for solving the flow
equations, twoway coupled to a finitevolume scalar
transport solver that keeps track of the integrity of the
network in the liquid, and an immersed boundary method to
deal with the solid particles. We then discuss drag force
results in static assemblies of spheres and the hindered
settling simulations. Finally we reiterate our main findings
in the Conclusions section.
Nomenclature
d, d
d
Db
FD F ,F",F"
f
k,,k,
kl, k2
L
Re,Re
S
U,
us, s
U
Greek letters
c++l
sphere diameter & "given" diameter
deformation tensor
Deborah number
Drag force, normalized drag forces
body force
thixotropy parameters (Eq. 1)
domain size
Reynolds number & Re based on
infinite shear viscosity
shear rate number dimensionlesss)
fluid velocity vector component
superficial fluid velocity
singlesphere settling velocity
viscosity ratio
Paper No
7, Y deformation rate & characteristic shear
rate
6 interparticle spacing
A, A, network integrity & in steady state
uao ,/u ,Au viscosity (apparent, zeroshear,
infiniteshear, steadystate).
v kinematic viscosity
P, p fluid & solids density
0, 0 solids volume fraction & at random
close packing
Thixotropy model
The thixotropy model we have adopted is based on early
work due to Storey & Merrill (1958), and Moore (1959),
more recently applied by Ferroir et al. (2004). In this purely
viscous model we keep track of a scalar A that varies
between 0 and 1 and indicates the integrity of the network
(A=0: no network; A=1: fully developed network). Its
transport equation reads
S+u = k + kj (1 A) (1)
at adx,
(summation over repeated indices implied) with u, the ith
component of the fluid velocity vector, and y = idd
1 fu" u,
a generalized deformation rate; d, = 2 au+'I is
2 ax, ax/)
the rate of strain tensor. The first term on the right hand
side of Eq. 1 indicates breakdown of the network due to
liquid deformation; the second term is responsible for
1
buildup of the network with a time constant 
k2
associated to it. In the model, the apparent viscosity /a is
linked to the network integrity according to
au, = A (1+ aA) (2)
with /u and a model constants.
In a homogeneous shear field with shear rate y, the
steadystate solution to Eq. 1 reads
k
k, = (3)
k y+kk2
The associated steady state viscosity is (combine Eqs. 2
and 3)
( k
lull =,. 1 +a k2 (4)
k kly+k2
The parameter tu_ can thus be interpreted as the infinite
shear viscosity. The zeroshear viscosity is
/0u = (1+ a) A typical representation of the
steadystate rheology (Eq. 4) for a > 0 is given in Figure
1. As can be seen, it represents a shear thinning liquid
making the transition from the zeroshear viscosity p0 to
infiniteshear shear viscosity u_ around a characteristic
71" International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
k2
shear rate a With the latter definition, Eq. 4 can
also be written as
/i =/ 1+ c+2/c)
7,r
O
0 Y
Figure 1. Steadystate rheology with infiniteshear viscosity
/. +/10
/_ zeroshear viscosity /,. At y= yc s =
2
Simulation procedure for static sphere assemblies
The flow systems we will be considering consist of random
solid sphere assemblies with liquid in between. All spheres
in the assembly have the same diameter d. The liquid
flows through the space not occupied by the spheres with a
superficial velocity us as the result of a body force f
acting on the liquid. The body force f is tuned such that in
steady state a preset value of us is achieved. The
boundaries of the threedimensional flow domain are fully
periodic. A force balance then implies that that the average
drag force FD on a sphere in the assembly is (Ladd, 1994)
FD = f ff d 3 (6)
6 )
with 0 the solids volume fraction of the assembly. This
procedure usually in the context of a latticeBoltzmann
scheme for solving the flow has been applied by a number
of authors to study drag forces due to Newtonian fluid flow
as a function of solids volume fraction and/or Reynolds
number (Kandhai et al., 2003; Van der Hoef et al., 2005;
Ladd, 1994; Hill et al. 2001). Note that Eq. 6 expresses the
convention for the drag force as e.g. expressed by Van der
Hoef et al. (2005) among others. The total average force by
the fluid on a sphere FfS then is the sum of FD and a
contribution from the body force (e.g. in its role of an
average pressure gradient): FS =FD +f 'd3. It can be
6
F
verified (Van der Hoef et al., 2005) that F, D
10
In Newtonian cases (with dynamic viscosity u ) the drag
force (Eq. 6) is commonly made dimensionless according to
F
F D Simulations are usually carried out with f
3ztpdu,
and thus FD acting in a specific direction with (as a result)
u, having the same direction. The scalar values FD and u
are the vector components in that direction. Computational
results for monosized spheres and low Reynolds numbers
Paper No
pu,d
(with Re = ) have been summarized in correlations in
/1
terms of FD, e.g. by Van der Hoef et al.(2005):
S 10 +(1 )2(1+1.5,) (7)
(10)2
The flow of thixotropic liquids (as defined in the previous
section) through monosized sphere assemblies can be
pinned down with a set of five dimensionless numbers. In
/1," .1
this paper these have been chosen as Re_ = ,
lu.
Db = S = a, and 0 (Db is the ratio of the
dk2 Us
1 d
time scale of the liquid and of the flow termed
k2 us
Deborah number; having a Deborah number does not imply
viscoelasticity). This fivedimensional parameter space we
only partly explore by only considering a single, low value
of Re =0.06 We also fix the viscosity ratio to
a+1 =16.
The particle configurations are created by randomly placing
a number of spheres in a cubic domain. Since we are
interested in dense suspensions (up to = 0.530) random
placement usually needs to be followed by compacting the
particle assembly. After compaction to the desired volume
fraction, the configuration is randomized again by giving the
spheres random velocities and letting the system evolve for
some time as a granular gas with the particles undergoing
fully elastic, frictionless collisions. Then the configuration is
frozen to be used in the fluid flow simulation procedure.
As in many of the earlier works on the subject of drag on
sphere assemblies, we used the latticeBoltzmann method
(LBM) (Chen & Doolen, 1998; Succi, 2001) to solve for the
flow of interstitial liquid. It has a uniform, cubic grid on
which fictitious fluid particles move in a specific set of
directions, and collide to mimic the behavior of an
incompressible, viscous fluid. The specific LB scheme
employed here is due to Somers (1993). The scheme can
accurately account for liquids with nonuniform viscosity as
recently demonstrated by Derksen & Prashant (2009). The
noslip condition at the spheres' surfaces was dealt with by
means of an immersed boundary (or forcing) method. We
have validated and subsequently used this method
extensively to study the interaction of solid particles and
Newtonian fluids (e.g. Ten Cate et al., 2002).
It should be noted, however, that having a spherical particle
on a cubic grid requires a calibration step, as first realized
by Ladd (1994). He introduced the concept of a
hydrodynamic diameter. The calibration involves placing a
sphere with a given diameter dg in a fully periodic domain
in creeping flow and (computationally) measuring its drag
force. The hydrodynamic diameter d of that sphere is the
diameter for which the measured drag force corresponds to
the expression for the drag force on a simple cubic array of
spheres due to Sangani & Acrivos (1982). Usually d is
slightly bigger than dg with ddg typically equal to one
lattice spacing or less. The issue of the hydrodynamic
diameter is particularly important in this study since the
71" International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
deviation of d from dg is a function of the viscosity of the
liquid (Ladd, 1994; Ten Cate et al., 2002) and since in this
study the liquid's viscosity varies over the sphere's surface
with an a priori unknown distribution. We have dealt with
this issue by assessing the resulting error in the drag force
and reducing it to an acceptable level. The relative error in
the drag force can be controlled by means of the spatial
resolution of the simulations. Since d dg typically equals
one grid spacing, increasing d relative to the grid spacing
dd
reduces the relative deviation and thus the relative
d
deviation in the drag force. In the next section the
quantitative details are discussed.
Having a thixotropic liquid requires solving a transport
equation in the network integrity parameter, Eq. 1. This we
do by means of a finite volume scheme on the same
(uniform and cubic) grid as the latticeBoltzmann
discretization. An advantage of employing a finite volume
formulation (over e.g. also using the LBM for scalar
transport modeling) is the availability of methods for
suppressing numerical diffusion. As in previous works,
TVD discretization with the Superbee flux limiter for the
convective fluxes (Sweby, 1984) was employed.
We step in time according to an Euler explicit scheme. The
deformation y= 2dd and the liquid velocity u
come from the LBM; the distribution of A is fed back to
the flow solver part by running the LBM with the apparent
viscosity field /u stemming from Eq. 2. This makes the
scalar transport and the LBM part of the simulation
procedure twoway coupled.
Verification
As explained above, the issue of the hydrodynamic diameter
versus the given diameter, and the dependence on the
viscosity of the former require careful consideration.
Calibrations according to Ladd (1994) were performed with
a Newtonian fluid at kinematic viscosity v = 0.04 (in lattice
units) such that a set of spheres was created with
hydrodynamic diameters of d=12, 16, 24, and 32 lattice
spacings.
14
FD
12
10
0 20 40
d (lu)
Figure 2. Dimensionless drag force FD for a Newtonian
simulation at Re=0.06 and 0 = 0.373 as a function of the
spatial resolution of the simulation expressed in the number
of latticespacings spanning the diameter d. Squares:
v =0.04 (the viscosity at which d was calibrated);
triangles: v=0.01;plusses: v=0.16.
Paper No
Spheres were placed in a random assembly in a cubic
domain with vertex length L = 6.25d at a solids volume
fraction 0 = 0.373 First a Newtonian liquid with v = 0.04
(i.e. the same viscosity with which the calibrations were
done) was forced through this assembly at the four
resolutions considered, all four at Re:
u d
d= 0.06.
v
The results are given in Figure 2 (square symbols). The
relative placement of the spheres is the same in the four
simulations; the only difference is the spatial resolution. In
terms of F the results at the three highest resolutions
differ less than 5% and can be considered grid independent.
Subsequently we performed simulations with a viscosity
being a factor of 4 higher, and a factor of 4 lower than the
default viscosity of 0.04 (still Newtonian fluid, still
Re=0.06). Now we clearly see the effect of the
hydrodynamic diameter being dependent on the viscosity:
At the lower viscosity F gets overpredicted (triangles); at
the higher viscosity underpredicted (plusses). As anticipated,
however, the deviations are a pronounced function of the
resolution. At the highest resolution (d=36) the deviations
are typically 3%, at the one but highest resolution (d=24)
+5%. These deviations are worstcase deviations. In our
nonNewtonian simulations the maximum ratio of highest
over lowest viscosity is a+1=16 with the center
11 +10
kinematic viscosity +/ equal to 0.085. Given the
2p
large parameter space and (consequently) large number of
simulations to be done it was decided to perform the
remainder of the simulations with a resolution of d=24, and
to keep the +5% error in mind when interpreting the
results.
30
20
10
0
0.3 0.4 0.5
Figure 3. F* as a function of solids volume fraction 0
for Newtonian flow at Re=0.06. Solid curve: Eq. 7
(correlation due to Van der Hoef et al., 2005); squares:
present simulations.
So far the performance of our simulations procedure has
been verified by internal comparisons (comparing
simulations done with the same computer code, however,
with different numerical and physical settings). Since the
immersed boundary method used in this paper (Ten Cate et
al., 2002; Derksen & Van den Akker, 1999) differs from
more common ways to impose noslip conditions at curved
surfaces in the latticeBoltzmann method (the defacto
standard being set by Ladd in 1994), drag force results have
also been compared with results from literature; more
specifically the results due to Van der Hoef et al. (2005) on
drag in random, monosized sphere assemblies as a function
71" International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
of the solids volume fraction 0 at creeping flow
conditions and with Newtonian liquids. They summarized
their results with the correlation given in Eq. 7. In Figure 3
that correlation is compared to our simulations at the default
L
settings (Re=0.06, resolution d=24, system size = 6.25,
d
v = 0.04 and the hydrodynamic diameter calibrated for this
viscosity) showing good agreement. The data points are the
average of three sphere configurations. Apparently the
resolution is sufficiently high, and the Reynolds number
sufficiently low to accurately reproduce the empirical
correlation (Eq. 7).
Drag force results
Now we turn to drag force results for nonNewtonian
liquids. As discussed above, of the five dimensionless
numbers governing the flow system, Re and the
viscosity ratio a+1 have fixed values, 0.06 and 16
respectively. Initially only shearthinning, timeindependent
liquids will be considered. Timeindependent liquids have
k2 o so that Db=0. The two degrees of freedom left are
the solids volume fraction 0 and S = All data points
u,
presented in this section of the paper are average values of
three statistically independent spherical particle
configurations.
o o o
Figure 4. xzcross sections through the flow domains, with
the horizontal (x) direction the direction of the mean flow.
Shearthinning, timeindependent liquids. Colors indicate
the apparent viscosity. From top to bottom =0.373, 0.459,
and 0.530. From left to right S=I, 4, and 16 respectively.
In Figure 4 we show contour plots of the distribution of the
apparent viscosity in a cross section through the suspension
after steady state has been reached. The cross sections span
the xzplane of the cubic periodic domain with the
xdirection the streamwise direction. The white circular
disks are cross sections thorough the spherical particles. The
higher S, the higher the apparent viscosity in the suspension
gets. Th i is not surprising. At higher S (and thus higher
y} ) the transition from zeroshear viscosity to fifiniteshear
Paper No
viscosity takes place at higher deformation rates (see Figure
1). Also the range of viscosities encountered in the
suspension is a function of S: if the characteristic shear rate
y~ of the liquid is of the same order as the shear rates
encountered in the interstitial liquid, a relatively broad range
of viscosities is anticipated.
The above observations are presented in a more quantitative
sense in Figure 5. The figure shows (for three solids volume
F
fractions) the doubly normalized drag force FD 
FDS
In these nonNewtonian cases FD is based on the
infiniteshear viscosity: FD D and F o (by
3f7p du s sO
definition) is the Newtonian drag force at the specific solids
volume fraction normalized with 3Trp.dus.
rms(p)
It :
iu iu i lu u d
Figure 5. Doubly normalized drag force FD* (top);
average apparent viscosity (au) (middle); and
rootmeansquare values of the apparent viscosity rms (/a)
as a function of S for three solids volume fractions: from
left to right: 0=0.373, 0.459, and 0.530. Shearthinning,
timeindependent liquids.
Figure 5 also shows the average apparent viscosity in the
suspension, and (as a measure of the spread in viscosities in
the liquid domain) the rootmeansquare (rms) values of the
deviations from the mean apparent viscosity:
rms() a) uu dV with Vf the volume
PVf
of (interstitial) fluid. Interestingly rms(/u,) goes through a
maximum with the location of the maximum dependent on
the solids volume fraction: the higher 0 the further the
maximum shifts to higher S. At higher 0 the space
between the spheres gets narrower and (since the superficial
velocity has a fixed value) the deformation rates in the
liquid increase. As a result, the distribution of viscosities
71" International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
gets wider for higher y, i.e. higher S.
The notion of the interaction between yc and deformation
rates in the suspensions proofs helpful in scaling the doubly
normalized drag force F' Critical regions in the
suspension are the waistshape gaps between neighboring
spheres. As a measure for the typical size of these waists we
take S d r j 1 with 0p the solids volume
fraction at random close packing. For the latter we took
rp = 0.62 (Wang & Ge, 2005). In Figure 6 we plot FD*
as a function of for all cases with timeindependent
us
rheology (Db=0) considered. The drag force behaves fairly
consistently over the wide range of solids volume fractions
considered (O from 0.330 to 0.530).
We now investigate the effect of thixotropy (timedependent
rheology) on the liquid flow and resulting drag force in the
suspension. The (now three) relevant dimensionless
u yd
numbers are Db = s S= and 0. As before we
dk2 Us
fixed Re = 0.06, and a+1= 16. The interpretation of
the thixotropic liquid results will be mainly by comparing
them with corresponding results with Db=0. In doing this,
we compare (per value of 0) simulations with the same
configuration of spheres only.
FD 0 330
S0.373 +
8 + 0420
x 0459
v 0 500
A 0 530 +7
A+
4
A+
0
102 10 10 101 102
us
Figure 6. Doubly averaged drag force as a function of
yc6/us with S defined in the text. The different symbols
relate to different solids volume fraction as indicated.
Shearthinning, timeindependent liquids.
Examples of steadystate results are given in Figure 7, that
shows a cross section through one of the suspensions in
terms of the apparent viscosity (note that the color scale in
Figure 7 is different from Figure 4). The most visible effect
of thixotropy is a smearingout of the viscosity fluctuations.
This effect sets in beyond Db=0.2 (the viscosity fields at
Db=0 and Db=0.2 are almost the same). The smearing out is
due to the time it takes to build up or break down the
network. In an infinitely fast (Db=0) liquid, locations where
the network is formed or broken down coincide with places
of respectively low (e.g. bigger voids in the suspension) and
high (shear layers at solid surfaces) deformation rates. If the
liquid needs time to respond to deformation conditions
Paper No
(Db>0) the breakdown and buildup processes are less
localized with a smoother apparent viscosity field as the
result.
,
0=0.420, S=4. Thixotropic simulations with (from left to
right and top to bottom) Db=0, 0.2, 1.0, 5.0, 25.0, 125.0.
Colors indicate the apparent viscosity.
In Figure 8 some of the results of thixotropic simulations are
displayed, with Db as the independent variable. The triply
normalized drag force is defined Th
smoothing effect on the viscosity field shows as nrms(, )
being a monotonically decreasing function of Db. The
average apparent viscosity in general slightly reduces as a
result of thixotropy (by some 20% max). The net effect is an
increase in the drag; in all cases considered: F** >1. The
effect of thixotropy on the drag force is not very big; the
maximum increase is approximately 20% (which still is
significantly more than our 5% uncertainty measure related
to calibration issues), occurring for high Db in situations
where the corresponding Db=0 system had a large
rms(,u ).
The increased drag force, also in cases for which the
average viscosity decreases demonstrates the relevance of
the viscosity distribution in the suspension. In the critical
regions between closely spaced sphere surfaces shear
thinning liquids reduce drag as a result of the locally high
deformation rates. Thixotropy largely decouples locations of
high shear with those of low viscosity thereby increasing the
stresses in the critical regions, and thus increasing drag.
Sedimenting spheres modelling approach
We now turn to simulations where the spherical, solid
particles are allowed to move in response to hydrodynamic
forces and hardsphere collisions. We again consider fully
periodic, threedimensional domains. The density of the
solid spheres is higher than the density of the liquid (in our
simulations A =4.1), and the spheres sediment through the
liquid. In order to globally balance forces, a body force acts
on the liquid in the direction opposite to gravity (in
fluidization this body force would be an overall pressure
gradient). The domain size is 20d in streamwise (is vertical)
rV? ko trv ? siuain hr heshrcl o
particles ~ (}Cllwd omoei' t*spos ohdoy
forces~, 7)d hadshr olsos eaancnie ul
periodi.K tredmn i tldmis hedniyo h
71" International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
direction, and 4d in the two lateral directions. The resolution
is such that d spans 16 lattice spacings. Noslip conditions
at the solidliquid interfaces are (as in the
staticspheressimulations) imposed by an immersed
boundary method. The forces required for maintaining
noslip are used to determine the hydrodynamic force and
torque on each sphere so that their equations of motion
(linear and rotational) can be numerically solved.
For suspensions involving Newtonian liquids, the
simulations procedure has been described in detail by
Derksen & Sundaresan (2007). That procedure has been
extended in order to be able to deal with thixotropic liquids,
according to the thixotropy model as described above: we
added a scalar transport solver to keep track of the network
integrity parameter A, and coupled A back to the flow
equations with an apparent viscosity according to Eq. 2.
:^L.
10 3 V
09
08
S02
1i 1
A 4+
+ 16
A 64v
V
a V A
V v8 1
1 1 " 4
10 3
09
101 100 101 102 103101 100 101 102 103
Db
k2d
Figure 8. Sample simulation results with thixotropic liquids.
From bottom to top: triply normalized drag force, rms of
viscosity, and average viscosity as a function of Db. Left:
0=0.530; right: 0=0.330. The different symbols denote
different values of S as indicated.
As discussed in Derksen & Sundaresan (2007), in the
movingparticle simulations the hydrodynamic forces
stemming from the LBM are explicitly supplemented with
lubrication forces. If two spheres move relative to one
another and their surfaces are in close proximity (typically
closer than one lattice spacing) the LBM does not resolve
the hydrodynamic interaction between the two spheres
anymore. We then add radial lubrication forces to the forces
acting on the spheres according to procedure due to Nguyen
& Ladd (2002). It was subsequently shown by Derksen and
Sundaresan that in their specific cases (that involved
Newtonian interstitial liquid) the momentum transfer due to
radial lubrication is small compared to collisional and
streaming stresses.
In the situation with spheres moving through thixotropic
liquids, the issue of lubrication forces needs to be
reconsideredd carefully. The analytical expressions for
lubrication forces (Kim & Karilla, 1991) are based on
Paper No
Newtonian fluids with the force proportional to the dynamic
viscosity. Here we have nonNewtonian fluids and an
apparent viscosity that changes in space and time. As a
preliminary approach we decided to check the sensitivity
with respect to lubrication model choices, i.e. we compared
simulations in which the viscosity as used in the lubrication
force expressions was equal to (1) the zeroshear viscosity
/t0 (2) the infiniteshear viscosity (3) the local
viscosity in a volume with size 0.375d centered at the point
of (near) contact. The results were compared in the way and
extent to which they developed planar wave instabilities.
To characterize the (single) physical situation we are
considering, the velocity scale is taken as the settling
velocity u_ of a single sphere in a Newtonian liquid
pu d
having viscosity /u We then have Re= =290;
Db=" =5.6, and S= =0.18. In addition the
dk2 u
viscosity ratio has been set to +1 = 41, the solids volume
fraction is =0.50, and the density ratio (see above) is 4.1.
Sedimenting spheres results
.34
0..1
0.1
0
Figure 9. Spheres settling in a thixotropic liquid.
Comparison between lubrication model choices. Left:
lubrication based on local apparent viscosity; middle: based
on pu_ ; right: based on /0 =u (1+ a) The bottom
panels show the local level of network integrity A.
Typical impressions of the flow configurations are given in
Figure 9. The snapshots in Figure 9 are taken 1.5105 time
steps (which corresponds to 170 times d/lu ) after start up
71" International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
from rest with the spheres homogeneously distributed in
space and A=0 everywhere in the liquid.
As for Newtonian liquids, the suspensions develops a planar
wave instability, somewhat visible in Figure 9. The onset
and strength of this instability, however, strongly depend on
the viscosity choice in the lubrication force model, see
Figure 10. There we show spacetime plots of the local
solids volume fraction with the inclined contours
representing waves. Where the simulation with the local
apparent viscosity in the lubrication force only develops a
weak solids volume fraction wave, the waves with either
/u_ or p0 in the lubrication expression are much stronger.
For the case with it0 this was expected: Since u, is the
highest viscosity in the system, lubrication makes it
relatively hard for spheres to approach and separate. Why
also the case with iu_ (the lowest viscosity in the entire
system) develops a stronger wave than the one with the
local viscosity is not at all clear, and is subject of current
investigations.
0.4 0.6 0.8 1 1.2
2I I I I I
320
240 
240 
t =7P
320 .
240 
a 160 X , _
160
80
20 40 60 80 100 120 140 160 180 200
tx 103
Figure 10. Solids volume fraction as a function of space
(the vertical coordinate x) and time. x and t are in lattice
units; x=320 corresponds to x=20d; t=2105 corresponds to
t=225 d/u ). From top to bottom: lubrication based on local
apparent viscosity; based on p ; based on
/0 = 'p (1+ a).
Conclusions
Direct simulations based on the latticeBoltzmann method
of solidliquid systems with the liquid having
nonNewtonian rheology have been performed. We
specifically looked into purely viscous thixotropic (i.e.
timedependent) liquids that show shearthinning behavior.
We first showed the impact of shearthinning on the drag
force on spheres in random, dense assemblies. The effect of
shearthinning is most pronounced if the dimensionless
number c (with S a measure for the interparticle
U,
spacing, and yk a characteristic shear rate of the liquid) is
of order 1 so that significant shear thinning occurs in the
Paper No
waists between nearby spheres. Thixotropy tends to increase
the drag on spheres in dense assemblies.
We also (preliminary) studied hindered settling with
particles moving relative to the liquid and relative with
respect to one another. Lubrication forces that we need to
add explicitly to the simulation procedure in order to
account for spheres in very close proximity have strong
impact on the way the sedimenting system develops
instabilities. The role of lubrication forces and their
generalization to nonNewtonian fluids therefore is an
important subject of further study.
Acknowledgements
Support from the Canadian Oilsands Network Research And
Development (CONRAD) within their Bitumen Production
Fundamental Research Group is gratefully acknowledged.
References
Chen, S., Doolen, G.D. Lattice Boltzmann method for fluid
flows. Annual Rev. Fluid Mech., Vol. 30, 329 (1998)
Derksen, J. & Van den Akker H.E.A. Largeeddy
simulations on the flow driven by a Rushton turbine. AIChE
J., Vol. 45, 209 (1999)
Derksen, J.J. & Sundaresan, S. Direct numerical simulations
of dense suspensions: wave instabilities in liquidfluidized
beds. J. Fluid Mech., Vol. 587, 303 (2007)
Derksen, J.J. & Prashant. Simulations of complex flow of
thixotropic liquids. J. NonNewtonian Fluid Mech., Vol. 160,
65 (2009)
Duru, P., Nicolas, M., Hinch, J. & Guazelli, E. Constitutive
laws in liquidfluidized beds. J. Fluid Mech., Vol. 452, 371
(2002)
Ferroir, T., Huynh, H.T., Chateau, X. & Coussot, P. Motion
of a solid object through a pasty (thixotropic) fluid. Phys.
Fluids, Vol. 16, 594 (2004)
Hill, R.J., Koch, D.L. & Ladd, A.J.C.
ModerateReynoldsnumber flows in ordered and random
arrays of spheres. J. Fluid Mech., Vol. 448, 243, (2001)
Kandhai, D., Derksen, J.J. & Van den Akker, H.E.A.
Interphase drag coefficients in gassolid flows. AIChE J.,
Vol 49,1060(2003)
Kim, S. & Karrila, S. J. Microhydrodynamics: Principles
and Selected Applications, ButterworthHeinemann (1991)
Ladd, A.J.C. Numerical simulations of particle suspensions
via a discretized Boltzmann equation. Part 2. Numerical
results. J. Fluid Mech., Vol. 271, 311 (1994)
Masliyah, J. Zhou, Z.J.,. Xu, Z., Czarnecki, J. & Hamza, H.
Understanding waterbased bitumen extraction from
Athabasca oil sands. Can. J. Chem. Engng., Vol 82, 628
(2004)
Moore, F. The theology of ceramic slips and bodies. Trans.
71" International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Br. Ceramic Soc., Vol. 58, 470 (1959)
Nguyen, N.Q. & Ladd, A.J.C. Lubrication corrections for
latticeBoltzmann simulations of particle suspensions. Phys.
Rev. E, Vol. 66, 046708 (2002)
Sangani, A.S. & Acrivos, A. Slow flow through a periodic
array of spheres. Int. J. Multiphase Flow, Vol. 8, 343 (1982)
Somers, J. A. Direct simulation of fluid flow with cellular
automata and the latticeBoltzmann equation. App. Sci. Res.,
Vol. 51, 127 (1993)
Storey, B.T. & Merrill, E.W. The rheology of aqueous
solution of amylose and amylopectine with reference to
molecular configuration and intermolecular association. J.
Polym. Sci., Vol. 33, 361 (1958)
Succi, S. The lattice Boltzmann equation for fluid dynamics
and beyond, Clarendon Press, Oxford (2001)
Sweby, PK. High resolution schemes using flux limiters for
hyperbolic conservation laws. SIAM J. Numer. Anal., Vol.
21,995(1984)
Ten Cate, A., Nieuwstad, C.H., Derksen, J.J. & Van den
Akker, H.E.A. PIV experiments and latticeBoltzmann
simulations on a single sphere settling under gravity. Phys.
Fluids, Vol. 14, 4012 (2002)
Van der Hoef, M.A., Beetstra, R. & Kuipers, J.A.M.
LatticeBoltzmann simulations of lowReynoldsnumber
flow past mono and bidisperse arrays of spheres: results for
the permeability and drag force. J. Fluid Mech., Vol. 528,
233 (2005)
Wang, J.W. & Ge, W. Collisional particle phase pressure in
particlefluid flows at high particle inertia. Phys. Fluids, Vol.
17, 128103(2005)
