7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Analyses of rotational diffusion of rigid fibers and orientation effect on rheology of
flexible fibers in shear flow with direct numerical simulation.
Asif Salahuddin and Jingshu Wu and C.K. Aidun
Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, GA 30318, USA
gth823e@mail.gatech.edu and wujingshu@gatech.edu and cyrus.aidun@me.gatech.edu
Keywords: latticeBoltzmann method, rotational diffusion, rigid and flexible fiber suspension
Abstract
The recent developed particle level simulation method (Wu & Aidun (2010a,b)) is implemented in this study to
investigate the orientation diffusion of rigid and the orientation and rheology of flexible rodlike fibers in semidilute
suspension of creeping shear flow, with a Newtonian fluid medium. For rigid fibers, an anisotropic, weak rotary
diffusion model proposed by Rahnama et al. (1995) is tested to evaluate the rotational diffusivity ratio, R () ,
from a one parameter fit of the model to the simulated orbitconstant Cbdistribution in semidilute regime. The
diffusivity ratios from simulation data compare very well with Stover et al. (1992)'s experiments and Rahnama et
al. (1995)'s hydrodynamicc interaction theory'. But, the simple anisotropic weak diffusion model can not describe
the violation of Stokes flow symmetry which is triggered by the presence of a small but detectable amount of
nonhydrodynamic (mechanical) fiberfiber interaction in creeping flow in semidilute regime and thus emphasizes the
need for a more sophisticated model to predict orientation diffusion. For flexible fibers, the decrease of fiber stiffness
(bending ratio BR) is shown to cause increasing asymmetry in pdistribution (here, p is the meridian angle in the
flowgradient plane). The impact of this asymmetry on first normal stress difference, NB is discussed in the light of
Batchelor (1971)'s theory.
Introduction
Fiber Orientation in Suspension
The understanding of orientation distribution and mi
crostructure of rigid and flexible fiber suspensions are
very important in papermaking and composite produc
tion. The wellknown Jeffery (1922)'s solution of the
motion of a single ellipsoid in suspension completely ne
glected effects such as fluid and particle inertia, Brown
ian motion and particleparticle interactions, all of which
may be present in industrially important semidilute to
concentrated suspension systems. By virtue of their abil
ity to simulate these effects in suspensions, particlelevel
simulations are useful and serve as complements to the
oretical and experimental approaches in research.
In the past, some notable investigations on diffusiv
ity of suspensions have mostly been performed (An
czurowski & Mason (1967a,b)) at concentrations lower
than aimed to be studied in this research. Also, the orien
tation diffusivity was considered to be isotropic (Folgar
& Tucker (1984)).
Jeffery (1922) showed that in simple shear flow, a single
fiber rotates indefinitely about the vorticity axis along
one of an infinite number of periodic, closed orbits.
For almost any rigid body of revolution including cir
cular cylinders, the orbit period (T) will be: T
27 (r + 1/re) /j, where j is the shear rate and a semi
empirical relation (Cox (1971)) replaces fiber aspect
ratio rp L/d with an effective aspectratio, re
1.24r (lnrp)1/2.
The spherical coordinate system used to describe
fiber orientation is defined in Fig. 1. The integration of
time evolution for 0 and p yields:
tan 6 2 1/
(r6 cos2 + sin2 9)1/
tan = re tan (27 + K)
(T
where K is the phase angle. The orbits of the ends of the
fiber are a symmetrical pair of spherical ellipses defined
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Voricity
Flow
Figure 1: The coordinate system with a fiber centered
at the origin. The unit vector, p is parallel to the fiber's
axis. Typical orbits for a slightly prolate spheroid (fiber)
are shown with solid blue lines.
C = tan00 C os2 0o + r sin2 yo
v r2
(3)
where, C is the orbit constant and subscript '0' denotes
the initial orientation. The value C 0 corresponds
to perfect alignment of the fiber with the vorticity di
rection of the flow, whereas C = o corresponds to
rotation in the flow gradient (xy) plane. It is advan
tageous to characterize the orbits with an orbit constant,
SC/(1 + C) [0, 1].
Fiber suspensions are broadly divided into three
regimes: in the dilute regime, nL3 << 1, the semidi
lute regime is defined by nL3 >> 1 but nL2d << 1.
And finally, in concentrated regime c, > d/L, (c
n7dd2L/4 for cylindrical fibers). Here, n is the fiber
number density, d is the fiber diameter and L is the fiber
length. The fiber bending ratio (BR) is a measure of
nondimensional stiffness of the flexible fiber defined by
Forgacs & Mason (1959) as:
Ey (In 2r, 1.5)
BR = (4)
2(/,)r
Here, Ey is the fiber's Young's modulus. It is to be
noted that, the fibers used in this study are noncolloidal
with P6clet number, Pe oc, so the effect of Brown
ian motion is negligible. The particle Reynolds number,
Re = L2 /v, is very small (the fibers are on the order
of microns) and therefore, the fiber inertia is also negli
gible.
Rotary Diffusivity
As a simplification to the detailed accounting of fiber
fiber interaction for rodlike fibers, the small changes in
orientation of one fiber because of the presence of it's
neighbors can be accounted for with a rotary diffusion.
Folgar & Tucker (1984) introduced a phenomenological
model which describes the effects of fiber interactions
in terms of an isotropic, rotary diffusivity, D,. Leal
& Hinch (1971) tried to develop a solution for steady
state Cdistribution by accounting for a weak scalar
rotary diffusivity but failed to reproduce Anczurowski
& Mason (1967b)'s experimental Cdistribution. Rah
nama et al. (1995) generalized Leal & Hinch (1971)'s
solution of the Burgers evolution equation by including
the orientation dependent anisotropic diffusivity tensor,
D,. So, the differential probability distribution function
2 (0, y, t) can be stated as:
Q +V, (pJ ) Vp.(D, VQ) V (Q2ph) (5)
where, Vp is the gradient operator in orientation space,
pij is the rotation rate of the fiber as described by Jef
fery (1922), and ph is the drift velocity. The original
eqn. (5) was reexamined to consider the consequences
of the action of a small amount of weak anisotropic dif
fusion in the limit as time, t oc and drift velocity
was neglected. In the (C, T) orbit coordinates (see Leal
& Hinch (1971) for this coordinate definition), where
T is the phase angle, the orientation distribution can be
separated into two parts:
2 (C, ) f (C) (C, T) (6)
In this expression f (C) is the unknown distribution
function describing the population amongst the Jeffery
orbits; g (C, T) is the distribution around an individual
orbit, C. In the large aspectratio limit, the integral
expression for f (C) was analytically solved to yield:
C) 4RC where R=( D) Do is the
7r(4RC2+1)3/2 D "
proportionality between gradients of 2 (0, y) in the 0
direction and the flux of probability in the 0direction.
A large value of R implies that the fibers are aligned
near the vorticity axis, as the effect of Doo is to push
fibers away from the flow direction and toward vortic
ity or toward decreasing values of orbit constant. The
f (C) can be transformed to give the differential proba
bility distribution function, p (Cb) for steady state orbit
constant distribution as:
p(Cb)
4RCb
(4R l Cb)]2 1 3/2
,4R [Cb/( C)+ (1
Ob)
(7)
where fo p (Cb) dCb = 1. Eqn. (7) would be referred to
as the 'Anisotropic Diffusivity Model'.
LBMEBF Method
In this study, we use a novel particlelevel numerical
simulation method (Wu & Aidun (2010a,b)) to analyze
suspension of both rigid and flexible fibers. The im
portant aspects of the implemented LBMEBF method
are: simulations include fiberfiber contact and lubrica
tion forces, the fluid and fibers are twoway coupled with
direct numerical simulation and the physical properties
of flexible fibers are directly related to simulation pa
rameters. In this method, the noslip boundary condition
at the fluidsolid interface is based on external boundary
force (EBF) method. The details of this method are pre
sented in previous papers (Wu & Aidun (2010a,b)).
Results and Discussion
Validation of the latticeBoltzmann approach with the
EBF method for fiber suspension is presented in previ
ous publications (Wu & Aidun (2010a,b)). The focus
of this paper is to investigate the anisotropic diffusiv
ity model Eqn. (7) for large aspectratio rigid fibers in
semidilute regime. The present report would also dis
cuss the effect of fiber stiffness on rheology and orienta
tion of flexible fiber suspension to some extent.
An unbounded shear domain is implemented based
on the LeesEdwards boundary condition (LEbc) (Lees
& Edwards (1972)) to improve computational efficiency
by removing wall effects. The computational domain is
5L x 5L x 4L and the suspending fibers have diameter
of d = 0.4 LBM unit lattice size.
Orientation Diffusivity of Rigid Fibers in
Suspension
The Cb at each timestep is divided into nno. of bins,
n
where ZdCb(i) 1.0. A MATLAB function dis
i= 1
tributes the fibers into i no. of bins based upon their
corresponding Cb(i) values. Then, the probability of
finding a specific fiber in ith bin is calculated as:
no. of fibers in Cb() ri / \ i
P Cb) ttal no of fibers The p ( b) is normal
ized with the integral area under the histogram. Sim
ilarly, the probability distribution, p (y), is calculated,
n
where, Jdo(i) 7r.
i= 1
Stover et al. (1992) experimentally measured rigid
fiber orientation in a cylindrical couette device using
an isorefractive suspension of matched densities, with
few observable opaque test fibers. They found that the
semidilute Cbdistributions were more uniformly dis
tributed than the dilute Cbdistributions of the exper
iments conducted by Anczurowski & Mason (1967b).
Anczurowski and Mason's experiments (in cylindri
cal counterrotating couette device) in infinitely dilute
regime heavily favored lower orbit constants. Accord
ing to Stover et al. (1992), the distribution of Cb in the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
nL3 Stover's Hydrodyn. LBMEBF
experiment Interaction simulation
Theory
rp 31.9 rp 31.9 r, 32.0
45 3.16 2.55 2.58
18 1.38 3.12 3.31
10 2.85 3.50 2.86
5 2.18 3.50 2.45
Table 1: Comparison of DiffusivityRatio, R (")
steady state has a peak in the range of 0 < Cb < 0.5
for both dilute and concentrated systems. Yamane et al.
(1994)'s numerical simulations (with only shortrange
hydrodynamic interaction) and Fan et al. (1998)'s nu
merical simulations (with both the short and longrange
hydrodynamic interactions) showed ili.ii sc.idv state Cb
distribution shifted towards higher values (0.42 ~ C ~
1.12 or 0.296 ~ Cb ~ 0.528) than observed in Stover
et al. (1992)'s experiments. The reason of this discrep
ancy is probably due to statistical errors in Fan et al.
(1998)'s numerical simulation induced by the choice of
smallsized domain (0.5 < 1 < 1.0, where, 12 is the do
main height in gradient direction) to save computational
cost which rules out the possibility of fiber orientation
perpendicular to the (xy) plane.
In this study, a set of LBMEBF simulations are
performed at nondimensional volume concentrations,
nL3= 5, 10, 18 and 45 with rigid fibers of rp=32.0. The
Cbdistributions with the LBMEBF simulations for this
nL3 range in semidilute regime demonstrate good pre
dictions of the Stover et al. (1992) experiments as can
be noted in Fig. 2. The Cbdistributions in Fig. 2 for
LBMEBF are smoothed over the last two orbit periods
simulated. The peaks of the Cbdistributions remain in
the range 0.15 < Cb < 0.4 for LBMEBF which is con
sistent with the report of Stover et al. (1992) experiments
in semidilute regime.
To determine the diffusivity ratio R= ) in
semidilute regime, Stover et al. (1992) fitted the
'anisotropic diffusivity model' (eqn. (7)) to their ex
perimental Cbdistribution with R as an adjustable pa
rameter. They also fitted the 'anisotropic diffusivity
model' to Anczurowski & Mason (1967b)'s dilute ex
perimental Cbdistribution. The fibers with rp=18.4
for Anczurowski & Mason (1967b)'s experiment show
a plateau at about ( ) =17 for the apparent dilute
regime (nL3 =0.016, 0.066). But, a steep drop to another
plateau, at about ( 2,,) =1.5 was observed for Stover et
al. (1992)'s semidilute regime experiments.
Rahnama et al. (1995) theoretically studied the ori
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
S 01 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Cb
SLBMEBF simulation
Best fitted 'anisoopic diffusivity
2 del'to BMEBPF imlim dta
x with R=2.86
x Stover et. al. (1992) experiment
x r =32.0
'
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Cb
25
x vert. al. (1992) experiment
SLBMEBF imulation
Best fitned 'aniso*tre cdifffaivity
0C 01 0.2 03 04 05 06 0.7 0.8 0.9
2 mdel'tIBMEBF im d&ta
1 .5.0
S 0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9
x SW .a. (1992) experen t
*LBMEBF simulation
2 Best fitted *amtsopop diffuvity
mdel'tomMF imlaim dita
IBrM=32.0
X
0 0.1 02 0.3 0.4 Cb0.5 0.6 0.7 0.8 0.9
Figure 2: p(Cb) with rp=32.0 for nL3=5, 10, 18 and 45
from top to bottom subfigure respectively. The exper
imental Cbdistributions are within 80% confidence in
terval. The solid curves demonstrate the oneparameter
fit (with parameter R) of 'anisotropic diffusivity model'
to LBMEBF simulation data. The bestfitted R values
are indicated in corresponding legends.
x Stovetal. (1992) expiment
* LBMEBF simulation
Best fitted 'anisotropic diffusivity
mdel'tolBMEBF aimdim ta
tth tR2. 4
=32.O
entation diffusivity in dilute to semidilute regime with
hydrodynamicc interaction theory'. In this theory hy
drodynamic, orientation diffusivity was obtained from
an ensemble average of the fiberfiber interactions. The
steadystate fiber orientation distribution is controlled
by the anisotropy and orientation dependence of the
diffusivity. The steady state and transient fiber ori
entation distributions are derived using a perturbation
analysis for weak hydrodynamic orientation diffusion.
For computational convenience, Rahnama et al. fitted
the 'anisotropic diffusivity model' to the Cbdistribution
calculated with hydrodynamicc interaction theory' and
obtained the bestfit value of R. Rahnama et al. used
an iterative solution with an initial guess of R (averaged
R=2.4 from Stover et al.'s experiments for rp=31.9 was
used as an initial guess) to calculate the Cbdistribution
with hydrodynamicc interaction theory'. We also fitted
the 'anisotropic diffusivity model' to the Cbdistribution
from LBMEBF simulation and predicted the bestfit R
values in semidilute regime. A nonlinear curve fitting
method is used for this purpose.
Table 1 summarizes R=() D values by Stover et
al. (1992) experiments, hydrodynamicc interaction the
ory' and LBMEBF simulations. The R values in the
semidilute regime from nL3=145 fall in the same range
for theory and experiment. The quantitative compari
son between LBMEBF simulation and experiments are
good, considering statistical uncertainties in the experi
ment. However, Stover et al. (1992)'s experimental val
ues of R do not show any systematic dependence on
the volume concentration, nL3 and aspect ratio, rp.
Whereas, the hydrodynamicc interaction theory' reveals
a dependence of diffusivity ratio, R, on nL3 and as
pect ratio, rp: for a fixed aspect ratio (rp=31.9), the
value of R decreases with increasing nL3. Physically,
this means that the fibers shift closer toward the flow di
rection with increasing nL3. According to Rahnama et
al. (1995), this shift resulted from the anisotropic hy
drodynamic screening incorporated in the renormalized
Green's function derived by Shaqfeh and Fredrickson
(1990). The LBMEBF predicted R values first increase
slightly from nL3=518 and then it decreases with in
creasing concentration following the theoretical predic
tion. The physical reason behind this can be explained
by studying change of (Cb) values with concentration.
From Fig. 3 it is seen that, the (Cb) for Stover et
al. (1992) experiments and LBMEBF simulations de
crease with increase of nL2d (the parameter nL2d is
used since, it is a more relevant measure of concen
tration in semidilute regime). Physically this means
the fibers are shifting towards the vorticity axis and
as a result R value also increases (diffusion occurring
primarily in flowvorticity plane) with concentration.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
C Stovr t al epetnment (1992), r31 9
LBM EBF smulation, r =32 0
Hydrodynac minteraction theory
Figure 3: Average value of orbit constant, (Cb) as a
function of nL2 d. The prediction of hydrodynamicc in
teraction theory' is indicated with a solid line.
But, as the concentration increases more (nL3>18 or
nL2d > 0.56253)), the fiberfiber hydrodynamic inter
actions dominate and also there can be presence of some
mechanical contacts among fibers. Both of these effects
cause extra fiber flipping and as a result the fibers move
closer to the flow direction ((Cb) increases) and the dif
fusivity ratio R decreases (i.e. diffusion occurring pri
marily in the flowgradient plane). So, the (Cb) values
of Stover et al. (1992)'s experiment, hydrodynamicc in
teraction theory' and LBMEBF simulation fall close to
each other for higher values of concentration, nL3 45
(nL2d z 1.41). Now, the apparent discrepancy be
tween theory and experiment at low concentrations can
be due to the lower accuracy of renormalized Green's
function at low concentrations (which is used to account
for hydrodynamic screening in hydrodynamicc interac
tion theory'). Actually the function shows better accu
racy at higher values of nL3 (Shaqfeh and Fredrickson
(1990)). Perhaps at lower concentration, a twofiber the
ory would work better. On the flipside, through a pri
vate communication with Dr. Koch, the authors came
to know that during Stover et al. (1992)'s experiments,
the group of researchers did not to pay too much atten
tion to these trends of diffusivity ratio, R not be
ing too sure if those were statistically ,igiilik.ii and
rather paid more attention just to the overall order of
magnitude of diffusivity ratio R. So, we can strongly
claim that, LBMEBF simulations produced excellent
results to verify the range of R in semidilute regime
and also shows encouraging results to physically ex
plain the trend of R in semidilute suspension of non
colloidal rigid fibers. The LBMEBF simulation has the
algorithm to count the number of mechanical contacts
present in the suspension during a flow situation. We
observed a small amount of mechanical contacts in the
semidilute regime which breaks the Stokes flow symme
try (results not shown here for conciseness). This phe
nomenon was also observed in Stover et al. (1992) and
Petrich et al. (2000)'s semidilute regime experiments
r,=16
c,=.o053
nL=17.3
BR=2940
 BR=1.47
*BR=0.15
0.2 0.4 0.6 0.8
~In
0 0.2 0.4 0.6 0.8
tIu
Figure 4: The pdistribution for flexible fiber suspen
sions with aspect ratio, r = 16 (top) and r = 32 (bot
tom) with different bending ratios, BR, and for a fixed
volume fraction, c, = 0.053 with LBMEBF simula
tion.
where the pdistribution showed asymmetry across the
flowvorticity plane due to nonhydrodynamic mechan
ical contacts among fibers and thus violated the Stokes
flow symmetry. This phenomenon can not be explained
by the weak diffusion considered in 'anisotropic diffu
sivity model' (Eqn. (7)) used to predict diffusivity ra
tio, R, since weak diffusion should preserve the symme
try of Stokes flow. This indicates that a complete de
scription of hydrodynamic interactions is more complex
and sophisticated than the simple 'anisotropic diffusiv
ity model' used to fit the orbitconstant distribution. Al
though the model was very good fit to the p (Cb) it can
not describe the asymmetry in p ( ).
Rheology and Orientation of Flexible Fibers in
Suspension
In this section, the change in fiber theological properties
with change of fiber flexibility is explained by study
ing the pdistribution, p (4) of the fiber suspension for
different aspectratios, rp and bending ratio BR and
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
those information are theoretically related to Batchelor
(1971)'s work.
Batchelor (1971)'s relation for the contribution of sus
pended fibers in dilute suspension without any Brownian
motion, yields the first normal stress difference, N1 as:
N B = jfiberf (KPI. )
fiberr (sin 4 sin 4 ))
*(<(si Ism14p))
Where, the superscript 'B' refers to Batchelor's the
ory and fiber is a function of concentration, orienta
tion distribution and fiber geometry. Eqn. (8) shows
that if the suspension has no direct physical contacts
between fibers and in the absence of fiberdeformation
(rigid fiber), the fiber orientation dispersion, p (p) would
necessarily be symmetric about (xz) plane. And in these
cases, first normal stress difference, Nf given by Batch
elor's theory would be zero, since it is an odd function of
/. So, if direct contact between fibers exists or if fibers
are deformable, Nf will not vanish.
The pdistribution in Fig. 4 with LBMEBF simula
tion clearly gives evidence that for decreasing bending
ratio, BR (i.e. by increasing fiber flexibility), the mean
orientation angle (p) will become slightly less than 7r/2;
and this small asymmetry of the fiber orientation distri
bution makes Nf nonzero (+ve). So, by reviewing the
Batchelor's theory above and by giving computational
evidence of asymmetry in pdistribution due to change
of BR, here it is proved that the first normal stress differ
ence is strongly dependent on the fiber flexibility. On the
other hand, from Eqn. (8), it is clear that the change of
0 is also very important for the first normal stress differ
ence. NA increases with 0, where 0 is directly related to
the orbit constant Cb (Eqn. (1)). But, the concept of orbit
constant is not welldefined for flexible fibers (Skjetne
et al. (1997)), since the geometry of the fiber vary with
time. For flexible fibers Jeffery's orbit become unstable
and C will tend to drift (for the most part) either to 0
or oo as evidenced in experiments (Arlov et al. (1958))
and numerical simulations (Skjetne et al. (1997)), de
pending on the initial C value, while intermediate values
also being observed. To the authors' knowledge, no ro
tational diffusion model exists which accounts for these
disturbed Jeffery's orbits for flexible fibers.
It is worth to point out here that, the effect of vol
ume fraction on relative shear viscosity of suspension
p,,e of flexible fibers has been discussed in previous Wu
& Aidun (2010a) paper with LBMEBF simulations. It
was shown that with increase of volume fraction the rel
ative shear viscosity, pre increases which can be at
tributed to the increased fiberfiber interaction due to
fiber crowding and fiber deformation with shear strain.
Conclusion
(. I /*' ))
The steadystate distribution of orbitconstant, p (Cb), of
rigid cylindrical fibers suspended in a Newtonian fluid
subjected to simple shear flow has been predicted for
several concentrations (nL3=545) in semidilute regime
with LBMEBF simulation. The 'anisotropic diffusivity
model' is fitted to the simulated p (Cb) which gives best
fit diffusivity ratios (R values) in the range of 2.45 ~
R ~ 3.31 and thus verifies the range of R observed in
experiment and theory. But, the small amount of asym
metry in p (4) as was observed in simulation and experi
ments questions the simplicity of 'anisotropic diffusivity
model' and emphasizes the need for the development of
a more sophisticated model. Also, in future, it would
be interesting to numerically investigate the time corre
lation function of Cb(t), defined as (Cb(t)Cb(t T)),
where, T is the delay time, t is the time measured dur
ing the simulation, and the angle brackets denote aver
ages over all available values. By analyzing the decay
of (Cb(t)Cb(t T)) as a function of T, the temporal
stochastic fluctuations of Cb(t) can be explored and it
would be a good test of the magnitude of hydrodynamic
rotary diffusivity.
The LBMEBF simulation with flexible fibers shows
fiber stiffness has strong impact on suspension rheology.
Also, the effect of fiber stiffness on the first normal stress
difference, NB is demonstrated based on analyzing the
pdistribution and Batchelor (1971)'s theory. It is proved
that the influence of fiber orientation due to change in
fiber stiffness is a major contributor to the variation in
theological properties.
References
Anczurowski, E. and Mason, S.G., The kinetics of flow
ing dispersions: II. Equilibrium orientations of rods and
discs (theoretical), Journal of Colloid and Interface Sci
ence, Vol. 23, no. 4, pp. 522532, 1967a
Anczurowski, E. and Mason, S.G., The kinetics of flow
ing dispersions: III. Equilibrium orientations of rods and
discs (experimental), Journal of Colloid and Interface
Science, Vol. 23, no. 4, pp. 533 546, 1967b
Arlov, A.P. and Forgacs, O.L. and Mason, S.G., Particle
motions in sheared suspensions IV. General behavior of
wood pulp fibers, Svensk Papperstidn., Vol. 61, no. 3,
pp. 6167, 1958
Batchelor, G.K., The stress generated in a nondilute
suspension of elongated particles by pure straining mo
tion, Journal of Fluid Mechanics, Vol. 46, no. 04, pp.
813829, 1971
Cox, R.G., The motion of long slender bodies in a vis
cous fluid. Part 2. Shear flow, Journal of Fluid Mechan
ics, Vol. 45, no. 04, pp. 625657, 1971
Fan, X.J. and PhanThien, N. and Zheng, R., A di
rect simulation of fibre suspensions, Journal of Non
Newtonian Fluid Mechanics, Vol. 74, no. 13, pp. 113
135, 1998
Folgar, F and Tucker III, C.L., Orientation Behavior of
Fibers in Concentrated Suspensions, Journal of Rein
forced Plastics and Composites, Vol. 3, no. 2, pp. 98 
119, 1984
Forgacs, O.L. and Mason, S.G., Particle motions in
sheared suspensions. IX. spin and deformation of thread
like particles, Journal of Colloid Science, Vol. 14, pp.
457472, 1959
Jeffery, G.B., The motion of ellipsoidal particles im
mersed in a viscous Fluid, Proceedings of the Royal So
ciety of London. Series A, Containing Papers of a Math
ematical and Physical ( li.In.ikic Vol. 102, no. 715, pp.
61179, Nov. 1, 1922
Leal, L.G. and Hinch, E.J., The effect of weak Brown
ian rotations on particles in shear flow, Journal of Fluid
Mechanics, Vol. 46, no. 4, pp. 685703, 1971
Lees, A.W. and Edwards, S.F, The computer study of
transport processes under extreme conditions, Journal of
Physics C: Solid State Physics, Vol. 5,no. 15 pp. 1921
1928, 1972
Petrich, M.P. and Koch, D.L. and Cohen, C., An exper
imental determination of the stressmicrostructure rela
tionship in semiconcentrated fiber suspensions, Journal
of NonNewtonian Fluid Mechanics, Vol. 95,no. 23, pp.
101133, 1994
Rahnama, M. and Koch, D.L. and Shaqfeh, E.S.G., The
effect of hydrodynamic interactions on the orientation
distribution in a fiber suspension subject to simple shear
flow, Physics of Fluids, Vol. 7, no. 3, pp. 487506, 1995
Shaqfeh, E.S.G. and Fredrickson, G.H., The hydrody
namic stress in a suspension of rods, Physics of Fluids
A: Fluid Dynamics, Vol. 2,no. 1, pp. 724, 1990
Skjetne, P. and Ross, R.F and Klingenberg,
D.J.,Simulation of single fiber dynamics, Journal
of Chemical Physics, Vol. 107,no. 6 pp. 21082121,
1997
Stover, C.A. and Koch, D.L. and Cohen, C., Observa
tions of fibre orientation in simple shear flow of semi
dilute suspensions, Journal of Fluid Mechanics, Vol.
238, pp. 277296, 1992
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Wu, J. and Aidun, C.K., A method for direct simulation
of flexible fiber suspensions using latticeBoltzmann
equation with external boundary force, International
Journal of Multiphase Flow, Vol. 36, pp. 202209, 2010a
Wu, J. and Aidun, C.K., Simulating 3D deformable par
ticle suspensions using lattice Boltzmann method with
discrete external boundary force, International Journal
for Numerical Methods in Fluids, Vol. 62, pp. 765783,
2010b
Yamane, Y. and Kaneda, Y. and Dio, M., Numerical sim
ulation of semidilute suspensions of rodlike particles in
shear flow, Journal of NonNewtonian Fluid Mechanics,
Vol. 54, pp. 405421, 1994
