Em:

‘fafi.

Fun-14”!

 

.‘.-~4
, - . .
1-..-

u.
I! "' ‘3
333
-.«
.cv
*

 

'3 Ag)“.§§a ’4 '3‘ " I ‘ I ' . I 7H “ . laa‘txfl‘

A :‘. ‘ s ‘ > . u ' ' v; T:
.45 55;! 2' . ‘: .A , _ A 2;: _ ‘ 2,.
5: - .'

hi :5"

1

g, i
I} a
f???

If 7
2:45
. ’if 8

 

100?

This is to certify that the
dissertation entitled

FLOW-INDUCED ALIGNMENT AND MIGRATION OF
PARTICELS IN SUSPENSIONS

presented by

LIPING JIA

has been accepted towards fulfillment
of the requirements for the

PhD degree in Mechanical Engineering

 

 

, \ 1 Major ProfessoiE’Signature
\
419 4‘] t 9 ; r200 <4“
,4

I“ “I Date

MS U is an Affirmative Action/Equal Opportunity institution

. ._.-.-u-c-o-oOI-.-u-v-3-.---0-I—o-I--I-c-0-o-.-----.--o-o---u-u-o-I.

-.—¢--I-l-O-‘-‘-

 

LIBRARY
Michigan State
University

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/07 p:ICIRC/DateDue.indd-p.1

FLOW-INDUCED ALIGNMENT AND MIGRATION OF PARTICLES IN
SUSPENSIONS

By

Lip‘ing Jia

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

2006

ABSTRACT
FLOW—INDUCED ALIGNMENT AND MIGRATION OF PARTICLES IN
SUSPENSIONS

By
Liping Jia

The alignment and migration of suspensions are important for industrial processes as-
sociated with composite processing, the fabrication of microelectronics devices, the manu-
facturing of products with micro- and nano-scale suspensions, the environment (pollutants
migration, particulates, and microbes), and the petroleum industry. In this work, problems
associated with the motion of a single particle are solved and models needed to describe
the orientation and migration of a large number of particles are developed;

The hydrodynamics of a single ellipsoid suspended in an unbounded homogeneous flow
was first investigated by Jeffery in 1922 [1]. Jefl'ery’s work deals with the problem of in-
compressible Newtonian fluid with constant viscosity, no-slip on the interface of solid/fluid,
and linear shear flow. Based on Brenner’s asymptotic method [2] analytical solutions are
developed to study the influence of other conditions on the motion of a single particle,
i. g. slip boundary conditions on the interface, other flow fields (a quadratic flow and
cubic flow), and viscosity. The results are partially validated by comparing with existing
solutions for some limiting cases of no-slip, perfect slip, sphere, and constant viscosity.
Equations describing the motion of a single particle under different conditions are derived.
A different method is used to study the influence of inertia forces on the motion of a single

particle, which is based on Burgess’ general solutions [3] of a viscous Oseen flow. Differ-

ent velocity fields of the fluid are found for the cases of translation motion of a sphere and
a deformed sphere with slip and no-slip boundary condition.

Equations describing the motion of ensembles of rigid particles of complex shapes are
studied next. Each particle is assumed to be non-axisymmetric, and its orientation is de-
scribed with three Euler angles. The geometry of such particles (e.g. ellipsoids) and their
interactions with the surrounding fluid are described by a third order tensor instead of the
single parameter often used for axisymmetric particles (spheroids). To compute the flow-
induced alignment of these particles, one must solve an evolution equation for the orienta-
tion distribution function but such computations are costly. Instead, an evolution equation
for the second moment of the distribution function, which forms a fourth order tensor, is
used in order to obtain the average orientation of the particles in homogeneous flows. A
closure model is introduced for the unclosed eighth-order tensor which satisfies six-fold
symmetry and six-fold projection properties.

In the last part of this work, models describing solid-liquid two-phase flows are devel-
oped using a continuum approach. A so-called Eulerian-Eulerian technique is adopted to
deal with the motion of the non-spherical particles and Newtonian fluid. Based on the mo-
ments of the distribution function, the evolution of the second moment of the orientation
tensor is used to govern the orientation of particles statistically. The concept of control
volume/control surface method is used to develop closure models for the stresses and inter-
facial force. The fully symmetric quadratic model (developed for axisymmetric particles)
is applied to close the problem associated with computing the orientation tensor. A finite
element code is developed to simulate the alignment and migration of dilute suspensions

of spheroids in a flowing liquid.

Dedicate this work to my parents and my husband.

iv

ACKNOWLEDGMENTS

I must first thank my dissertation advisor Dr. André Bénard, who was never tired and
always with his kindly leading causes to find the way in my dissertation.

I am grateful for Dr. Charles Petty, who not only gave me advices, but also kept tracking
with all the detail of my research work.

I would like to thank to Dr. Alejandro Diaz, who helped me to join in my research
group and spent lots of time in reviewing my dissertation.

I want to acknowledge Dr. Peter Bates’ help for his serving my dissertation defense
committee, when he was extremely busy.

I would like thank my husband, Zhijian Huang, who accompanied me all the way
through this journey. Without him, I would not be able to finish the research on time.

I want to thank all my instructors at Michigan State University, who taught and encour-
aged me to reach this point. I would like to thank department staffs for their secretarial
support on my study and life at Michigan State University.

Also I want to acknowledge all my friends and family members who gave me support
through my years of study, enabling me to accomplish the work.

I gratefully acknowledge partial financial support of this work by the National Science
Foundation through the following education and research programs at Michigan State Uni-

versity: NSF/CTS 0083229, NSF/EEC-0331977, and Dissertation Completion Fellowship.

TABLE OF CONTENTS

LIST OF TABLES ix
LIST OF FIGURES x
1 Introduction 1
1.1 Background on alignment and migration of particles ........... 1
1.2 Dilute suspension theory .......................... 2
1.2.1 Hydrodynamics of a single particle ................ 2
1.2.2 Fundamental equations of creeping flows ............. 4
1.2.3 Particle material tensor ....................... 5

1.2.4 Jeffery’s solution of the motion of an axisymmetric particle in
a linear field ............................. 7
1.3 Non-Dilute Suspensions ........................... 11
1.4 Objectives of this work ........................... 12
2 MOTION OF A DEFORMED SPHERE WITH SLIP 15
2.1 Introduction ................................. 15
2.2 Basic formulations .............................. 18
2.3 Relation between pn, rpm and X" with X", Y" and Zn ............ 22
2.4 Stokes’ resistance of a spheroid ....................... 25
2.4.1 Functions to describe the shape of a deformed sphere ...... 25
2.4.2 Uniform streaming flow past a stationary deformed sphere . . . 26
2.5 Uniform streaming flow past a rotating deformed sphere ........ 28
2.6 Flow-induced motion of a spheroidal particle with slip .......... 30
2.7 Summary ................................... 47
3 MOTION OF AN ELLIPSOID IN QUADRATIC AND CUBIC FLOWS 48
3.1 Motion of an ellipsoid in a quadratic flow field .............. 48
3.2 Motion of an ellipsoid in a cubic flow field ................. 55
3.2.1 The hydrodynamic resistance ................... 55
3.2.2 Motion of a deformed sphere in a simple cubic flow ....... 57

vi

POWER-LAW MODEL OF A DEFORMED SPHERE 62

4.1 Introduction ................................. 62
4.2 Power-Law model for the Non-Newtonian viscosity ........... 62
4.3 Solution to flow problems using a Power-Law model ........... 65
MOTION OF A SPHEROID IN AN OSEEN FLOW 70
5.1 Introduction ................................. 70
5.2 Analytical structure of the Oseen flow ................... 71
5.3 Applications of Burgess’s solution ..................... 73
5.4 Summary ................................... 79
THE FLOW-INDUCED ORIENTATION OF RIGID PARTICLES IN DI-

LUTE SUSPENSIONS 80
6.1 Introduction ................................. 80
6.2 Prediction of orientation of axisymmetric particles ............ 82
6.3 Predictions of orientation of non-axisymmetric particles ......... 85
6.4 Algebraic restrictions on averaged orientation tensors .......... 88
6.5 Symmetry operator ............................. 90
6.6 Construction of the eighth-order orientation tensor ........... 92
6.7 Conclusions .................................. 95

PREDICTION OF FLOW-INDUCED ORIENTATION AND SPATIAL MI-

GRATION OF PARTICLES . 96
7.1 Introduction ................................. 96
7.2 Hydrodynamics of ensembles of particles ................. 99
7.2.1 Theory of ensemble averaging ................... 99
7.2.2 Averaged balanced equations .................... 101
7.3 Equations of motion and orientation for a dilute suspension ....... 102
7.4 Stress model ................................. 103
7.5 Interfacial force ............................... 106
7.6 Summery ................................... 109
SIMULATION OF THE FLOW INDUCED FIBER ORIENTATION AND MI-
GRATION USING A FINITE ELEMENT METHOD 110
8.1 Governing equations for 2-dimension problems .............. 110
8.1.1 Basic assumptions .......................... 110
8.1.2 Governing equations ........................ 112
8.1.3 Boundary conditions ........................ 118
8.2 Mixed finite element model ......................... 118
8.2.1 Weak form .............................. 118
8.2.2 Finite element model ........................ 120
8.3 Simulation of a plane Poiseuille flow .................... 123

vii

9 SUMMARY AND CONCLUSIONS

APPENDICES

A. Tensor notation used in this dissertation ..................
B. Matrix form of the weak form of the governing equations ........

BIBLIOGRAPHY

viii

132

136
136
138

148

2.1

3.1

LIST OF TABLES

Predicted periods of the particle motion induced by a simple shear flow.
The periods of the motion of the particle increase with the increasing of the
deformation of the particle and decrease with the increasing of the parame-
ter y ....................................... 39

Predicted periods of the particle motion induced by cubic flows and simple
shear flows. .................................. 60

ix

1.1

1.2

1.3

2.1

2.2

2.3

2.4

2.5

2.6

LIST OF FIGURES

Schematic of the coordinate system used to represent the orientation of a

particle. ....................................

An axisymmetric particle, i.g. a spheroid suspended in a simple shear flow.
The geometry of the particle is a spheroid and the aspect ratio of the particle

- -2
lSdp-a. ...................................

Jeffery orbits are shown for spheroidal particles at different aspect ratios a p
and different flow fields V.» = 7 ye x. The motions of the spheroidal particle
suspended in a simple shear flow are periodic and the periods depend on

the shape of the particle and the constant 7 of the flow field. ........

The slip length L s is defined for a simple shear flow in the presence of slip

boundary conditions at the interface of solids and liquids ...........

Illustration of lengths Ls used to describe different slip cases. When L5 = 0,
no slip occurs on the interface; when Ls = finitenumber, finite slip occurs

on the interface; when L S = Infinity, perfect slip occurs on the interface. . .

Illustration of three Euler angles ¢, 6, and ifi used to describe the coordinate
systems of a particle. x’, y’ and z’ are the reference coordinate system and

x, y and z are the rotating coordinated system. ................

Influence of the slip coefficient ,6 and the particle aspect ratio a on the
fluid/particle coupling coefficient A. xi is always positive when a < O and

negative when a > 0. .............................

Illustration of a simple shear flow surrounding a deformed sphere. .....

The influence of the fluid/particle coupling coefficient 11 on the rotation
of a spheroid in the flow/shear plane is shown by plotting the cosine of
the rotation angle a of a particle for various values of ,1. When lin < 1,
the induced motion of the particle is periodic; when Ill 2. 1, the induced

motion of the particle is steady. .......................

8

9

16

17

31

35

36

2.7

2.8

2.9

2.10

2.11

2.12

3.1

3.2

3.3

3.4

4.1

5.1

5.2

6.1

The figure shows the influence of the fluid/particle coupling coefficient xi
and the particle aspect ratio a on the temporal response of a spheroid to a
steady simple shear flow at positive values of as ................ 40

The above figure shows the influence of the fluid/particle coupling coefli-
cient and the particle aspect ratio on the temporal response of a spheroid to
a steady simple shear flow at negative values of as. ............. 41

These graphs show the velocity field around of a spheroid with steady mo-
tion. The particle reaches a special orientation after some time and then

keep this orientation inside the flow field. .................. 42
These graphs show the pressure field around of a spheroid with steady motion. 43
These graphs show the velocity field around of a spheroid with periodic

motion. .................................... 44
These graphs show the pressure field around of a spheroid with periodic

motion. .................................... 45
This figure shows a spherical coordinate system ................ 50
Illustration of a cubic flow field surrounding a deformed sphere. ...... 58

Evolution of the cosine of the rotation angle of the particle induced by
simple cubic flows at different constant Ks. ................. 59

Evolution of the cosine of the rotation angle of the particle induced by

simple shear flow at different constant Ks. .................. 60
Schematic of a typical viscosity variation w.r.t to the shear rate. ...... 64
Description of the surface velocity of the particle. .............. 74

Velocity vectors of the surrounding fluids with different boundary condition
on the surface of a particle. (a) no-slip boundary condition applying on a
sphere; (b) slip boundary condition applying on a spherew = 0.1); (c) no-
slip boundary condition applying on a deformed sphere(s = 0.2); (d) slip
boundary condition applying on a deformed sphere(s = 0.2, B = 0.1). . . . . 75

Description of the Euler angles used in this chapter. ............. 85

xi

6.2

8.1

8.2

8.3

8.4

8.5

8.6

Mapping procedure of a vector associated with the particle between the

reference configuration and the current configuration. ............

Quadrilateral elements used for the finite element model. (a) A nine-node
biquadratic element is used for the shape function of velocities. (b) A four-
node continuous-bilinear element is used for the shape function of the pres-

86

sure of fluids. ................................. 120

Domain and mesh for a plane Poiseuille flow with particle suspensions.
(a) Geometry, computational domain, and (b) the finite element mesh used
for the analysis of the slow flow with particle suspensions between parallel

plates. ..................................... 124

Contour plots of the principal eigenvalues 11mm of the orientation tensor
superposed with corresponding eigenvecotors for the problem of spheroids
suspended in a plane Poiseuille flow. The results are shown for three dif-

ferent times ................................... 127

Contour plots of the principal eigenvalues Tpmax of the particle stress su-
perposed with corresponding main eigenvecotors for the problem of spher-

oids suspended in a plane Poiseuille flow. .................. 128

Contour plots of concentration of the particles or for the problem of spher-

oids suspended in a plane Poiseuille flow ................... 129

Contour plots of the fluid pressure p f for the problem of spheroids sus-

pended in a plane Poiseuille flow ....................... 130

xii

CHAPTER 1

Introduction

1.1 Background on alignment and migration of particles

Alignment and migration of suspension of solids or droplets in a continuous medium at
low Reynolds numbers are important phenomena to various fields in associated with ma-
terial processing (composites), the fabrication of microelectronics devices, the manufac-
turing of products with suspensions (micro- and nano-scale particulates), the environment
(pollutants dispersion, particulates, microbes), and the petroleum industry [4, 5]. Short-
fiber composites are widely used in automobile bodies, business machines, and customer
parts. When a short-fiber reinforced polymer is molded, the mold filling flow changes the
orientation of the fibers. The fiber orientation in turn affects the physical properties of the
composite, including stiffness, strength, thermal expansion, and electrical conductivity. For
example, the composite is stiffer and stronger in the direction of greatest orientation, and
weaker and more compliant in the direction of least orientation. The most common types

of fibers used are glass, carbon, and aramid. Injection molding, extrusion and compression

molding are common manufacturing methods for polymer matrices.

1.2 Dilute suspension theory

The basic assumptions employed in almost all dilute suspension models are

(1) The volume fraction, (1) of fibers is so small that hydrodynamics interaction between

fibres or between a fiber and a flow boundary may be ignored.
(2) The particle size is small compared with the macroscopic characteristic length.
(3) The aspect ratio a p of the fiber is uniform.
(4) The suspending liquid is incompressible and Newtonian.
(5) The effects of inertia and external body force may be neglected ‘

(6) Nonslip boundary condition is applied on the interface of the particle and the fluids.

1.2.1 Hydrodynamics of a single particle

A single particle motion and hydrodynamic forces acting on the particle are fundamentally
important in the nature. Comprehensive information about the interaction between the
particle and the fluid in low—Reynolds-number flow is required for many practical systems
and industrial processes. Much is known at present about a single particle or two particles
in a creeping flow. Lamb’s classic treatise [6] on hydrodynamics contains much historical

and technical information on the development of solutions for creeping flows. Happel and

Brenner [4] developed the theoretical calculations of the Stokes resistance of a particle to
translational and rotational motions in an unbounded fluid. The motion of a rigid ellipsoid
in a uniform simple shear flow at a low Reynolds number is solved completely by Jeffery
[1] and verified accurately by the experiments of Trevelyan and Mason [7]. By using
Jeffery’s method [1], Bretherton [8] investigated theoretically on the orbit of a particle of
a more general shape in a non-uniform shear flow. The motion of non-neutrally buoyant
spheroidal particles is investigated with considering the effect of inertia by Broday and his
coworkers [9]. The resistance functions for two unequal spheres are derived by Jeffery [10]
and extended by Keh [11] to the slip problem on the interface of the particle and fluids.
Wetzel [12] set up an analytical model for the deformation of an ellipsoidal Newtonian
droplet, suspended in another Newtonian fluid with different viscosity and zero interfacial
tension. Transient wake flow patterns and dynamics forces acting on a rotating spherical
particle with non-uniform surface blowing are studied by Niazmand [13] for moderate
Reynolds numbers.

Hydrodynamics of a single particle include the relations between the hydrodynamic
force F, the torque T, and the stresslet 1: exerted by the fluid on the particle. Two kinds of
problems are classified in this area. One is from the viewpoint of mathematical boundary
value problems. The velocities of the particle and the surrounding flow field are fixed,
which supply for the suitable boundary conditions. Then calculate the force F, the torque
T, and the stresslet 1: . This kind of problem is called resistance problem defined by Brenner
[14,15]. The other problem is inverse to the first one, which is defined by Batchelor [16,17]
as the mobility problem. For this problem, the force F, and the torque T are given and the

relative motion of the particle through the fluid is to be determined.

1.2.2 Fundamental equations of creeping flows

pdlUl
p

Introduced in [14], the particle Reynolds number is defined as in the case of trans-

p d2 lwl

 

 

lating bodies or streaming flows, and , in the case 'of rotating bodies; U being the
translation velocity of the particle; d a characteristic particle dimension and to the angular
velocity. At small particle Reynolds number, the convective term pv - Vv in the Navier-
Stokes equation is very smaller in comparison with the viscous terms, szv. Neglecting

the influence of the convective terms in the Navier-Stokes equation, the dynamic and kine-

matic equations of motion of a viscous, incompressible fluid can be written as

6v
pE+Vp=pV2v (1.1)

and

V-v=0 (1.2)

where v, p, p, p, and t are respectively the local fluid velocity, the fluid density, the pressure,
the viscosity, and the time.

The local acceleration terms p% in the equation of motion is equal to zero for steady
problems. However, this term can also be ignored at a small particle Reynolds number
even for unsteady problems. The dynamic equations of motion of the fluid are therefore

simplified as

szv = Vp (1.3)

Consider a rigid particle immersed in an unbounded quiescent flow. The undisturbed am-

bient flow field is composed of the uniform streaming velocity U°° and the linear field

(constant velocity gradient), described by
v=U°°+Q°°xx+S°°-x. (1.4)

where x is the position vector of a point relative to the origin at 0, 9°° is the rotation of the
flow field, and S°° the rate-of-strain of the flow fluid. The motion of the particle induced
by the fluid has translational velocity U at a point 0, which is regarded as the origin of this
particle, and angular velocity u). If no-slip is applied at the interface of the particle and the

fluid, the instantaneous velocity of the fluid at the particle surface is
v(x)=U—U°°+(ro-fl°°)xx—S°°-x, xeSp (1.5)

in which 5 p is the surface of the particle.
The force, torque and stresslet exerted by the fluid on the particle about the origin of
the particle are F, T, and 1: respectively. The relations between the quantities F, T, and 1:

with U — U°°, (o — Q”, and S°° are to be determined.

1.2.3 Particle material tensor
The resistance tensor

When the specified quantities are the velocities of the particles and of the prescribed flow,
the linearity of the Stokes equations (1.3) permits the expression of the forces, the torques

and stresslets [5] in the form

F A E E uw—U
T =p B c E QW—m 06)
1: GHM S°°

The square matrix in the above equation is called the resistance matrix, in which A, B, and
C are second-order tensors, G and H are third-order tensors, and M a fourth-order tensor.

According to the reciprocal theorem of Lorentz [4], the resistance matrix is symmetric, i.e.,
Air = A1,, Cij = Cji’ Mijkl = Mkuj.
Bi} 2 Bjiv Gijk = Gkij, Hijk = Hkij (1-7)
The mobility tensor

When the particle forces and torques are prescribed in a known ambient flow, the so—called

mobility problem [5] satisfies the following relation

U°°—U a b ‘g’ %F
Qw—w = b c If fiT (L8)
Ill" g h m S°°

in which the square matrix is the mobility matrix. Similarly, the a, b, and c are second-order
tensor, g and h third-order tensor, and m the fourth-order tensor. The mobility matrix is

also symmetric as the consequence of the Lorentz reciprocal theorem:

aij = ajia Cij=Cjia mijkl=mklijs

~

bij = bji, gijk=§kija hijk=7fkij (1-9)

The resistance problem and mobility problem are inverse to each other physically. So
the resistance matrix defined by (1.6) and the mobility matrix defined by (1.8) are inverse

to each other. The transformation matrix between them can be found in [18].

1.2.4 Jeffery’s solution of the motion of an axisymmetric particle in a

linear field

In 1922, Jeffery [1] investigated the flow induced motion of an ellipsoid in an unbounded
flow field. Some assumptions are made in Jeffery’s model: ( i ) the particle is rigid, neu-
trally buoyant, and large enough to neglect Brownian motion, ( ii ) the ambient fluid is
Newtonian, ( iii ) the inertia forces of the particle and the fluid are negligible and the
motion of the fluid is governed by Stokes’ equations, ( iv ) the particle is immersed in a
homogeneous flow, which means the velocity gradient of the flow field is constant, ( v )
no-slip boundary conditions are applied on the interface of the particle and the fluid. Under
these assumptions, the time evolution of a spheroid orientation is expressed in the form
as [1,19]

0=9°°><P+4IS°°'P(P'P)-S°°:ppp] (1.10)

in which p is a unit vector representing the orientation of the particle. p can be characterized

with the angles 6 and (b as shown in Figure 1.1

sin 0 cos ¢
p — sin 6 sin (15

cos¢
(1.11)

length

 

xi is a function of of the particle aspect ratio, a p = , and is expressed as

diameter

2
p-l

2
ap+1

a

 

(1.12)

xi > 0 represents a prolate particle, xi = 0 a sphere, and xi < 0 an oblate.

 

 

 

Figure 1.1. Schematic of the coordinate system used to represent the orientation of a parti-

cle.

For a special case of a neutrally buoyant torque—free axisymmetric particle in the shear

field v°° = «y y ex shown in Figure 1.2. The ambient rotation and rate-of—strain can be given

by
010
Q°°-—l'e S°°—7- 10 o (113)
_ 2y 2’ —2 .
0 0 0

Substituting (1.11) and (1.13) into (1.10), the orientation of the particle is governed by the

coupled differential equations

 

 

_ ap—l y .

9 = — 2 —srn263in2¢ (1.14)
ap+1 4

(fl — 2y (agcosz¢+sin2¢) (1.15)
ap+1

Geometry of the particle

 

Shear flow: y

”=. \ ‘ 2 2
v ryex \ L+y_2+

2

2

 

 

=1

 

a
Q
GNIN

 

. C
Aspect ratio. a p = Z

 

 

Figure 1.2. An axisymmetric particle, i. g. a spheroid suspended in a simple shear flow. The

geometry of the particle is a spheroid and the aspect ratio of the particle is a p = g.

The above equations can be solved exactly and yield periodic trajectories known as the

 

Jefiery orbits [5],
Cap
tan 6 =
(a% 0052 q) + sin2 (12)”2
't
tan¢ = —aptan[—L——l] (1.16)
(1,, + a;

From these equations it can be seen that the motion of the particle is periodic. The period
T = 27r(ap + a;1)/y is proportion to the 7" and becomes longer with increasing particle
nonsphericity. The constant C is known as the orbital constant determined by the initial
angle of the particle. The exact shape of the orbit can be determined by the the particle
aspect ratio via the Bretherton constant B [8]. Jeffery orbits are shown in the Figure 1.3 for

different shapes of particles suspended in different shear flows. In these cases, the initial

 

 

 

 

 

 

 

 

 

 

 

Cos (15 C034)

1 l
0.5 0.5

/
t 10 20 t
—0.5 —0.5

-1 . -1 >

Cos ¢ (£05 ¢

1

0.5 » 0.5 * /

10 20 t 10 20 t
-0.5’ —0.5-
._1 ~

Figure 1.3. Jeffery orbits are shown for spheroidal particles at different aspect ratios a p and

 

 

different flow fields V,» = y yex. The motions of the spheroidal particle suspended in a
simple shear flow are periodic and the periods depend on the shape of the particle and the

constant «y of the flow field.

10

angles of the particle are (p = 0, and 6 = O. Suspended in the simple shear flow V,» = 7y ex,
the axisymmetric particle rotates along the z axis. It can be seen the period of the motion
of the particle with the aspect ratio of a p = 4 is longer than that of the particle with aspect
ratio of a p = 2 when the constants y are same for the two cases, and the period of a same

particle increases with decreasing the constant 7.

1.3 Non-Dilute Suspensions

Let v be the number of particles per volume. Solutions for rod-like suspensions, length L

and diameter b are distinguished by [20]

t 1
v < — dilute

L3
1 1 S . d'l
{ B <<V<< EL? emr— IUIC (117)

 

1
v = 0 (—) Concentrated
. dL2

Work has been done theoretically and experimentally on dilute suspensions [21—23]. With
the development of new materials, it is necessary to model fiber-fiber interactions and de-
scribe quantitatively the effect of the interaction to properties of the suspensions and the
final products.

The rheological features of nondilute fiber suspension have been observed experimen-
tally. A summery of the previous experimental studies on the nondilute suspension can be
found in [24]. Interesting results stated in those experiments include Blakeney [25], who
found that the effective viscosity of a steady-state suspension was approximately equal to

the viscosity of the solvent; Maschmeyer [26] showed a Power Law dependence of vis-

ll

cosity on shear rate; Miles’s experiment illustrated that the viscosity was dependent on
the volume fraction of fibers but independent on the shear rate and the length of fibers for
glass-fiber suspensions [27]; The shear-induced corrugation of free interfaces in concen-
trated suspensions was observed by Loimer [28], in which it was shown that the roughness
of the surface disturbances depended on the particle size, the particle concentration and the
fluid surface tension.

To model the interactions between suspended fibers, Batchelor [29] treated the semi-
dilute regime by using a cell model in the case of aligned particles in Newtonian fluids.
This model was applied to the non—Newtonian fluids for aligned fibers in extensional flows
[30]. In Batchelor’s work, the interaction among randomly oriented dispersed fibers was
modelled as a Newtonian drag, which was exerted on a fiber segment in motions relative to
the bulk suspension [24]. Similar ideas to model suspensions of non-Newtonian matrices

were used later in [31].

1.4 Objectives of this work

Hydrodynamics of a single ellipsoid suspended in an unbounded homogeneous flow was
investigated by Jeffery in 1922 and the flow induced motion of the ellipsoid was applied
later on to the analysis of dilute suspensions. Jeffery’s work dealt with the problem of in-
compressible Newtonian fluid with constant viscosity, no-slip on the interface of solid/fluid,
and linear shear flow. Recent advances in materials processing make it relevant to under-
stand the influence of a slip boundary condition, other flow fields, viscosities, inertia forces

on the drag force and rotary motion of a single particle as well as the dynamics of ensemble

12

particles.

The objectives of this work are:

(1)

(2)

Determine the influence of the slip boundary conditions, flow field, viscosity, and

inertia force on the drag force and rotary motion of a single particle suspended in an

unbounded flow field.

No-slip on the interface of the fluid/solid is an idealization in the transport process.
The slippage phenomena on the solid/fluid interface has been confirmed both exper-
imentally and theoretically. Unfortunately no equations describing the flow-induced
motion of an ellipsoid with slip are available so far. Some results are available for a
deformed sphere with the slip [32, 33]. These results are limited on the special case
of a fixed particle in the space. With the successful application of the motion of a
single particle to the analysis of dilute ensemble particles, it is an important step to
get the dynamic motion of an ellipsoid with slip to fully understand the microstruc-
ture of dilute and concentrated suspensions. A translating or rotating non-spherical
particle suspended in a slip flow is one of the main topics in this dissertation. Both
of the creeping flow and the Oseen flow are investigated to determine the influence
of physical parameters, e. g. geometry parameters of particles, fluid viscosities, and
the interfacial slip coeflicient, on the drag, torque, and rotary motion of a deformed

sphere in an unbounded homogeneous shear flow.

Develop a new closure model for Non-Axisymmetric particles required in the mo-

ment equation.

The shape of particles is another factor to the microstructure of suspension systems.

13

(3)

Previous studies of the orientation of particles have restricted attentions on the ax-
isymmetric particles. The current orientation of a nonaxisymmetric particle is spec-
ified by a rotation matrix from the reference configuration to its present orientation.
Closure approximation is needed in the governing equation of the second moment
orientation equation of nonaxisymmetric. This dissertation is also devoted to study

the microstructure of non-axisymmetric particles suspensions.

Predict the flow-induced alignments coupled with the migration of multiparticles.

Orientations of particles can affect the rheology properties of a suspension system,
and mechanical properties of composite materials. In order to treat the microstruc-
tural kinematics of a suspension system, lots of research has been done to predict the
orientation of ensemble particles statistically. For the suspension of axisymmetric
particles the orientation statistics can conveniently be given in terms of a unit vector
p along the axisymmetric axis and a second-order tensor a =< pp > is often used
to describe the orientation state of particles. in the previous research, uniform sus-
pension of particles is assumed. It is known that nonuniform suspension, shear force
and inertia force will cause particle migration. The influence of flow field parameters
on the orientation and fiber stress is studied for a specific two-dimensional problem
using finite element calculations. Therefore, to study particles orientations coupled

with the migration is another goal in this dissertation.

14

CHAPTER 2

MOTION OF A DEFORMED SPHERE

WITH SLIP

2.1 Introduction

With the evolution of micro- and nanoscale systems, there has been a recent wave of interest
in challenging the idea of a ”no-slip” boundary condition on liquid/gas flows. Slip has been
confirmed experimentally and theoretically by using sensitive force measurements [34, 35],
visual techniques [36] and molecular dynamics simulation data [37—41]. Situations for
which slip may occur can be encountered in solid particles suspended in rarefied gases, and
can also be encountered in liquid/solid systems such as polymer melts or water flowing in
thin hydrophobic capillaries [42—45]. Contrary to macroscopic flows, a small amount of
slip can strongly influence the transport phenomena and serious consequences may occur
for miocro- and nanoscale flows. On the design of nanoscale flow devices it is necessary

to understand the fundamental aspects of interfacial phenomena and particularly accurate

15

 

 

Vx(z) LIqUId

 

 

LI

Vslr'p X
a Solid

 

 

 

Figure 2.1. The slip length L, is defined for a simple shear flow in the presence of slip

boundary conditions at the interface of solids and liquids.

prediction of the fluid transport through tiny structures [46,47].

The degree of slip is normally characterized by an extrapolated slip length Ls defined as
the distance from the surface within the solid phase to where the flow velocity vanishes. The
definition of the slip length L 3 is explained in Figure 2.1. For most practical situations, with
simple fluids (composed of small molecules with a diameter d), small slip lengths d ~ L s
are generally expected. There are many factors that can affect the slip length including the
degree of hydrophobicity [48,49], the substrate topography and surface roughness [50—57],
the presence of interstitial lubricating layers [56, 58, 59], the polymer molecular weight
[60—62], and the applied shear rate [63—68]. Three levels of slip corresponding to the slip

length can be distinguished by no slip, finite slip, and perfect slip (infinite slip) as shown

16

in Figure 2.2. The relation between the slip length and the surface slip coefficient can be

 

 

 

 

 

 

 

 

 

 

 

A
2
Liquid

4

p

b

D

Vslip
.3
14 Solid
a. No Slip b. Finite Slip c. Perfect slip
L3 = 0 L3 = finite number ‘ L, = 00

Figure 2.2. Illustration of lengths Ls used to describe different slip cases. When L; = 0,
no slip occurs on the interface; when L, = finitenumber, finite slip occurs on the interface;

when L S = Infinity, perfect slip occurs on the interface.

obtained by [37]

L5 = — (2.1)

in which a is the viscosity of the fluid and ,8 is the slip coefficient on the interface.
Analytical solutions for Stokes flow past non-spherical particles with slip on the solid-

fluid interface are limited to flows around spheroids fixed in space [32,33] by using a stream

17

function formalism. In this approach, the slip boundary condition can only be satisfied on
the surface of a sphere particle even though solutions for flow around a slightly deformed
sphere are given in those previous work. Equations describing the motion of spheroidal
particles in creeping flows with slip are not available, unlike the case of ellipsoidal particles
with no-slip surfaces. The latter equations are available in the classical work of Jeffery
[1]. In his paper, the behavior of an ellipsoid suspended in a uniform shear flow field is
analyzed based on Stokes’ equations of motions. Jeffery’s work is extended by Bretherton
[8] to investigate the orbit of a particle of a more general shape in a non-uniform shear
flow. A series of research papers on the intrinsic hydrodynamic resistance of particles of
arbitrary shape based on the application of singular perturbation techniques are studied by
Brenner [2, 69—71] presented. Jeffery’s equation has been applied to the analysis of dilute
suspension problems [19, 21, 39, 72, 73]. Corresponding equations for the induced motion

of an ellipsoid with slip boundary conditions are however not available.

2.2 Basic formulations

The slow motion of a slightly deformed sphere in an unbounded Newtonian and incom-
pressible flow is considered. Apart from the disturbance produced in the immediate neigh-
borhood of the particle, the motion of the fluid is assumed to be quasi-steady, and variable
in space on a scale which is large compared with the dimensions of the particle. The fluid

is allowed to slip over the surface of the particle. At small translational and/or angular

roIUI rglml

 

and

 

Reynolds numbers, , respectively, the fluid motion is governed by Stokes

18

equations,

,uVZV = Vp, v - v = o (2.2)

Representation of the general solution for the Stokes equations in terms of solid spherical

harmonics and surface spherical harmonics is given by Lamb as [6]

 

 

_ 0° (n + 3) _ n
v 7 n :2.» [V X (”f") + W” + 2(n + 1)(2n + 3)/1 VP" rm + 1)(2n + 3)ppn (2'3)
P = Z Pn (2.4)

Viz-'00

where X”, 43,1, and Pn are solid spherical harmonics of degree n and can be determined
by suitable boundary conditions. Lamb’s general solution to the Stokes equations assumes
the velocity field vanishes as r —> 00. In the event that the velocity field is required to be
prescribed at v00, the specified velocity V,» and corresponding pressure p,>0 are added to the
right side of (2.3). Naturally vDO must itself satisfy Stokes equations.

The surface of a deformed sphere is assumed to be described by an equation of the form

r = r0 (1 + ef(6,¢)) (2.5)

in which (r, 6', (b) are spherical coordinates and the origin is located at the center of an
undeformed sphere with a radios of r = r0; | sl is a small dimensionless parameter, and
f (6, p) is an arbitrary function which can be approximated by a series of surface spherical

harmonics, fk(6, ¢)~ Hence, the surface of the deformed sphere can be written as
r = r0 [1 + 82 mam) (2.6)
k=0
The slip boundary condition at the solid-fluid interface is modeled as [6, 11]
u = v — EU" (2.7)

19

where u and v are the velocities of the particle and the fluid on the interface respectively, ,8
the slip coefficient. on is the tangential traction on the surface of the particle. The traction

on the particle can be decomposed by two parts

0r :

Q

l'
r

‘NI :l

['1'
-or+(I-—,)-or
V ELF—J
O'rr 0'rt (2-8)

in which or, is the normal traction and on the tangential traction. The motion of the particle

on the surface can be expressed as
u=U+toxr onSd (2.9)

in which U is the translation at the center of the particle, to is the angular velocity of the
particle around its axes. Sd represents the surface of the deformed sphere. The boundary

condition at infinity is described by
v = voo(x) as r —9 oo (2.10)

where v,>0 is an arbitrary flow field.
For the problem of slip flows pasting a deformed sphere, solutions are sought by as-

suming that the velocity and pressure can be expanded in powers of s in the form

00

v = Z Siva) (2.11)
i=0

p = sip“) (2.12)
i=0

0,, = Zsiog) (2.13)
i=0

20

00

V,» = 2 avg? (2.14)
i=0
m u .
u = s'um . (2.15)
i=0

Due to the linearity of the Stokes equations both of the individual perturbations v“) and p“)

still satisfy Stokes equations.
pvzv“) = me, V - v“) = o (2.16)

Substituting (2.11) into the boundary conditions (2.7) and (2.10), and equating terms in-

volving the orders of s, it can result:

. . 1 .
u") — v00 = v“) — Bug) on 5,, (2.17)
v“) = v52 at r z oo (2.18)

(i)

Approximating v“) and or,

using Taylor series expansion about r = r0, boundary condi-
tions on the interface can be written as:

00 i i 00 i i 1 i
2qu = 28 {[14 i), z .0 — gins). ._. .1

i=0 i=0

i -_ i-D
1 . . . arrive D 1 6(Do(
+ E:_81,ij(g ¢) [(__) _ _[__Cf_ (2,19)
. O 9
i=1]! 6r") r=ro '8 ‘9’”) r=ro

atr=r0

It is difficult to obtain solid spherical harmonics X", pn, and ¢n by directly substituting the
boundary conditions (2.17-2.18)into (2.3). Three identities are instead introduced to the

slip boundary conditions for the each power of s [2]:

:_ (O)_ (0) :E_ (0)_l (0) 220
r [u 1’00],er r v .30" r=r0 (O)

21

—rV - [n(o) — v9], = r0 = —rv-

 

1
van - 4153)] (2.21)
B r : r0

r-vx[u(0)—vf,g)] r-vx

 

(0)‘ 1 (0)]
_. v — —o (2.22)
r - r0 :3 n r = r0

and

 

(0)
l' 1) 612(0) 1 0'" l‘ l 1 (1)
- . (u(1))— vfx, (r, 9, ¢) [( )— — {—H} = — - [v( i— —o (2.23)
r { 6r ,8 are) r=ro r ,8" r=r0
(0) (0)
(1)_ (1) 6__v _1 0n _
r = r0

(0) (0)
<1) (1) 6V -1 Se.
r V X {0' ) v00 (r 6 m i(_ 6r _) '8 (6r(fl]i}r

 

v(1)_ 1 (1) (2. 24)
30, r— - r0

r-Vx[v(l)— [130" (I) (2. 25)

zm r=m

2.3 Relation between pn, (bn, and X" with X", Y,, and 2,,

Substituting (2.3) into right hands of the three scalar identities (2.20-2.22), it yields,

:- (i)__ 10(1') (r7 0)" (0+ "(in)" (i) 226
r [v ,B 0" r=r0 =n__oo[2—_(2n+3),u p" +r0 r " (° )

. 1 - n(n + 1)r0( r_())" (i) n(n - 1)( r O (D
_ V (l) (0 —— 2
r [V on] — r0 :i2(—2__n+3)l1(r)pn + r0 (2. 7)
1 2n(1-n2)iu(r_o)" (i) n 2("+2))('o)" (x)
+1:an '02 r ¢n 2n + 3 r p"

22

. (0-1 (1) (i)_ 11109-1)" ro)" x“)
r VX[v fio”ir=x0= n=-oo[n(n+l)( —"")x),, 48 r0 (—) (2.28)

When the velocity field is prescribed at the surface of the particle, the left hand sides of

(2.20)—(2.25) can be represented by surface spherical harmonics

ll
M8
:8

E. (0) _ (0)

r [u v00 ]r = m (2.29)

n=l

00
—r V {11(0) - V2900, 69¢)lr = r0 = 2 n0) (2’30)

n=1

0 0° 0

r - v x [n(O) — 8,380, a, ¢)]r = r0 = Z z}, ) (2.31)

and

 

 

' (0) {0(0M] 00
E, (1) _ (1) 6" _l “r: = (1)
r {(u ) voo(r.6.¢) ( ax) BUM» } 2X, (2.32)
L 3 r=r0
lav<0> 1‘0“») oo 1
_ )
( 6r )73 arm ‘Zy'i (2'33)
k l. [‘er

(o) lio(0)il oo
. (1) _ (1) 3V __ rt = 1
1' V><{(u ) voo(r,9.¢)L( ) warm} _ 2,25.
- r—ro

 

 

—r v - {(u‘”) — vél’m 6. c»)

v

 

 

(2.34)

 

 

 

 

By replacing n by —(n + 1) in those terms of (2.26)—(2.28) for which n is negative, the sum
can be made to extend form n = l to 00 rather than —00 to 00. Hence the following three

relations are obtained

 

0° (1') _ °° ' "to (r; 0)" (i) _(ro)" (1')]
X _ __ — (2.35)
"A: " §L2(2n+3)ju pn+ r0 r n
+ EM ' ‘0'“) pi() "+1 r ’0'”) (i)
n—1.2('2"+1)# ro AW” r0 r0 p ("m

23

(i)_ n(n + 1)r0(__ (1') n(n — 1) Q n (i)
trim Y" Z 1i2——_(2" + 3).U(_ )npn + r0 ( r ) ¢n ] (2‘36)

4;; — We) ——(-’£’-)"p%°]

+2 —n(n + ___1_)r0(_ ("+1) (1') + (n + l)(n + 2) _r_ —("+1)¢(,~)
2(2n — 1))u( ” -(n+1) r0 -(n+1)

 

 

 

 

 

r0
‘2 2(n+ l)(n2 +2n),u r(”+”¢(i) _ (n+1)2(n-1) L "M” (i)
r3 70’ —(n+1) 2n _ 1 ,0 p—(n+l)

°° - °° l (n2 — 1)n r ’1
225:" = Z ["0” ”l 70),!“ (0 51—70—09) *9] (237)
n=1 n= 1

oo ’ —(n+l) -(n+l)

r (i) 1 ,un(n + l)(n + 2) r (i)
+ n(n + 1)(—) )(_ — — — X
"21) t ,0 (n+1) fl ,0 r0 -(n+l)

For an exterior problems, conditions that the fluid should be at rest at infinity require
that

pg)- - S): XSL — 0 for n 2 0, (2.38)

so only the negative harmonic functions survive in (2.35)-(2.37). Due to the orthogonality
of surface harmonics of different orders, the left— and right—sides of the resulting equations
are equated term-by-term under the summation signs. The resulting equations may then be

solved simultaneously for the three harmonic functions. When n 2 0, it is found that

(1') (2n — 1)73a(2r0,8Xg) + ronfixg) + rQBYS) + 4n].rX,(,i) + 2n2pX,(,i))

 

 

 

= (2.39
p ("+9 (n+1)r("+1)(r0)8+,u+2n/.r) ‘ )
$8) _ r3"*2’p<2rorXfP war}? —2qu." +2n2uX§i5 (240)
‘0’“) 2(n + 1)r("+1)(r0l3 + p + 2np) '
(n+2) Z(i)
. r ,8
ngmq) = 0 n (2.41)

n(n + l)r("+1)(r08 + 11 + 2n/r)

24

Substituting these relations into Lamb’s solution (2.3), the velocity and pressure field can

be written in the form

(1‘) _ (i) (i) _ (n - 2) '(i) ('1 +1) (0
Z[VX(YX_(H+1)) +V¢—(n+l) 2,10"- 1)#’2VP—(n+1)+rn(2n_ ”pp-(n+1)
(2.42)
and
00
' _ (0
p0) _ Z p-01“) (2.43)
n=1

The hydrodynamic force and torque exerted by the fluid on the deformed sphere can be

written as [2]

F = 2 air“) = Fm) + 3F“) + 0(3) (2.44)
i=0
T - 2 8‘1““) = Tm) + 8T“) + 0(3) (2.45)
i=0
in which
F“)- — —47rV(r3p 9,) (2.46)
and
T(')— — —8an(r3X(_‘)2) (2.47)

2.4 Stokes’ resistance of a spheroid

2.4.1 Functions to describe the shape of a deformed sphere

The flow past a spheroid is considered here and the hydrodynamic force and torque exerted

by the fluid on the surface of the spheroid are determined. The shape of the spheroid can

25

be described by

 

2 2 2 P 1 t h 'd h < 0
x :y + z :1 roaesp eror wens (2.48)
r0 r30 — (5)2 Oblate spheroid when a > 0

To the first order in the deformation parameter s, (2.48) can be written in a polar form as

in [74]
1 2 2
r = r0 1 — e §p0(cos 6) + §p2(cos 6) + 0(8 ) (2.49)
Comparing (2.49) with (2.5) results in
1 2
f(6) = — {3p0(cos 6) + §p2(cos 6)} (2.50)

where cos 6 = 5. p0(cos 6) and p2(cos 6) are Legendre functions given by
r

1
p0(cos6) = 1, p2(cos6) = -2-(3c0526 — 1) (2.51)

2.4.2 Uniform streaming flow past a stationary deformed sphere

The undisturbed uniform streaming flows past a stationary spheroid with velocity v.>0 =
U e x + Vey + Wez is considered in this section. The expansion of u and v00 can be expressed

as

u“) = o i 2 0 (2.52)

vi?) = Uex + Vey + WeZ v32 = 0 i z 1 (2.53)

26

The hydrodynamic force and torque acting on the surface of the spheroid as a whole are

 

 

 

 

therefore
F : 6me(y + 2 _ 8(272 +117+ 6))e
y + 3 5(y + 3)2
2 2
+ 67rr0/1V(y + - 8( 72 + 117 + 6)) , (2.54)
7 + 3 5(7 + 3)2
2
+ 67rr0pW(y+ 2 _ 8(7 + 87+ 18)) z
'y + 3 5(7 + 3)2
and
T = 0 (2.55)

in which y = rig is a dimensionless parameter.
For a streaming flow parallel to the spheroid axis with velocity v = vooeZ , the hydrody-
namic force and torque exerted on the spheroid by the fluid can be obtained as

67rr0voop (y2(—5 + e) + y(-25 + 8.9)): + 6(—5 + 38)p2)
— 5(y + 3)2 '

 

ez (2.56)

and

T = 0 (2.57)

For the limiting cases of the perfect slip and non-slip boundary conditions, the hydrody-

namic forces are found to be

i

471' ro v00 ,u(1 — %)ez when ,6 —+ 0
F = t (2.58)

 

\67rr0voop(1-§)ez when fl—>oo

For an undeformed sphere, for which a = 0, the hydrodynamic force with slip is

 

F = 67rr0voo/r (1 — 7’: 3)eZ (2.59)

27

which is exactly same as Basset’s solution [75].

The comparable results for streaming flow perpendicular to the spheroid axis, v00 =

 

Vooex, ar€ _
2 11
67rr0voop [72(1 - 755) + y(5 - —SE) + 61120—2”
(7 + 3)2
and
T = 0 (2.61)

For special cases of perfect slip and non-slip boundary conditions, the forces can be ob-

tained as e
i
47rr0vooa (1 — 3)ex when ,8 —> 0

F = i (2.62)

 

l 67rr0voop(l — %) eJr when ,3 -—> 00

Equations (2.58) and (2.62) for the hydrodynamic forces of non-slip boundary condition

are exactly same as Brenner’s expression [2].

2.5 Uniform streaming flow past a rotating deformed

sphere

Let the spheroid rotate in a uniform streaming flow about its axis by (t), which can be
decomposed as

to = wxex + myey + wzez (2.63)

in which w,- is the angular velocity component along the i direction. Boundary conditions

appropriate to the rotation of a spheroid are given by
r
usz-r on Sd (2.64)
r

28

On the surface of the particle, each of the perturbation of the spheroid velocity can be

written as
um) = m x Fr-ro . (2.65)
u“) = m x Er0f(6, ¢) (2.66)
u"? = 0 i2 2 (2.67)

The undisturbed motion of the fluid in the neighborhood of the spheroid is given as

V,» = Uex + Vey + Wez (2.68)

The hydrodynamic force and torque can be obtained as

 

 

1:5
ll

y+2 3(272+11y+6)
2+3_ 56+3V ix
y+2 8(2y2+11y+6)
7+3_ 5(y+3)2 )y
y+2 8(2y2+8‘y+18) -
y+3— My+oz )2

67rr0/.rU (

(2.69)

+

 

67rr0,uV (

+

 

 

67rr0,uW (

and

a
ll

 

 

4 _ 8 488(r0j8 + 4p))
”'03”“”‘( mg + 3p + 5(r06 + 3p)2 8"

4 _ 8 483(r0/3 + 4p)
"’Wwyi me + 3n + 5668+ 311V )°’

4 _ 8 248(7‘03 + 411))
mofiflwzi roB + 311 + 5(r66 + 3102 ez

+

(2.70)

 

 

+

 

 

At the first order of s, it can be seen that the translational motion of the fluid relative to
the spheroid determines the hydrodynamic force, while the rotation motion of the particle

determines the torque on the spheroid.

29

For a spheroid rotating about the symmetry axis in a quiescent flow, to = wzez, the

hydrodynamic force and torque exerted by the fluid on the spheroid can be written as

F = 0 ' (2.71)

 

T = 8nrgap[y(-5 + 38) + 3(—5 + 4e),u]wzez (2.72)
xy+mz

of which the limiting cases are

{0 when ,B=0

T = 8
gm'3(—5 + 3e),uwzez when )8 —> oo

(2.73)

Analogous results for rotation about the equatorial diameter to = wxex in a quiescent flow

 

are
F = 0 (2.74)
and
87rr4 rim—5 + 68) + 3(—5 + 88)].l]w
T = OB xex (2.75)
5(rOB + 3102
For the limiting cases of the perfect slip and non-slip, the torques are
0 when 6 = 0
T = (2.76)
gnr3(—5 + 68)pwxex when )8 ——> 00

For the non-slip case, the torque on the deformed sphere is same as Brenner’s solution [2].

For a sphere rotating in a quiescent flow with slip the torque is given by

Sargflwxex
T = —— when 8 = 0 (2.77)
y + 3

2.6 Flow-induced motion of a spheroidal particle with slip

A rigid deformed sphere is considered in this problem to be suspended in a homogeneous

flow. The problem is expressed in two coordinate systems. The first system rotates with

30

 

 

 

Figure 2.3. Illustration of three Euler angles ¢, 6, and i]! used to describe the coordinate
systems of a particle. x’, y’ and z’ are the reference coordinate system and x, y and z are the

rotating coordinated system.

the particle and is denoted by x, y and z. The second coordinate system is fixed in space
and denoted by x’, y’ and z’. The position of the rotating coordinate system with respect to
the fixed system can be described in terms of three Euler angles shown in Figure 2.3.The
angle 6 is simply the angle between the z axes of both coordinate systems. The angle ()5 is
the angle between the x axis of the reference coordinate system and the projection of 2 into
the x’ , y’ plane. Finally, ip is the angle between the y axis and the line of nodes. These two

systems are connected through the transformation

x cos¢cos¢—cos6sin¢sin¢ coswsin¢+cos6cos¢sin¢l sin6sinil/ x’

y = —sinzl/cos¢-cos6sin¢cos¢ —sin¢/sin¢+cos6cos¢cosr/I cosr/zsin6 y’

2 sin 6 sin (b — sin 6 cos (12 cos 6 z’
(2.78)

31

Rotations of the particle around the body axes may be described by

to = (cox coy LUZ)

( sin i/l sin ¢6 + cos W} cos t/J sin (06 — sin ([16) cos (156 + (b ) (2.79)

A homogeneous flow field far from the particle can be defined in the rotating coordinate

system as
a h g 0 —{,’ n x
Voo = h b f + g 0 —.f y (2.80)
8 f c -n 6 0
a h g 0 —{ I]
in whichS = h b f is the distortion of the fluid andW = g 0 —§ is the
g f c -n E 0

rotation of the fluid. The parameters are constant in space, but can be a function of the

time, since the particle will rotate under the influence of the flow field. Let

G=S+W . (2.81)

The velocity field at infinity of the fluid can be written into a tensor form as

V,» = G - r (2.82)
and

vi? = G - Ero (2.83)

v2) = G . £r0f(6) (2.84)

v52 = 0 when i 2 2 (2.85)

The hydrodynamic force and torque applying on the surface of the spheroid are

F = 0 (2.86)

and

T = Txex + Tyey + Tzez (2.87)

in which

_ may {f8(3 + y)(32 + 38y + 5%) — m + 5)[-5(7 + 3) + 68(7 + 4)](6 — on}

X'—

 

 

 

5(7 + 3)2(7 + 5)
(2.88)
T _ _ 87rrgp {g8(3 + 'y)(32 + 38y + 5%) — y(y + 5)[—5(y + 3) + 68(y + 4)](77 — op}
y ‘ so+3>2<~y+5>
(2.89)
Z 2 _ 87rr3p[y(—5 + 38) + 3(—5 + 4a)](r — (oz) (290)

fly+®2

If a particle is subjected to no external forces except those exerted by the fluid upon its
surface, then the resultant torque on the particle will be vanished at each instant. Let the
three components of the resultant torque equal to zeros. The angular velocities of the

particle can be solved to be
(fix: _f/l+§, wy=gd+7la (02:4,, (291)

where

= (3 + y)(32 + 38y + 5y2).«.~
7(5 + 7)[-5(3 + 7) + 68(4 + 7)]

 

(2.92)
and the dimensionless parameter 9/ = r%B- is introduced. If y —> 00, the no-slip boundary

condition is satisfied on the interface, and (2.92) can be simplified to

58
-5 + 68

 

z -8 _ 282 + 0693) (2.93)

33

Actually xi can be decomposed into two parts as

,1 _ 5.9 + s(—480 — 3557 — 65y2) + 82(576 + 276y + 48y2)
’ —5 + 68 7(5 + y)(—5 + 6s)[—5(3 + y) + 68(4 + 7)]

 

 

(2.94)

In equation(2.94), the first term of xi is the no-slip part and the second term is the additional
part due to slip on the interface.

The evolution of the orientation of a spheroidal particle can be described with [76]

P=—W'P+4(S'p-S:PPP) (2.95)

in which p is the orientation of the particle which is a unit vector along the long axis. This
equation is identical to Jeffery’s equation for the orientation of the particle with exception
of the definition of parameter xi. With consideration of the slip on the particle surface, it

is convenient to group xi into parameters related to the geometry of the particle, the slip
2
a

 

coefficient, and the viscosity. In Jeffery’s solution [1], xi J = 1 is only related to the

a
geometry of the spheroid, i.e. the aspect ratio the ellipsoid app(1ength/diameter), and the
range xi J is [—l, 1]. xi J = 1 represents for long fiber with infinite aspect ratio, xi J = 0
indicates for sphere, whereas xi J = —1 corresponds to a disk. The relation of xi with the
dimensionless variable y when 8 = 10.1 and e = $0.2 are shown in the Figure 2.4. It can
be seen that xi is always positive when a < 0 and negative when 8 > 0. For the case of
e = —0.1 only one positive root y = 0.1239 exist to let xi = 1 shown by the dash line in
the Figure 2.4. Jeffery [1] found the equations describing the motion of an ellipsoid in an

unbounded fluid with nonslip boundary conditions. When simplified to a slightly deformed

sphere, it is found that

can = -f41 +6, wa = gxij + 27, sz = 4', (2.96)

34

     
 

 

5 .|\ —_' €=—O.2
4 _| (0.1239,1) ___. £=-O.I
3 . (0.2407,1) 8:0.1
2

 

 

(0.1786,-l) (0.5746,- 1)

 

    

Figure 2.4. Influence of the slip coefficient 6 and the particle aspect ratio a on the
fluid/particle coupling coefficient xi. xi is always positive when a < 0 and negative when

£>0.

where

’1 8(8 — 2)
’ .92 — 28 + 2
= —s — $32 + 0(83) (2.97)

By the present approach, for the limiting case of nonslip boundary condition, the error of
the present approach is 0(82), and this is also observed after comparing (2.91) and (2.93)
with Jeffery’s solution [1], i.e. (2.96) and (2.97).

The case of a simple shear flow imposed far from the particle is considered next and
described with vfx, = ( 0 0 Ky’ ) where K is a constant. The velocity field at infinity is

shown in Figure 2.5. From (2.78) the distortion and rotation of the undisturbed fluid in the

35

 

Figure 2.5. Illustration of a simple shear flow surrounding a deformed sphere.

rotating coordinate system are found to be

0 0 0
1
S = 0 K cos a sin a -2-K cos 2a (2.98)
0 éKcos 2a —Kcosasina
0 0 0
K
W = 0 0 3 (2.99)
K
0 —3 0
The motion of spheroid can be described with
cox = (I, Lay = 0, 612 = 0 (2.100)
Using (2.91) the rotation of the particle is
, K
(1),; =0: —-2—(xicos2ar+ 1) (2.101)

36

a ---- x=—2.0
--- A=—0.3
21:60
x=o.5
A=3.0
A=1.0

 

 

 

 

 

 

 

Figure 2.6. The influence of the fluid/particle coupling coefficient xi on the rotation of a
spheroid in the flow/shear plane is shown by plotting the cosine of the rotation angle a
of a particle for various values of xi. When Ixil < 1, the induced motion of the particle is

periodic; when |in 2 1, the induced motion of the particle is steady.

Integrating (2.101) over the time domain, the angle a can be expressed in terms of the

dimensionless time 1' as

 

 

2 _ (IT:—
a = arctanl :_ 1 1 tanh (- ’12 17]] (2.102)

in which 1 = Kt. The evolution of angle a is shown in Figure 2.6. It can be seen that when
|in < 1, the motions of the particle are periodic; when MI 2 1, the particle rotates to a fixed
angle and then reaches a steady state. The period increases with the absolute value of xi.
In Figure 2.4, the horizontal lines xi = 1:1 separate the steady motion and periodic motion
of the particle. Between these two lines, the motion of the particle induced by the simple
shear flow is periodic. For the cases of |xi| 2 l, the steady state of the particle means that

the orientation of the particle does not change with time. Taking the time derivative of the

37

orientation 0 equal to zero, the orientation of the particle reaching a steady state can be

obtained as

6,, = i (2.103)

xi-l
— — xi>1
arccos 2’1 _

 

L
From Figure 2.6 it can be seen that when xi < —1, the angle ass is positive, and when xi > 1,

the angle ass is negative.

The phenomenon of a periodic motion or a steady state of a force-free particle inside
a simple shear flow can also be explained physically. Eq. (2.95) stems from a balance
of angular momentum. The underlying assumption in the Jeffery analysis is that there
is no external torque acting directly on the particle. Therefore all the intrinsic couples
(torques) between the particle and the fluid must balance and the angular momentum of
. the particle is constant. As a consequence, the angular velocity must satisfy (2.95), which
requires the orientation of the particle to change continuously if |in 3< 1 (i.e. , periodic
behavior). Eq. (2.95) shows that the angular velocity is caused by two processes: 1) a
coupling of the particle with the vorticity of the imposed flow field; and, 2) a coupling of
the particle with the strain rate of the imposed flow field. A vector produced by a pure
rotation of the orientation vector represents the first effect: —W - p. This vector is always
orthogonal to the orientation vector. If the external flow field has vorticity and if xi = 0
(sphere), then the sphere must rotate continuously to balance the torque induced by the
external hydrodynamic field. For a neutrally buoyant system (i.e., density of the fluid and
particle are equal), the sphere translates with the local velocity and a drag force is exerted

on the particle, which is countered by the external force needed to sustain the flow field

38

Table 2.1. Predicted periods of the particle motion induced by a simple shear flow. The
periods of the motion of the particle increase with the increasing of the deformation of the

particle and decrease with the increasing of the parameter 7.

 

period 7 = 0.2 y = 1.0 y -—> oo Jeffery
s = 0.1 30.2 14.1 13.3 13.3
s = 0.2 steady 17.4 13.2 13.1
e = 0.3 steady steady 14.5 12.7
e = —0.1 16.5 13.2 12.5 12.5
e = —0.2 steady 14.3 13.1 13.2
e = -0.3 steady 15.6 13.4 13.5

 

 

 

 

 

 

 

 

 

 

and the particle velocity. If xi 9t 0, then the particle can also couple with the strain rate
of the external flow to produce a torque. (2.95) shows that the resulting angular velocity
is proportional to the coupling coeflicient xi and a vector produced by a rotation-stretch-
projection operation on the instantaneous orientation vector: xi(S - p — S : ppp). This vector
is also orthogonal to the orientation vector. If lxil < 1 , then the torque produced by the
strain rate is too weak to balance the torque produced by the vorticity field. Consequently,
in order to satisfy the torque balance, the particle must rotate (Jeffery orbits). If |in > 1,
the torque on the particle due to the coupling with the strain rate is now large enough to
balance the torque due the vorticity field in the absence of tumbling! Most significantly,
the particle is not aligned with the flow field. This surprising and interesting result is due
to the slip phenomena.

To illustrate with a simple problem, the parameters affecting the motion are selected

to be r0 = 0.01, ,u = 0.02, and K = 0.5. e and y are changed to study the influence of

39

  
  
 
  

"" 7:0.2
-—'y=l.0

._.... 7:00

 

Jeffery

 

 

 

 

 

 

Figure 2.7. The figure shows the influence of the fluid/particle coupling coefficient xi and
the particle aspect ratio a on the temporal response of a spheroid to a steady simple shear

flow at positive values of as.

40

“' ' y=0.2
—- y=1.0
..._... 7:00

— Jeffery

 

 

 

 

 

 

y=0.2
- x" ' - '\ —7=1.o
-_.. 7:00

_ Jeffery

 

 

 

Figure 2.8. The above figure shows the influence of the fluid/particle coupling coefficient
and the particle aspect ratio on the temporal response of a spheroid to a steady simple shear

flow at negative values of as. 41

r=0.2

          

 

1 b

P f r‘ll‘vI‘Il."ll

 

fiffi‘

DDDD\
bL‘V
.|\\.
10"
(5"

IIIIIIIIIII‘ITTTIIJ
TOI““‘DII”‘I”‘1‘ITT

v--\D\.‘\ ‘1””"“‘.

 

"4
r77a§
v97).
7"!)‘
r""’I
'U"”’a'.’

10"“1'1'1'101'

'i'

I'J'I' l'l'l' I'JI'I'II'I'i'I"|'I‘I‘

ff fl'fl'

 

 

'

 

 

 

 

0015
0.0
00

 

 

4 1
‘11““‘ II A ‘

III-1““‘IOII

TTDI“D\‘.II

riIITTIIOIOJ
‘1"1r’lr'PL
11”""0.
.Illl‘dl

V CII‘L

' ’4‘;

.K '44“

\‘d

 

1" " U I ’1'
i'l'l'l'l"l'l'rl'
T'JI'I'I'I'I'I'l'

iii

 

I'I'l‘l'l'i'I'l‘fn‘4

 

.‘il‘ ‘l‘ ‘1‘

 

 

 

<15 0 5
0. 0. m.
0 0

O

-0.015—0.01—0.005 0 0.005 0.01 0.015

—0.015—0.0l—0.005 0 0.005 0.01 0.015

r=0.6

   

 

‘1

 

4 4
FI“OI.\II ‘ ‘

.01T““D\DIIT”"‘I‘I”
—

ll'll‘IJ

 

NDPTDRRRTQFIHIIFF“.

Ila.

 

 

"‘
L 1
1st;
-\!;
K“.

T"""“
*‘I'U'il'r'l'l'

 

\“‘
\\“‘

fr‘.‘“t‘".'

 

 

iffffl'

 

. piflfibiPlPliD .b Iiiiiilvlitil} lulPIlliPibLif hilt

5 m
m. 0.
n_U _

 

fif’ff’l",

 

‘lfi‘l‘l‘

I“““0i‘lI‘l’T-OIOI‘I‘J‘I
T‘I““‘TIT""‘7TTT

.0“.“\.‘”’I””“

ll'lll“

  

vvlllx

rillA
«via.

"‘

 

T'P'VII’
1'.".'."1'i'

7"“r'l'r'l'l'

 

.I'I'I'I'l'l'l'l‘i‘

l\‘j
.fi“;
\“_
\“'

\\“‘

 

It‘t‘““'-.

T'1‘\\‘\~L‘u‘t'

 

 

 

5 O
m. 0.
n_u _

—0.015—0.01-0.005 0 0.005 0.01 0.015

—0.015—0.01-0.005 0 0.005 0.01 0.015

    

t=0.8

w m w
0 .
0. 0 0.

 

4

 

iii 4 ‘i

Iii“ ‘IIIII‘lf‘lI 0101011

TDIDI‘“‘DIIT”‘I‘I’OITT

.‘f‘b-“~!‘.’””“-

Illlllll
'IIIIL
will.
Afl‘.
vll

 

'U'Fix’n'n'

l'n'il'l'r'l'l'

I. I'I'II'I' I'I'J

,.I\‘l
.\‘1
8“.
\“‘
\\“‘
Iii
T~\‘““"

I'\‘\‘\‘t‘u“.‘1'

I'I'I'I'l'l'i‘l‘l‘

 

 

 

.i‘l‘l‘l‘l‘l‘fi

 

 

 

II‘lIiIiiii
I‘l““‘.|‘l
01.1““‘0101
.b-‘b|.\.|.bui
bb\b\ .i,ii1
,1 n s »
u- \ \
1 D ‘4
F b ‘)

{fff’f‘l‘lf
‘I‘I‘I‘I‘r‘u‘v’f
’I‘.’""‘.
I’lll“J

' ' (ll‘L

   

will.
..A(444
'0‘

 

0"}

V'I'Ilrl

 

.r'"".'.'r'r'

.
W‘l'r'r'l'l'l'l'l'

.
T...l"l.'i'rl'il‘t

  

.I'I'l'l' 1'1'1'1‘1‘

(\‘j
.“l
\\‘.
\“‘
\‘\“
Ir‘x“\“"

I'i‘i‘t‘x‘l‘.“r'

 

 

 

O

5 10..
m. 0.
0. _

—0.015—0.0l-0.005 0 0.005 0.01 0.015

—0.015—0.01—0.005 O 0.(X)5 0.01 0.015

Figure 2.9. These graphs show the velocity field around of a spheroid with steady motion.

The particle reaches a special orientation after some time and then keep this orientation

inside the flow field.

42

0.015
0.01 i

0.005 i

    

—0.005 -0.005[ ‘3

-0.01

 
    
    

—o.01 i _

 

 

 

—0.015—0.01—0.005 0 0.005 0.01 0.015 —0.015—0.01—0.005 0 0.(X)5 0.01 0.015

1:04 1:06

0.015 0.015 _

    
  

0.01 .

0.005 E . 0.005_

 

 

0 0‘
-0005 i —0.005
—0.01 ..._ —o.01
—o.0150.01—0.005 0 0.005 0.01 0.015 —o.015-o.01—0.005 0 0.005 0.01 0.015
:08 _
0.015 I . 0.015 T—l'o

0.01 if

0.... 7‘)
o

—0.005

—0.01

 

 

- i

—0.015—0.0l—0.005 0 0.005 0.01 0.015 ~0.015—0.01—0.005 0 01115 0.01 0.015

Figure 2.10. These graphs show the pressure field around of a spheroid with steady motion.

43

ii...x.i. ii.

 

 

I‘I‘IIIIII‘lI‘i‘iI‘ii

tv\\\\\\\\
i\\\‘\\\\|\

viiii'ifiVi'irViiivAAriiiii

 

 

 

 

 

 

 

—0.015—0.01—0.005 0 0.005 0.01 0.015

 

 

 

-0.015—0.0L—0.005 0 0.005 0.01 0.015

 

milliii. .131.

. irisirlt

 

11\\~\\\|\Iu

Wiiiiitii

All\l\\|
.

I‘llllill‘ililrf

 

 

 

 

 

itlrfllffi\\\\\‘\

 

'\\\1\1I\\‘.\

II‘l'i‘r'i'r‘r‘

 

 

 

m.
b b

—0.015—0.01—0.(X)5 0 0.005 0.01 0.015

—0.015-0.01—0.(X)5 0 0.005 0.01 0.015

1:0.8

 

ilialoiiliii.

 

I I ‘II‘III

 

II‘I‘I‘I’I‘: 41

 

 

 

 

 

Wirxxx;;x
""’Ili'

wili"..rrl.'i'

 

   

{I'III'lvifilIi'ii'u‘iiii

V\\'\1~\\\1\

 

 

 

0015
0.01
0005

O

—0.005
—0 01

—0.015-0.01—0.005 0 0.(X)5 0.01 0.015

—0.015-0.01—0.005 0 0.005 0.01 0.015

Figure 2.11. These graphs show the velocity field around of a spheroid with periodic mo-

tion.

-0.005

-0.01

 

 

 

 

—0.015—0.01—0.005 0 0.005 0.01 0.015 —0.015—0.01—0.005 0 0.005 0.01 0.015

1:=0.4 r=0.6
0.015

0.015

0.01
0.005
-o.005 7 -o.005

—0.01 . —0.01

 

 

 

   

—0.015—0.01—0.005 0 0.(X)5 0.01 0.015 —0.015—0.01—0.005 0 0.005 0.01 0.015

0.015 1:=O.8 1:=1.0

0.01
0.005
0
—0.005

—0.01

 

 

   

 

—0.015—0.01—0.005 0 0.005 0.01 0.015 —0.015—0.01—0.005 0 0.005 0.01 0.015

Figure 2.12. These graphs show the pressure field around of a spheroid with periodic

motion.

45

slip coefiicient and the geometry to the motion of the particle. Different motions of the
particle are shown in Figure 2.7 and Figure 2.8. The predicted period of each motion is
shown in Table2.1. It can be seen that the motion of a particle is periodic and the period
becomes longer with the decreasing of the slip coeflicient. When the slip coeflicient goes to
sufficiently small, i.e., at y = 1.0, the particle rotates to a fixed angle and achieves a steady
state shown in Figure2.7 and Figure2.8 by fine dashed lines. For oblate spheroids(e > 0),
as the sliding friction decreases relative to the viscous friction (slip coefficient), the period
increases for a fixed aspect ratio; as the aspect ratio increases, the period increases for a
fixed slip coeflicient. For prolate spheroids (s < 0), as the sliding friction decreases relative
to the viscous friction (slip coefficient), the period increases for a fixed aspect ratio; as the
aspect ratio increases, the period increases for a fixed slip coefficient. As mentioned before,
the solutions with the slip by the present method is accurate to 0(52). For the slightly
deformed sphere the motion of the particle for when y —> 00, shown in Figure 2.7 by dashed
lines, agrees well with Jeffery’s [1]no-slip solutions. So the present method is effective for
slightly deformed sphere with consideration of slip boundary conditions. Figure 2.9 and
Figure 2.10 show the velocity and pressure fields of the case of steady motion of the particle
when )7 = 2.5 and s = 0.2 at different time. Figure 2.11 and Figure 2.12 show the velocity
and pressure fields of the case of periodic motion of the particle when y —> co and e = 0.2

at different time.

46

2.7 Summary

A perturbation method is used to solve for the motion of a fluid influenced by the presence
of a deformed sphere. Slip is assumed at the surface of the particle.The hydrodynamic force
and torque exerted by the fluid on the deformed sphere are expressed explicitly for a fixed
and rotating particle in a uniform streaming flow. Solutions to the limiting cases of non-slip
and perfect slip are identical to the existing solutions. The motion and orientation evolution
of a spheroid induced by a homogeneous flow are derived. Errors in the angular velocity
calculated by this method are of the order of 009). The period of rotation of the spheroid
is found to be longer as a dimensionless parameter that incorporates the slip coefficient
becomes small. When the slip coefficient becomes sufliciently low, the deformed sphere

rotates to a fixed angle and reaches to a quasi-steady state in the flow.

47

CHAPTER 3

MOTION OF AN ELLIPSOID IN

QUADRATIC AND CUBIC FLOWS

The influence of the physical parameters, such as the geometry of the particle, the viscosity
of the surrounding fluid, the slip coefficient, to the drag force and the motion of a deformed
sphere suspending in homogeneous shear flow has been studied in Chapter 2 and Chapter
3. The influence of the velocity field to the drag force and the motion of the particle will be

studied in this chapter.

3.1 Motion of an ellipsoid in a quadratic flow field

First, quadratic flow field is considered and the hydrodynamic force and torque exerted on
the surface of the spheroid by the surrounding fluid are determined. The same particle is

studied as that in the chapter 2 and chapter 3. The surface of the particle can be described

48

by the same function in the polar form as
r = r0[1 + ef(6)] + 009) (3.1)

in which
1 2
f(6) = — {3120(0089) + 519260860} (3.2)

The velocity of the particle can be decomposed by two parts, the translation velocity

and the rotation velocity.

ll = [10 + (r) X r
uxO wxz — only
= uyo + wzx — wxz (3.3)
uzo wxy - wyx

in which no is the translation velocity of the center of the particle, and to is the angular

velocity. The velocity field of the fluid far away from the particle is defined as

ill)

xy

vxo G11 G12 G13 I A111 A112 A113 A123 A122 A133 x2

V00: vyo + G21 G22 G23 y + A211 A212 A213 A223 A222 A233
V20 G31 G32 G33 Z A311 A312 A313 A323 A322 A333 yz

 

In spherical coordinates, denoted as (r, 6, (b), the position vector shown in (3.1), is ex-

pressed as

r = rsin6cos (flex + rsin6sin (trey + rcos 6ez (3.5)

Using the following definitions
x(0) = r0 sin 6cos (b (3.6)

49

 

 

 

Figure 3.1. This figure shows a spherical coordinate system.

y(0) = r0 sin 9 sin (b (3.7)

z(()) = r0 cos 6 (3.8)

and substituting (3.6-3.8) and (3.1) into (3.4), the velocity field can be decomposed as

 

 

W0 011 G12 G13 x(0)
(0)
V.» = VyO + G21 G22 G23 M0) (39)
Vzo G31 G32 G33 Z(0)
(
1%»
x(0)y(0)
A111 A112 A113 A123 A122 A133
x(0)z(0)
+ A211 A212 A213 A223 A222 A233
y(())z(())
A311 A312 A313 A323 A322 A333 2
Yb)
2
\ z(0) J

50

and

G11 G12 G13 x(0)
v2!) = f(6) 021 622 023 y(0) (3.10)
G31 G32 G33 Z(0)

 

 

i i
x30)
x(0)y(o)
4111 A112 A113 A123 A122 A133
x(0)z(0)
+2f(9) A211 A212 A213 A223 A222 A233
y(0)z(0)
A311 A312 A313 A323 A322 A333 2
y(0)
2
\ Z(0) i
and
v2? = 0 when i2 2 (3.11)
To satisfy the continuity equation of the fluid, it is required that
Tr(Vv) = 0 (3-12)
The following relations can then be obtained
011+ 022 + 022 = 0 (3.13)
(24111+4212 +A313)x = 0 (3-14)
(A112 + 2A222 +4323)x = 0 (3-15)
(4113 +4223 + 24333)x = 0 (3-16)

In (3.13-3.16), x, y, and 2 can be any arbitrary number. Hence the coefficient of x, y, and 2

should satisfy the following conditions

G33 = -(Gn + G22) (3-17)
A +A
4111: -—--—m 2 313 (3.18)

51

A112 +A323

A222 = ———-2-——— (3.19)
A +A
A333 = _ 113 2 223 (3.20)

If no-slip is assumed to be satisfied on the interface of the‘solid/fluid, then the appropriate

boundary conditions are

v=u on Sd (3.21)

and

v = voo on r —> oo (3.22)

where S d denotes the surface of the deformed sphere.

By using a similar method in chapter 2 [2], the three-scalar identities for a sphere with

no-slip are
°° nro L0)" 0') 1(1))" (1') :5. (i) 2
n;m[2(2n+3ifl(’ p, +rO r ¢,. r v (3. 3)
n(n+l)r()( :0)" (i) "(n‘1)("0)nn(i) (i)
__ _ V .24
:[2(2n+3)—_p(r) p" + to r _., v (3 )
Z n(n+1)(rr— °)xg)=r°VXV(i) (3'25)
nz—OO

When the velocity field is prescribed at the surface of the sphere, v(r0, 6, (1)) is a given func-
tion and each of the functions appearing on the rightside of (3.23- 3.25) may be expanded

in a series of surface harmonics as

E -v(ro. 6. ¢) = Z] w. 4) (3.26)
—rV - v(r0,6,¢) = Z Yn(6,¢) (3.27)
n=l

52

r - v x v(r0, 6, ¢) = Z 2,,(9, ¢) (3.28)
n=l

For exterior problems, the following relations can be obtained for n 2 0

2n—1

 

r0n+l
p_(,,,,, = (n + 1)r07 [02 + M. + m (3.29)

 

_ irfirfil
¢—(n+1) — 201+ l) r [an + Yn] (3.30)
1 r0n+l
X—(n+1)— n(n+ 1)}- Zn (331)

Substituting (3.29)-(3.31) into Lamb’s general solution for Stoke’s equation (2.3), the ve-

locity field of the fluids surrounding the particle may be written in the form

°° (n—2) 2 (n+1)
= V _ V _ —— V _ —— _
V ":1 ><(I‘X (n+1))+ ¢ (n+1) 2n(2n—1)pr P (n+1)+rn(2n_1)pp (n+1)
(3.32)
and
p = 226.1) (3.33)
n=l

By using (2.46) each component of the hydrodynamic force on the particle with con-

sideration of the no-slip on the interface can be obtained as

Fx = 67Tfo#(vx0 - ux0)(1 - +%8) + maria/1122 + 2A133 - A212 - A313)
+%7rr(3)}1(—14A112 - 12214133 + 714212 + 39A313)8 (3.34)
and
Fy - 67Tr0#(v)0 - “)0)(1 - %8) + mam-4112 + 24211 + 24233 - A323)

+31—57rr8p(7/1112 — 14A211 — 122/1233 + 39A323)8 (3.35)

53

and

Fz = 67"011(sz — uz0)(1 - %8) + ”ram-Ar 13 - 2A223 + 21“311 + 2A322)

+3l—57rrg/r(—2A1 13 — 211223 + 32A31j + 32A322)5 (3-36)

The hydrodynamic torque exerted on the particle by the surrounding fluids can be obtained
by using (2.47)

§[(—5023 + 5032 — 10wx) + 8(11023 - G32 + 12wx)l
TT = 4mg). %[(scl3 — 5631 — 106),) + s(—1 1013 + G31 + 12%)] (3.37)
(-Grz + G21 - 2%) + %8(012 - G21 + 240:)

If the particle is subjected to no external torques except those exerted by the fluid on its
surface, the torque will be cancelled out at each instant. From (3.37), the angular velocities

can be solved as

: 5(023 — G32) + (G32 - 11023»;

 

 

x -10 +123 (3'38)
_ 5(031-013)+(11013 - 031).»:
a), 7 —10 +128 (3'39)
1
(dz = 50012 + G21) (3.40)

From the equations of force and torque on the particle, it can been seen that the quadratic
term and the constant term in the velocity field of the fluids only result in force on the
particle while the hydrodynamic force is resulted from the linear term of the velocity field.

Compared with the homogeneous velocity field at infinity of the fluid in the chapter 3,

G can also be decomposed by the symmetric and antisymmetric part, i. g. S and W,

011 G12 G13 0 h g 0 -€ 77
G21 G22 G23 = h b f + 4’ 0 g (3.41)
G31 G32 G33 8 f C ~77 if 0

54

From (3.41), the angular velocity of the particle can be written as
w. = -f4 + 6. wy = g4 + (f. a». = r (3.42)
in which

_ 58
_ —5+68

 

(3.43)

This equation is exact same as (2.91) for the no—slip problem. The flow-induced rotation
motion of a particle suspended in quadratic flows is same as that the particle in homoge-

neous shear flows (constant gradient of the velocity field of the fluid).

3.2 Motion of an ellipsoid in a cubic flow field

3.2.1 The hydrodynamic resistance

To study the influence of cubic flow on the motion a particle, a cubic flow field far from a

deformed sphere is investigated here, i.e.

A1111 A1112 A1113 A1123 A1122 A1133 A1222 A1223 A1233 A1333
v80: A2111 A2112 A2113 A2123 A2122 A2133 A2222 A2223 A2233 A2333 344)
A3111 A3112 A3113 A3123 A3122 A3133 A3222 A3223 A3233 A3333
T
(x3 xzy x22 xyz xy2 x22 y3 yzz yz2 23)

Using (3.6)-(3.8), if Isl is a small number, the velocity field can be expanded by

A1111 A1112 A1113 A1123 A1122 A1133 A1222 A1223 A1233 A1333
A2111 A2112 A2113 A2123 A2122 A2133 A2222 A2223 A2233 A2333 345)
A3111 A3112 A3113 A3123 A3122 A3133 A3222 A3223 A3233 A3333

.53) =

T
2 2 2 2 3 2
(rig) x(0)y(0) x(0)z(0) x(0)y(0)z(0) x(0)y(0) x(0)z(0) y(o) y(20)z(o) y(o)z(0) Z3»)

55

A1111 A1112 A1113 A1123 A1122 A1133 A1222 A1223 A1233 A1333
(1).
Voo -3f(9) 42111 A2112 A2113 A2123 A2122 A2133 A2222 A2223 A2233 A2333
A3111 A3112 A3113 A3123 A3122 A3133 A3222 A3223 A3233 A3333

T
2 2 2 2 3 , 2 3
(R30) x(0)y(0) x(0)z(0) x(0)>’(0)Z(0) x(0)y(o) Homo) no) rimzm) )(0)z(0) Z(0)) (3.46)
v52 =0 when 1'22 (3.47)

For incompressible flow, to satisfy the continuity equation, it is required that

A2112 = -(341111+A3113) (3-48)
A1122 = -(3A2222+A3223) (3.49)
A2233 = -(3A3333+A3113) (3.50)
A3123 = —2(A1112+A2122) (3-51)
A2123 = -2(Arir3+A3133) (3.52)
A1123 = -2(Az223+A3233) (3.53)

With a no-slip boundary condition, the hydrodynamic force and torque exerted on the

particle are

F = 0 (3.54)

4
7",, = Enrng—S + new] 13 + (—5 + 11.9)A2223 + (—15 + 638)A2333 (3.55)

+ (5 - (9)1431 12 + (15 + 38)A3222 + (5 - 118)A3233 + (—50 + 608)w—:]

’0
4 5
Ty = —§§7rr0p[(-5 +118)A1113 + (-5 +118)A]223 + (—15 + 638)A1333 (3.56)
(U
+ (15 — 38)A3111 + (5 - 8)A3122 + (5 — 118)A3133 + (50 - 6089);]
0

56

4
TZ = Enrng—S + 3s)A1112 + (—15 + 9e)A1222 + (—5 + 138)A1233 (3.57)

+ (15 — 98)A2]11 + (5 — 38)A2122 + (5 — 138)A2133 + (-50 + 308)—w2z]
r
0

One possible induced angular velocity of the particle if no external torque applying on the

particle can be obtained

,2
a), = ——°—[(-5 + new-113 + (—5 + 11mm, + (—15 + 638)A2333 (3.58)
50 -— 608
+(5 - 8)A3112 + (15 + 38)A3222 + (5 - 118)A3233l
,2
o, = Eff—“Ema + 11.9)A1113 + (—5 +118)A1223 +(-15 + 63s)A1333 (3.59)
+(15 - 38)A3111 + (5 - 8)A3122 + (5 -118)A3133l
r2
(oz '2 3-()—:()3—Og[(—5 + 38)A1112 +(-15 + 98)A1222 + (-5 +138)A1233 (3.60)

+(15 - 98)A2111 + (5 - 3642122 + (5 -138)Azr33l

3.2.2 Motion of a deformed sphere in a simple cubic flow

The motion of a deformed sphere in simple cubic flow, shown in Figure 3.2, is investigated.
Described in the fixed coordinate system (x’, y’, z’), the velocity field far from the particle

is assumed to be

v;.=(o 0 Ky’3) (3.61)

Induced by this cubic flow, the deformed sphere rotates around x (x’) by angle a. In
Figure 3.2, (x’, y’,z’) is the fixed coordinate system and (x, y, z) the rotating coordinate

system attached on the axis of the particle. Between these two coordinate systems, the

57

 

 

 

 

 

Figure 3.2. Illustration of a cubic flow field surrounding a deformed sphere.

following relation can be satisfied

x’ i l 0 0 x
y’ = 0 cos a — sin a y (3.62)
z’ 0 sin (1 cos a z

According to the (2.78) between these two coordinate systems, i.e., the fixed coordinate
system and the rotating coordinate system, the velocity of the surrounding fluid far away

from the particle can be described in the rotating coordinate system as

. 1 0 0
Voo = O cosa sina vfx,
0 —sina cosa
1 0 0 0
= 0 cosa sina 0 (3.63)
0 —sina cosa Ky’3

58

Cosa v2=Ky3, 5:0.1, r0 =0.01, u=0.02 — K=1-0

 

1
- K=2.0
0.5 , --- K=5.0
tx 10000

 

 

Figure 3.3. Evolution of the cosine of the rotation angle of the particle induced by simple

cubic flows at different constant K s.

From (3.62), then

0

Voo = K sin cr(cos3 ay3 — 3 cos2 or sin ayzz + 3 cos a sin2 aryz2 -‘ sin3 az3) (3.64)

cos a(cos3 oy3 — 3 cos2 a sin ayzz + 3 cos a sin2 aryz2 — sin3 az3)

Comparing with (3.46), the following is obtained

0 0 0 0 0 0 0 0 0 0
Voo = K 0 0 0 0 0 0 sinacos3a —3cos2asin2a 3cosasin3a —sin4ar
(O O 0 O 0 0 0084a —3cos3asina 3cos2crsin2cr -cosasin3a
T
(x3 xzy ”22 3% xyz XZZ )3 y2z yz2 Z3) (3.65)

From (3.58)-(3.60), components of the angular velocity of the particle are

3Kr(2)[5 — 8(11 — 10cos26)]

 

w" 7 “7 —50+60e (3‘66)
wz = 0 (3.68)

59

COSO’ v'z: K)’, 620.], I0=0.01, ”=0.02 _ K‘l 0

- - - - K=2.0
t --- K=5.0

 

 

 

Figure 3.4. Evolution of the cosine of the rotation angle of the particle induced by simple

shear flow at different constant K s.

Integrating (3.66), the angle a is obtained

_ / 5-8 3Kr0 file—5'.
a—arctanl _5+218tanh[10 _5+8 (3.69)

The evolution of cosine of the angle a is shown in Figure 3.3 at different Ks. It can be seen

 

 

that the induced motions of the particle by a cubic flow is periodic. It is known that the
induced motion of a deformed sphere in a homogeneous shear flow field is also periodic.
In order to compare the difference of the two motions induced by a cubic flow and a simple
Table 3.]. Predicted periods of the particle motion induced by cubic flows and simple shear

flows.

 

v2 = Ky’3 v’ = Ky’
K=1.0 K=2.0 K=5.0 K=1.0 K=2.0 K=5.0
period(s) 6000 3100 1250 12.50 6.40 2.56

 

 

 

 

 

 

 

 

 

 

 

 

 

 

60

shear flow, the period of each motion at different Ks are listed in the Table3.1 when other
parameters are e = 0.1,r0 = 0.01, and p = 0.02. The periodic motion of the particle
induced by the simple shear flow is shown in the Figure 3.4. From the Table 3.1, it can be
seen that at the same K, periods of the induced motion. by a cubic flow are much longer
than those by simple shear flows. For both cases the period is proportional to the inverse of

the coeflicient K.

61

CHAPTER 4

POWER-LAW MODEL OF A

DEFORMED SPHERE

4.1 Introduction

Newtonian fluids are investigated in the first four chapters, in which the viscosity is as-
sumed to be a constant. In several industrial problems, polymeric liquids are involved and
their viscosity is often dependent on the shear rate, temperature, pressure, etc. This chap-
ter is devoted to studying the influence of a non-Newtonian viscosity on the motion of a

particle.

4.2 Power-Law model for the Non-Newtonian viscosity

One of the earliest empiricism for Non-Newtonian fluids is based on the modification of

Newton’s law of the viscosity in which the viscosity is allowed to vary with the shear rate.

62

For example, for an arbitrary incompressible Newtonian fluid, v = v(x, y, z), the constitutive

equation is modelled as :
t = w? (4.1)

in which a is a constant for a given temperature, pressure, and composition, 7 the rate-
. 1 . . . . .
of-strarn tensor ~2-[Vv + (Vv)T]. To include the 1dea of a non—Newtonian vrscosrty, the

constitutive equation is modified by
t : —7n'r (4.2)
in which n is a function of the scalar invariants of "y.

There are three invariants for a second order tensor 7 as the following:
I = Z 7’17 (43)
i

11 = Z Z i’iji'ji (4.4)
I

i
m = Z Z 2 2mm.- (4.5)
i j k
For an incompressible fluid I = 2(V - v) = 0. For shearing flows III turns out to be zero;
Because (4.2) should be used only for shearing flows, or at least flows that are very nearly
shearing, omitting [H from any further consideration is not a serious restriction. Hence I)
is taken to depend only on H. In practice, 7, the magnitude of the rate-of—strain tensor 9, is

often used instead of II, i.e.

 

. 1 . . 1
7= Jigg'fifl’ji = (fill (46)

Experimentally, a typical viscosity vs shear-rate curve is shown in Figure 4.1. It is
composed of two regions, a zero-shear-region and a power—law-region. In the zero-shear-

region (low shear rate), the shear stress is proportional to y, and the viscosity approaches a

63

Zero-shear-region Power-law-region

A A

T

V

 

n(Pa-s)

 

 

V

44)

Figure 4.1. Schematic of a typical viscosity variation w.r.t to the shear rate.

constant value 770, the zero-shear-rate viscosity. At the power-law-region (the higher shear
rate), the viscosity decreases with increasing shear rates.

Two models for the viscosity 27(7) are often used to describe the viscosity in term of
shear rate [77]. The first one is the Carreau-Yasuda model. The Carreau-Yasuda model is a
five-parameter model, which has sufficient flexibilities to fit a wide variety of experimental

77(7) curves. The model is

’i — 77 . _
—°i =11+(4y)“1<" 1”“ (4.7)
’70 7 7700
Here 770 is the zero-shear—rate viscosity, 7700 is the infinite-shear-rate viscosity, xi is a time
constant, it is the ”power-law exponent”(since it describes the slope of 5,7647% in the

”power-law region”, and a is a dimensionless parameter that describes the transition region

between the zero-shear-rate region and the power-law region. In most industrial problems,

64

the descending linear region (the ”power-law-region”) shown in Figure 4.1 is more impor-

tant. The titled straight line can be expressed simply by

,7 = myt-l ~ (4.8)

which contains two parameters: m (with units of Pa - s"), and n (dimensionless). (4.8) can
be regarded as the limiting expression for high shear rates obtained from (4.7) with 7700 = 0.
When n = 1 and m = p, a Newtonian fluid is recovered. When n < 1, the fluid is said to
be ”pseudoplastic” or ”shear thinning”. When n > 1, the fluid is called ”dilatant” or ”shear

thickening”.

4.3 Solution to flow problems using a Power-Law model

Since the viscosity depends on the strain rate of the flow field, the Stokes’ equations to
govern the motion of the fluid will be nonlinear. It is not feasible to get an analytical
solution for such nonlinear equations. Assume that the viscosity of the fluid is locally
constant. Solutions obtained from the linear Stokes equations (with constant viscosity)
can be used to get the hydrodynamic force and torque on the particle approximately (see

Chapter 2 and Chapter 3). The hydrodynamics force on the surface of the sphere is given

szfr-ds (4.9)

where 1: is the stress dyadic and (18 is a directed element of surface area parallel to the outer

by

normal direction. By using the divergence theorem, the hydrodynamic force exerted by the

65

fluid can be written in the form

F: f PrerQ (4.10)
51

where P, = 1: - r/r is the stress vector acting across the surface, r = constant and
d9 = sin 6d6d¢ is an element area on the surface of a sphere of a unit radius, S1. The

hydrodynamic torque applying on the surface of the sphere is
T=ert-r/rdfl (4.11)
Substituting (4.2) and (4.8) into (4.10) and ((4.11), the force and torque will be
F = me-ly - mm (4.12)
T = fr x "ii/'71? - r/rdo (4.13)

(1) A fixed sphere in a uniform stream flow
A sphere is fixed in a space and a uniform streaming flow passes by the sphere at the
velocity U ex. By using (4.12) and (4.13), the hydrodynamic force and torque on the
sphere are

( 9m7rU3
8r0

 

81an5
1 320,3

 

" = 5 (4.14)

2187an7
[ 35840rg

 

and

T = 0 (4.15)

66

(2) A fixed spheroid in a uniform stream flow
When a spheroid with deformation coefficient a is fixed in a uniform flow at velocity
Uex. The geometry function of the spheroid is same as (2.49). By using (4.12) and

(4.13), the hydrodynamic force and torque can be obtained when n = 3

3 3
F: 9an +13.4181’"”U
8r0 8r0

 

 

s + 0(32) ex (4.16)

and

T = 0 (4.17)

If the velocity field of the fluid at the infinity is Wez, when n = 3 the hydrodynamic
force and torque are

_ 9m7rW3 mer3

F — - 0.9366
8m 870

 

 

e + 0(52) eZ (4.18)

and

T = 0 (4.19)

(3) A free spheroid in a homogeneous shear flow
The far away velocity field of the fluid is assumed to be simple shear, e. g. v’ = K y’ez,
which is described in the fixed coordinate system. Suppose that initially the long
axis of the particle is parallel to the z’ axis [seen Figure 2.5]. Due to the nonlinear
property of this problem there may be multiple solutions for the induced motion of
the spheroid. In order to make this problem solvable a hypothesis is introduced here,
which is that the particle can rotate around the x axis and the other angular velocity

coy 2 (02 = 0. According to this hypothesis the distortion and the rotation of the fluid

67

described in the rotating coordinate system which is attached on the particle can be

written in the form

0 0 0
S: 0 b f (4.20)
0 f —b-
and
0 0 0
W: 0 0 —§ (4.21)
0 g 0
in which
(4.22)

1 1
b = -2—Ksin2ar, f = §Kcos2a, andtf = —K/2

The hydrodynamic force and torque exerted by the surrounding fluid on the spheroid

can be obtained as

 

"”"68 3 2 2 2 2 2
F x = 18018[(292328K +569868K wx+1287K(253wx—149wy)+34749wx(wx+wy)]
(4.23)

mzrr2

F, = j—Q[388K25wy + 3(o3, + obese), + 14wx10g 2)

+K(4028(oxwy + 56(6)}, + o§)16g 2] (4.24)

F - —mm‘2’8wy {+8316[(7897088+214865 ) 2 102400 2]
Z 7 441548800 ’r a" “y ”x

+K2(158195474432 + 572984890570

+924K (222636032 + 706159570} (4.25)

and
[—bm7rr33(33938801875b2 + 3(14079425275f2— (4.26)

—7702408035f(§ — 6.)) +3976823046(§ — wx)2)] /386260875

68

7",. = 0 (4.27)

T2 = 0 (4.28)

If there is no external torque applying on the spheroid, two possible induced motions

of the spheroid are

 

w, = -O.9684 f :1: 1.6132 ,/(—1.093032 - f2) + g (4.29)
(oz = 0 (4.31)

Substituting (4.22) into (4.29), it can be obtained

 

(ox = —0.9684 f :1: 0.8066 J(—1 - 0.093 sin2 26)— K/2 (4.32)

From above equation it can be seen that the induced angular velocity of the particle is a
complex number. The reasons to result in a complex angular velocity may be from two
hypophyses: one hypophysis is that the local viscosity is assumed to be constant when
solving the Stokes’ equation by Brenner’s method; and another one supposes that the par-
ticle just rotates along the axis x. Due to the complexity of this problem, some other model

to describe the property of Non-Newtonian flow will be sought in the future work.

69

CHAPTER 5

MOTION OF A SPHEROID IN AN

OSEEN FLOW

5.1 Introduction

It is well-know by the Stokes’ law, that a force of 67rrrWa is required to maintain a uniform
velocity W of a sphere of radius a moving through a liquid with viscosity 11. However,
this result can not apply to the cases which arise in practice to deal with, for example, the
motion of ships. A mathematical incompleteness of the solution is because the advective
terms are not negligible compared to the viscous terms at large distances. From Chapter 2,

the large viscous term is of the order

. . Ua
vrscous force 2 stess gradrent ~ ’1—3 as r —> 00 (5.1)
r

while the largest inertia force is

. . (9a U 2a
mertra force ~ pur——6- ~ p
6r r2

 

as r —9 00 (5.2)

70

Therefore

inertia force pU a r r
. ~ —- = Re— as r ——> 00 (5.3)
vrscous force it a a

 

It can be seen that the inertia forces are not negligible for distances larger than r/ a ~ 1/Re.
The neglected terms become arbitrarily large at sufficiently large distances, on matter how
small Re may be [78].

Several attempts have been made to correct this error. In 1893, Rayleigh introduced
some additional forces to improve the accuracy of the Stokes’ solution. In 1911, Oseen
proposed a modified system of equations, in which the inertia terms are partly taken into
account, and obtained the solution for flow past a fixed sphere using these equations. Os-
een’s solution is satisfied at infinity; it also gives good approximation near the sphere if the
velocity or the radius of the sphere is small. Lamb [6], used a different method to present
Oseen’s solution, but it still can not overcome the restriction of Oseen’s. Burgess [3] used
a simpler method than that of either Oseen or Lamb. The solution from Burgess method
can satisfy the boundary conditions at infinity and the boundary conditions on the sphere
can be satisfied to any desired degree of approximation. Afterwards, solutions of Oseen’s
equation have been extended to the spinning sphere [79], circular cylinder [80—83], circular

and elliptic cylinders [84—86] and flat plate [84, 87].

5.2 Analytical structure of the Oseen flow

Derived by Burgess [3], the motion equation of a viscous flow in the cylinder coordinate
system is

(—W—6— — vD)D¢ = O (5.4)
62

71

in which W is the velocity of the particle v = 11 /p is the kinetic viscosity of the fluid, and
62 1 a 62

D a — — —— + —.
6r2 I“ 6r 622
Due to the commutativity of this equation, a solution can be given as
ii! = W + ill" (55)
where
Dt/x’ = 0 (5.6)
and
(9 H
(vD + W—)(// = 0 (5.7)
6z

Solving above partial differential equations by using a perturbation method, the stream

function of the Oseen flow is obtained by Burgess as

. - 1 3k 3
ill-'7- e7"p(1+C°SQ) C+Dcos6+b1sin26(k+;)+b2sin26cos6(k2+—+-—2-)+---
p p
M sin26 Nsin26cos6 Psin26
+ 2 + 3
p P .0
= sin71

. . _ Z r
lIlWhIChpZ Vr2+z2 andgzcosl 2 2 __

function can be fixed by suitable boundary conditions. From the stream function the radial

(500826—1)+--- (5.8)

 

 

 

Lcos6+

 

. Constants in the stream

velocity and tangential velocity of the fluid can be easily obtained as

 

 

vp : p2 slin 6% (5.9)
and
= _p sin 62—: (5.10)
The normal and shear stress in the fluid are
6vp
O'pp =2l1'a)‘ (5.11)

72

 

 

16v v
Ugg=2y(- 672+ 5) (5.12)
1 3V¢ vacot6
:2 .l
0"” ”(psin06tfi +1 p p+ p ) (5 3)
a v 16v
Tp6: 2”(pap —(— 6’)+ p 78;) (5.14)

Two kinds of boundary conditions can be applied on the interface of the particle and

fluids.
(1) No-slip Boundary conditions on the particle
vp = 0, v6) = 0 at Infinity
vp = 0, v49 = 0 on the surface of the particle
(2) Slip Boundary conditions on a symmetric particle
vp = 0, v6) = 0 at Infinity

v,D = up and V9 — 149 = firpg on the surface of the of the particle

5.3 Applications of Burgess’s solution

A simple form of stream functions of Oseen flows is used as

M sin2 6
p

 

(j, = e7kp(l+°039)b — 0(1— cos 6) + Lcos6 + (5-15)

A symmetric particle moves in a quiescent fluid with a velocity W in the z direction with
different boundary conditions and different shapes of the particle. Described in Figure 5.1

the radial and tangential velocity velocities on the surface of the particle are
up = Wcos6,u9 = —Wsin6 (5.16)

73

 

 

V

 

Figure 5.1. Description of the surface velocity of the particle.

Applying slip and nonslip boundary conditions and taking account of different shapes of
the particle, several solutions are obtained as the following:
Case I Uniform motion of a sphere in Oseen flow with nonslip boundary condition

Solutions of this case have been already given by Burgess as [3]

3aW _ Wa3

W'M — 4 (5.17)

L=b0=-

in which a is the radius of the sphere, k = 2p/p. The velocity field of the fluid around the
sphere with no-slip boundary conditions is shown in the Figure 5.2(a).
Case 11 Uniform motion of a sphere in Oseen flows with slip boundary conditions
If the production of ak is a small number, i.e. a is small or/and k is small, the term

(“MM“)8 9) can be expanded by series with neglecting the higher order terms of ak as
(“Wm“) z 1 — ak(1 + cos (9) (5.18)

74

 

 

4444
4444

4
Abltpgp

444

1444
0444(1‘11

‘

b

wriWI'

4
4

.b bib

5.5.8.. r4.4.4-l4444 5.185
DDLLDD: v4llllIO‘ArALDD
LDDDDDim m‘444444‘AfADLD
DDDDDLLO .444444445DLDD
DDDDDDA. '44411444A_DLDD
DELL D05” wrl“““ll DDDD
Abbbhheo. .4444444‘JAAAK

.‘Dfil ~444“.‘If‘11‘1“D‘
‘1...“ 4‘1 .‘0-

l
l

Ibbbbggpg‘

4

fi.bh§}\"

‘Abbssst

4|||§|§;.
0111201114

0

 

 

rurrtbtl

"W-VrYfi—w

 

 

’1sstirril

f
g
'0
(b)
9
1
A
f
f

 

 

 

-0.004 -0.002
(d)

 

 

 

0.(m’.\§;ltttl

0.1); (c) no-slip boundary condition apply-

 

 

44444444
‘4414414
44441441l

I
4
4

‘O44A .I-‘ ‘O‘A
4444.. vthlrviibiblbl 64.4.4444
4444:” .AAADAAAA4624444IA
44444430. “tbtbhbbbA44444444
4444441._ vbfibbhbbhl.444444“
444444;. .lblllbblx.|4444444
444444;“. wlbbbhllll.44444444
4““‘:$ TD‘DDDDDDA-‘44O‘l‘l
444444: vbbbfibbbfixl44444444
m m m m o m m
0. 0. 0 0. 0.
_ _ a _
b b b b I b e 4 4.4 l 4 4 4 4 .- .AVJL IiHib I. L
bbhlblx.w #444444441filhbbllbl
hlbbbb: . w444444441 Abhbbbbb
bbbbblto r444444441.lllbbblh
bbbbbh: .444444445blbbbbbh
bhllb|1.m .4444444445‘5hhh5b
.
5351050. _4444444.4.4.§.|.§lbbt|
DOD-I4. .434fi‘lr‘l‘l 133R‘ADL
-4031 v 0.4“

i
4/

0.2); (d) slip boundary condition applying on a deformed
75

0112 0.(X)4

0
(c)

 

 

| A A 1
t A A 4
D |

slbl
ikbbl

 

 

"W-vvv—v—v
’rsstt

b

(a)

l

l
f ,
it

,

 

 

 

)

r 4 4
=0.1).

—0.(I)4 —0.(X)2

 

444144444

 

 

0.CX)4 .

‘I‘r ‘. ‘ 4 4 I: _.I 4
64.4.4444... .Lliiiibiblbl 4.4.
‘ ‘ ‘ 5% .D D D D. A A L. D. 4 4.. 4. ‘ 4
4 4 4 l ‘ l 50. v A D L D L D D O A 4 4 0 4
4 4 4 ‘ ‘ J A. _ v D A D D L D D L 1.4 4 4 4
4 4 4 4 4 4 41 -D D L D D D D A A A 4 4 4
4 4 4 4 4 4 :m wk 5 b b h I l A 1. 4 4 4 4
4 l l l l l :O— v D L D D D D D D 1- 4 4 4 l
l 4‘ l l l l A. V D D U | D b D | LI ‘ ‘ l ‘
m m m m o m
4 x. o o. ._.

Figure 5.2. Velocity vectors of the surrounding fluids with different boundary condition
on the surface of a particle. (a) no-slip boundary condition applying on a sphere; (b) slip

boundary condition applying on a sphere(fi
ing on a deformed sphere(e

sphere(e = 0.2,,8

—0.004HH“‘“

Substituting (5.15) into (5.9), (5.10), and (5.14), the radial and tangential velocities, and
the shear stress of the fluid on the surface of the particle with neglecting the higher order

terms of ak can be written as

 

 

L 2Mcos€ b
v,, [Fa z —a—2 + ——a3—— + 5%(1- 2akcos 6) (5.19)
MsinO b ksin6
V.) [W z 03 + 0 a (5.20)
6 MsinO
Tp6lp=a z "—# a4 (5.21)

Applying the slip boundary conditions on the sphere to determined the constants in the

stream function, the following relations can be obtained

 

L b 2M 2bk
Wcosaz-—+—°+(————°—)cose (5.22)
02 02 a3 a
M ' 0 b k ' 0 6 M ' a
8;“ + 0 5‘" +Wsin9=-——” :1“ (5.23)
a (1 ,Ba

Equate the constant terms and the coefficients of cos 0 and sin 6 and each side of (5.22) and

 

 

(5.23 ), then
f M M
—3 + b—Olf + W = -6fl’84
a a a
i L : b0 (5.24)
2M
_ _ 2b0k = W
t a3 a

 

Solving the above equations, the constants of the stream function are

_3(a2W,B + 2am.) _ _ a4w,6
4k(aB + 3p) ’ — 4(a3 + 3p)

 

L = b0 = (5.25)

When ,8 —> oo, slip boundary conditions change into the nonslip boundary conditions.

From the solution above letting 6 ——> 00 will result in the following constants

3aW _ Wa3

L: :——— _
1’” 4k ’ 4

 

(5.26)

76

These constants are exactly same as those solved by using nonslip boundary conditions.
The velocity field of the fluid around the sphere with slip boundary condition is shown in
the Figure 5.2(b).
Case III Uniform motion of a slightly deformed sphere with no slip on the interface
A spheroid, regarded as a deformed sphere, moves in a viscous fluid with the velocity

W in the z direction. The shape function of a spheroid is

xz+y2 :2

a2 + m = (5.27)

If s is a small number, the shape of the spheroid can be described in a polar form with

neglecting the terms of 0(82) as
p = a(1 — ecosz 6) (5.28)

For the nonslip boundary conditions by using Burgess’s solution, it can be obtained

 

 

 

 

L 2 b0 = —%a:—IQ;T? (5.29)
and
3
M = _a4(_V:’(—1_;8_;:) (5.30)
Expanding these constants by series of a, the constants can be expressed as
Lzb0:_3aW_3aW8_3aW32+0(83) (5.31)
4k 4k 2k
and
M ___ _a3W _ 3a3We _ 3a3W.‘:2 + 0(83) (5.32)

4 4 2

From (5.31, 5.32) it can be seen that when a —> O, the constants solved for this case are

exactly same as the constants solved for the nonslip sphere problem (Case I). The velocity

77

field of the fluid around a deformed sphere with no-slip boundary condition is shown in the
Figure 5.2(c).

Case IV Uniform motion of a deformed sphere with slip on the interface

Applying the slip boundary conditions on the interface of the spheroid, the constants

can be obtained

 

 

 

 

 

 

2

L : b0 = 214—234:V fig—2g: )+ 43p), (5'33)

and

3
=.:J.‘;::i:;ii;f::. (532
When ,8 —> 00,

L = b0 = “213—313;” = 3:3] — 322,8 — 312:2 + 0(83) (5.35)

and
M 2 a3W(1 +8) _ a3W _ 3a3W8 _ 3a3W£2 + 0(83) (5.36)

2(—2 + a) ‘ 7 4 8 16
Because both ak and 8 are small numbers, the errors between these solutions and the solu-

tions for case III are negligible.

Let ,8 —> co and s = 0. The constants will be

3aW Wa3
“b0: ‘17”: 4

 

(5.37)

which are the same constants for the first case. The velocity field of the fluid around a
deformed sphere with slip boundary conditions is shown in the Figure 5.2(d). It can be
seen that the difference of the velocity field around the particle with considering slip and

no-slip boundary conditions is small.

78

5.4 Summary

Uniform motions of a particle through a viscous flow are solved analytically by using
Burgess’ general solution for the Oseen flow. Nonslip and slip boundary conditions are
considered on the interface of the particle and the fluid respectively. TWO kinds of geom-
etry of the particle, e.g. a sphere and a slightly deformed sphere, are studied. Four cases
are calculated respectively according to different boundary conditions on the interface and
the shape of the particle, e.g. (1) the motion of a sphere with nonslip, (2) the motion of a
sphere with slip, (3) the motion of a deformed sphere with noslip, and (4) the motion of
a deformed sphere with slip. Both of the solution of the case (2) and case (3) can recover
the solution of case (1) when letting the slip coeflicient )6 and deformation coefficient 8
equal to zeros correspondingly. The error between the solution of case (4) when the slip
coeflicient goes to co and that of case (3) is negligible if when the velocity and the diameter
of the particle are small. The boundary condition at the infinity are satisfied well and the
boundary condition at the interface can be approximated satisfied if the length dimension

of the particle or the velocity of the particle is small.

79

CHAPTER 6

THE FLOW-INDUCED

ORIENTATION OF RIGID

PARTICLES IN DILUTE

SUSPENSIONS

6.1 Introduction

The prediction of the flow and orientation of suspensions during the processing of com-
posite materials is important to understand how the processing conditions influence the
mechanical properties of the final part. A variety of models and algorithms have been
made to derive the relationships between processing conditions and orientation of fibers

[39, 73, 88—97]. If the particle is axisymmetric and the diameter is much less than the

80

length of the particle, a unit vector p, collinear with the long axis of the particle, can be
used to represent the orientation of the particle. Only two Euler angles are related with this

vector p by

p = (sin 6 cos ¢ sin 6 sin (12 cos 6)T. (6.1)

Two methodologies are often used to describe the flow-induced alignments of particles.
One is based on an orientation distribution function N, which is a function of Euler angles
and time, and the other one uses an ensemble average orientation tensor < pp > or < R >.
For spheroids, the orientation distribution function N can be defined over a state space of
orientation vectors p and N is also related to the aspect ratio 'of the particle, i. c. 2. For
ellipsoids, N can be defined over a state space of rotation operators R, and N depends on
two aspect ratios 2 and S as well. Comparing with the distribution function method, the
averaged orientation tensor method is often preferred but a closure problem arises due to
the averaging procedure [98,98—103].

Efforts have also been focused on the motion of non-axisymmetric particles suspended
in a slow viscous flow. The classic work pertaining to the motion of a particle, e.g. an
ellipsoid, in a uniform shear flow can be tracked back to the work of Jeffery’s [1]. In his
paper, the behavior of an ellipsoid suspending in uniform shear flow field is analyzed on
basic of Stokes’ equations of motions. Some previous authors [8, 69—71, 104] extended
Jeffery’s work to a more general shape particle.

The orientation distribution function for rigid ellipsoidal particles in a simple shear

flow was given by Workman and Hollingsworth [105] in terms of a series of spherical

harmonics by expressing the coeflicients in these series in the Feenberg perturbation form.

81

In both Brenner [71] and Workman’s (1969) work [105], Euler angles are used to represent
the orientation of the particle. Rallison [106] proposed a rotation matrix to represent the
orientation of the particles instead of using the Euler-angle representation. Considering the
influence of Brownian motions, Rallison [106] derived the time evolution of the orientation
probability distribution on the basic of Fokker-Planck equation and obtained the form of the
orientation probability distribution for small departures from an isotropy state. The second-
order moment of the probability distribution was also developed in Rallison’s paper [106].

In this chapter, a new closure model is developed for the motion of rigid particles of
complex shapes. Each particle is non-axisymmetric and its orientation is described with a
second order tensor < R >. An evolution equation for the second moment of the distrib-
ution function, which forms a fourth order tensor < RR >, is used in order to obtain the

average orientation of the particles in homogeneous flows.

6.2 Prediction of orientation of axisymmetric particles

Due to large amount fibers encountered in many applications, an approach based on a
distribution function is preferred to predict fiber orientation. The distribution function is
¢(p,t) where p is a unit vector along the axis of each fiber to represent its orientation.
The governing equation for (Mp, t) depends on the conservation of fiber orientations. In

homogenous flows with neglecting Brownian motions, it is given by [95]

Dill__i. .
Dt — 6p (W) (6-2)

The flow-induced motion of a single spheroid in the uniform shear flow is derived by

82

Jeffery [1] as
I’J=—W'P+4(S°P—85PPP) (6.3)

where W is the vorticity tensor and S the strain rate of the flow fields. With considering
the rotary diffusion due to the particle-particle interaction, the following hypothesis is used
[77,95,107]

6w

(P - 151W = 43116—1)- (64)

in which DR is the rotary diffusion coeflicient. Combining (6.2) - (6.4), the evolution

equation for the distribution function can be derived as

D¢_ a . a 6.1.
—— 6p [W p+/I(S p S.ppp)]+ DRap

Dt 6p (6.5)

This approach can represent the exact and full solution for fiber orientations [108, 109].
However, this partial-differential-equation is not analytically solvable for a general prob-
lem. Furthermore, using numerical calculations this method is, still too complex when
solving three-dimensional fiber orientational in complex geometries.

An alternative approach is to use a more compact description of the distribution func-
tion. A second-order tensor a is often used to represent the fiber orientation state at any

point in the material. This tensor is defined as

a = fpw<p>dp (6.6)

where the integral is taken over all possible directions of p. A fourth-order tensor can be

defined in the same way as

A = f ppppmmdp (6.7)

83

The second-order orientation tensor is symmetric (aij = a jg), and it has a unit trace (“ii =

1). Applying the similar symmetry conditions to the fourth-order tensor A, it is shown that
Aijkl = Ajikl = Akijl = Alijk =‘Aklij, etc (6-8)

The fourth-order—tensor provides complete information about the second-order—tensor be-

cause

a.‘ j = Aijkk (69)
Combining (6.6) and (6.5) gives the evolution equation for the orientation tensor,

g;=—W-a—a-WT+.1(S-a+a-S—2S:A)+2DR(6—3a) (6.10)

The advantages of using a tensorial approach [73] are that a is independent of the co-
ordinate systems and easily transforms between coordinate systems; it can be measured
by direct experiments; the tensor representation is extremely compact, and computational
efliciency. However, a disadvantage of this approach is that an unknown fourth-order ori-
entation tensor is introduced in the momentum equation since only the second-order tensor
is used to represent the orientation state. Closure approximations must be used to describe
the fourth-order orientation tensor.

One simple approach is to approximate the fourth-order tensor in terms of the second-
order tensor. Several closure models have been developed to predict the orientation of
ensemble particles, such as the quadratic model [98], the linear model of Hand [99], Hinch
and Lea] model [98], Hybrid model of Advani and Tucker [100], Orthotropic closure mod-
els [101], and the fully symmetric quadratic model [102, 103]. All these models can be

applied to fibers, spheroids, and other axisymmetric particles. The time evolution of the

84

 

 

 

 

I'D

Figure 6.1. Description of the Euler angles used in this chapter.

second-order orientation tensor < pp > is solved to obtain the average orientation of the

particles.

6.3 Predictions of orientation of non-axisymmetric parti-

cles

For rotating particles of arbitrary shape suspended in the fluid domain, a vector associated
with the particle can be mapped from the reference configuration to the current configura-

tion by using a rotation matrix R, which is defined by

R = mm“ + qmq" + r<t>r° (6.11)

85

}, R(¢,6.w) ‘4
,L j .4: (0° q. .°).__" (9° q. .»)*<' 3*
R7049”) ' '1" Y

Figure 6.2. Mapping procedure of a vector associated with the particle between the refer-

ence configuration and the current configuration.

where p, q, and r are three orthogonal axes of the particle in the current configuration and
p0, qO, and r0 are the corresponding axes in the reference configuration shown in Figure
6.1. Since

R-erl win

a vector can also be mapped from the current configuration to the reference configuration.
The mapping procedure is shown in the Figure 6.2. This rotation matrix R depends on the

three Euler angles by the following relation

—sinr// +cosr/1 0 +0036 0 —sin6 +sin¢ —cos¢ 0
R: —cosr// —sinr// 0 0 +1 0 +cos¢ +sin¢ 0 (6.13)
O 0 +1 +sin6 0 +cosl9 0 0 +1

The ranges of the three Euler angles are

OS¢SZN (6M)
os¢szn (aw)
OSBSH (6%)

Instead of using Euler angles, the orientation of the particle is represented by a rotation

matrix R, which can be written in the component form as Ric, as well. A convention is

86

employed to use the Greek suffixes to evaluate a tensor in the reference configuration and
Latin suffixes in the current configuration.

For the rigid particles suspension, an orientation probability distribution function
N (R, t) is introduced so that N dr gives at time t the fraction of particles whose orientation
states lie within a small region of orientation space dr(E sin 6d6d¢d¢). N (R, t) satisfies

the normalization condition:

f N (R, 0dr = 1 (6.17)
orientation

According to the conservation law in orientation space, the orientation states of particles is

governed by the continuity equation

¥+V-T=O (6.18)

where T is the probability flux vector in orientation space and V is the gradient vector in
that space. It is shown in Rallison’s paper [106], if f is any scalar function of orientation,

then

0f

(Vf) = EkiniabE; (6-19)

The probability flux

7 = N u) (6.20)

in which n) is the angular velocity of particles.

In studying force-free particles suspended in a flow field with no external couples ex-
erted on them, there are two separate contributions on the angular velocity of the particle.
One is the hydrodynamic contribution from the surrounding fluid straining motion and one

is from the Brownian couples [106]. Ignoring the influence of the Brownian motion of the

87

particles, the angular velocity of the particle is
H 1
(02m =S2+-2-B:S (6.21)

where Q = —%e : W is the angular velocity of the‘flow field at the infinity, S is the
strain rate, W is the vorticity tensor of the fluid, and B is a third-order material tensor
introduced by Bretherton [8]. The geometry of such particles (e.g. ellipsoids) and their
interactions with the surrounding fluid are described by the third order tensor B instead of
the single parameter often used for axisymmetric particles (spheroids). By using (6.19) it
can obtained that

V - B = —2F (6.22)

in which

1

F i j = 5 (801.3ka + sklekil) (6.23)

Substituting the angular velocity (6.21) into the continuity equation (6.18), the probability

conservation equation can be obtained as [106]

6N 6N 1 6N
E- + Winjaa—Ric; + EijqEqu‘kiniam " NEijFiJ’ = 0 (6°24)

6.4 Algebraic restrictions on averaged orientation tensors

At an any given point in a particle-filled system, the state of orientation can be described

with a forth and an eighth-order averaged orientation tensors, respectively defined by:

< RiaRJ-fi >= f RiaRJpN(R, t)dr (6.25)
orientation

< RiaRjfiRkleé >= f RiaR jBRkleéN (R, t)dT (6.26)
orientation

88

It is noticed that the average orientation tensor state does not change if switching each pair

of indices in (RiaR 13) and (RiaRjflRkyRM)
(RiaR 1,3) = (Rjfikia) (6.27)
and

(RiaRjflRk'lecS) = (RJBRiaRk'yRkS) = (RkijflRiaR16> = (RléRjBRkyRia)
= (RiaRkyR 13165) = (RiaRkayR 1,3) = (RiaRjBRMRky) (6.28)

Furthermore, since RT - R = R - RT = I, the following relations can be obtained

(R31R31) = 1 - (<R11R11> + (R21R21>) (6-29)
(R32R32> = 1 - ((R12R12) + (R22R22)) (6-30)
(R33R33) = 1 - ((R13R13) + (R23R23>) (6.31)
(R31R32) = - (<R11R12> + (R21R22)) (6-32)
(R31R33> = - (<R11R13) + (R21R23)) (6.33)
(R32R31) = - (<R11R13> + (R22R23)) (6.34)

Using the symmetry and projection conditions, 39 independent entries of a can be listed
by (R11R11). (R11R12>, (R11R1'3). (R11R21>, (R11R22>. (R11R23>. (R11R31), (R11R32).
(R11R33). (R12R12>. (R12R13>. (R12R21). (R12R22>, (R12R23). (R12R31). (R12R32).
(R12R33). (R13R13). (R13R21>. (R13R22). (R13R23). (R13R31). (R13R32). (R13R33).
(R21R21>. (R21R22). (R21R23>. (R21R31). (R21R32). (R21R33). (18221922). (18221823).

(R22R31). (R22R32). (R22R33). (R23R23), (R23R31). (R23R32). and (R23R33). Since
n(RT . R) =.tr(R . RT) = tr(I) = 3 (6.35)

89

the eighth order orientation tensor should also satisfy the following projection properties

(RiaRiaRkyRkS) = (RiaRJBRiaRM) = (RiaRjfiRkyRia)

= (RiaR 13R jfiRlcS) = (RiaRjflRkijfl) = (RiaR'jflRk-yRk-y> = 3 < R > (6.36)

Similarly to the analysis of suspensions of axisymmetric particles, it is convenient from
a computational view to use an alternative approach based on the moments of the proba-
bility distribution function. The derivation of equation (6.24) can be rewritten in terms of
average orientation tensors. An evolution equation for the second moment of the proba-
bility distribution may be obtained by multiplying (6.24) by (RR) and integrating over all
orientations. The time derivative of the fourth order orientation tensor can be expressed
as [106]

a R R
$8771) ‘ Wm (Rme) + (RWRWS) Wm‘

= _%80W83,86<R5#R16RP5R€19>S pq ‘ $80116ng (RthpRpfiRqél) S M (6-37)
(6.37) gives the evolution of (R) in terms of (RRRR). To solve the second moment
equation describing the evolution of the particle orientation (6.37) it is necessary to know
the eighth-order orientation tensor (RR). Closure problems are introduced in (6.37) due to

the unknown eighth-order orientation tensor (RRRR). A fully symmetric quadratic model

for the eighth-order tensor (RRRR) is constructed in the following section.

6.5 Symmetry operator

According to the above discussion the only information available about the structure of the

eighth-order tensor are its symmetry and projection properties. This indicates that higher

90

order tensors can be reconstructed form combinations of lower order tensors and the unit

tensor. It is customary to make the following hypothesis
(RRRR) 2 IF ((RR)) (6.38)

The operator 1F should satisfy the six fold symmetry and six fold projection properties. To
ensure the symmetry of the eighth-order tensor, a symmetry operator S is introduced with

six terms in the form of
S(X, Y) = XY + XYT23 + xvT24 + XY + YXT23 + YXT24 (6.39)

where XYTxy indicates a switch of the xth pair and y‘h pair of indices of the dyadic X, Y
and X and Y are symmetric forth-order tensor ( X = XT12 and Y = YT12 ). Satisfaction of
the six symmetry properties implies that the permutation of any pair of indices of a dyadic
must give the same result. Thus, S(X, Y) is a six-fold symmetry operator. Permuting all six

combinations of indices can prove the symmetry property of S(X, Y) by

S(X, Y) : XiajflkalcS + Xiakijfild + Xial6Yk7 1B + Yiarjflxkyld + Yiakyxjfilé + Yiarldxk'yjfl

I'm-+13
= X jfliakalé + X )pkinazcs + X jflléYk'yia + Y jfliaxkylé + Y jfikyxialé + Y jflldxkyia

iaHk'y
= Xkyijiala + XkyiaY 11316 + XkyzaYiajp + kajflxiarlé + kaiaxjfilé + kazaxia 1'13

('06-’15

: XléjfiYkYia + Xlék'ijflia + Xhfiiakajfl + Yléjfixk'yia + Yldkyxjfiza + Ylézaxkyjfi

jfiHk‘)’
= XiakyY j,616 + XiajBkald + Xialdyjfiky + Yiakyx j,816 + Yiajpxkyza + Yialéxjfiky

JflHM
= Xialékajfl + Xiaksza 113 + XiajfikalcS + Yial6xkyjfi + Yiakyxldjfl + Yiajpxkyza

k'yHld
= XiajBkaza + XiazaYmky + XiakleéjB + xiaffiYk‘r’5 + X‘0’5Yffi"? + X’WY’W

(6.40)

iO‘HjB . . . . . .
where = means to swrtch the rndrces of 1a wrth 1,8.

91

6.6 Construction of the eighth-order orientation tensor

A fully symmetric quadratic closure model is constructed based on the hypothesis of equa-

tion (6.38). Four operators are needed to represent a closed form of (6.38):

P1 = s (8, 8) (6.441)
= 2 (6111 1735/0/16 + (Sta/k7)
P2 = s (8, (RR)) (6.42)
= [<RiaR 1'5) 5ky16 + (RtaRky) 5j,316 + (RiaR18> 5ky 113+
(Sin/1'13 (RkyR18> + 5.61.7 (RJBRRS) + (Sic/16 (Rk'yR 173)]
P3 = S ((RR). (RR)) (6.43)
= 2 (5.8113 (Rth8) + (Siak’y (RjBerS) + 5.618 (RkyR 173))
P4 = s (8, ((RR) : (RR)T12)) (6.44)

= (Sm-,3 (RkyRm><Rm¢R15) + 8,“, (Rf/312m) (RMRM)
4.8% (RkyRm) <Rm¢R 1,3) + (RiaRM) (RMRJ-fl) 61,715
+ (RiaRm) (RmRky)6J-,315 + (RiaRM) (RMR15)6k7j/3
where P1 is a constant term, P2 is a linear term of (RR), and both P3, and P4 are

quadratic terms of (RR). 6 is a 4th order symmetry tensor with unit entries that trans-

poses an arbitrary tensor under double dot product, i.e. 6 : (RRRR) = (RRRR)T12, and

92

(RRRR) : 6 = (RRRR)T34. The components of 6 is shown as the followings:

 

 

H100 010 001
000 000 000
[000 000 000
(000 000 000
6: 100 010 001 (6.45)
(000 000 000
r000 000 000
000 000 000
((100 010 001)

 

The FSQ model is constructed out of a linear combination of two different terms. The
first term contains the constant term and the linear terms of (RR), the second contains

constant term, and the linear and quadratic terms of (RR)

(RRRR>1 = (111131 + 012132 (6.46)

(RRRR)2 = aZlPl + 023133 + (124134 (6.47)

and the closure model is given by

(RRRR) = lF((RRRR)) = (1 — C2) (RRRR)1 + C2 (RRRR)2 (6.48)

(RRRR)1 and (RRRR)2 constructed by (6.46) and (6.47) are symmetric and the coeffi-
cients are selected so that contraction property (6.36) should be satisfied. After contraction
of the linear closure and the quadratic closure that satisfy projection and symmetry proper-

ties the following coefficients can be obtained

0'1] = —§§6 (6.49)

93

3

012 = E (6.50)
3(sz — 112)
._. ____ . 1
“2‘ 14311 (6 5 )
a - -i (6 52)
6

in which J], and 12 are the first two invariants of the fourth-order tensor (RR) [110] with
the following expressions
J 1 = -tr((RR)) (6-54)
and
J2 = -% ((RiaRia> (ij1813} - (1%]?th (Rm¢Ria>) (6-55)
Substituting these coefficients into (6.46) and (6.47), the two terms of the closure model

can be rewritten as

(RRRR)1 = —2—86S(6, 8) +— 33,S(6 a) (6.56)
3(an 42) 3 6
(RRRR)2 = ———14—3jl—s (8, 8) — —s (a, a) + HS (8, (RR)) (6.57)

Substituting (6.56) and (6.57) into (6.48), the closure model for the fourth moment orien-

tation tensor can be obtained as

< RRRR>- _ (1 —C2)[— 2—86S(6, 8)+ 336(6, 3) (6.58)
3(212-12) _3_2 6

 

which preserve both the six-fold symmetry and six-fold projection properties.
For suspensions of spheroids, to satisfy the following suflicient condition for realizabil-

ity of the orientation dyadic (i.e., microstructure), i.e.,

z-<pp>-220 (6.59)

94

C 2 has been shown to be

_ 8 + 45111.,
2 ' 18(1 + 91111,)

 

(6.60)

in which b = < pp > —- %6 is the anisotropic orientation tensor for spheroids and IIIb =
tr(b-b-b) the third invariant of b [11 1,112]. For suspensions of non-axisymmetric particles,
C 2 is still in development. The tensor calculation associated with this chapter can be found

in Appendix A.

6.7 Conclusions

As suggested by Rallison [106], it may be practical to predict the microstructure of a sus-
pension of rigid , non-axisymmetric particles by using the rotation operator as a state vari-
able rather than the Euler angles. This research has identified a closure for the 4th-order
moment of the orientation distribution function in terms of the 2nd—order moment that sat-
isfies all six-fold symmetry and projection properties of the exact 4th-order moment. This
result may provide a means to improve the accuracy of estimating the rotary diffusion co-

efficient from retum-to-isotropy experiments.

95

CHAPTER 7

PREDICTION OF FLOW-INDUCED

ORIENTATION AND SPATIAL

MIGRATION OF PARTICLES

7 .1 Introduction

Particle migrations in suspension flows are of importance to a variety of scientific and en-
gineering applications, e. g. the transport of sediments, chromotography, composite materi-
als processing, sequestration processes in porous media, and secondary oil recovery tech-
niques. For suspensions with micron size particles, for which inertia effects and Brownian
diffusion can be neglected, the interaction between the particles adds a random compo-
nent to their motion that is additional to the deterministic translation along streaming in
the slow viscous environment. This random component results in migration of particles,

which was first identified by Leighton [113]. Some valuable information on the many-

96

particle interactions have been simulated by Stokesian Dynamics and boundary element
methods [72, 114—118]. The experiments of Segré and Silberberg [119, 120] have large
influence on fluid mechanics studies of migration and lift of particles. They studied the
migration of dilute suspensions of neutrally buoyant 'spheres in a pipe flow at Reynolds
numbers between 2 and 700. The particles migrate from the wall and centerline and ac-
cumulate at about 0.6 of a pipe radius from the centerline. Karnis et al., verified the same
phenomenon and observed that particles migrate faster for larger flow rate and closer to the
axis for the larger rigid sphere. Aubert et a1. [121,122] found that there were some forms
of migration in all flows, curved or uncurved; however, in parallel flows he found that there
could be no rrrigration perpendicular to the direction of flow; that is, the polymer just lags
or precedes the flow along a single streamline. He found further that cross-streamline mi-
gration occurs in curvilinear (e. g., circular Couette) flow when he approximated the curved
flow as a quadratic. Flow-induced polymer migration is treated in a rather intuitive though
mathematically formal way by Sekhon et al. to study the effects due to the hydrodynamic
interaction as well as flow geometry [123].

Two different approaches have been used to develop models capable of describing
various multiphase flow regimes. The first case is known as the Dilute phase approach,
also called the Lagrangian approach, in which the fluid phase is treated as a contin-
uum and the particle trajectories are calculated for the equation of particle motions. La-
grangian method is used in modelling the dynamics of a single particle or a dilute sus-
pension [10, 14, 124, 125]. The second approach is known as the Dense phase approach,
sometimes also called the Eulerian (or two-fluid) approach. In this approach, each phase (or

component) directly influences the motion and the behavior of the other phase and the par—

97

ticle phase also is treated as a continuum. This method is used widely in fluidization [126],
gas-solid flows [127], pneumatic conveying [128], and suspensions [129].

Two continuum theories are developed in the dense phase approach: Mixture theory and
Averaging theory [130]. Both theories are based on the assumption that each phase may be
mathematically described as a continuum. The ideas of Mixture theory can be traced back
to the branch of mechanics [131—134]. The fundamental assumption in Mixture theory
is that at any instant time, all phases are present at every material point. In contrast, the
Averaging method directly modifies the classical transport equations to account for the
discontinuities for ’ jump’ conditions at moving boundaries between the phases [135,136].
The modified balance equations must then be averaged in either space, time or statistical
to arrive at an acceptable local form [23, 137, 138]. Some difference of the equations of
two-fluid by Mixture theory and by Ensemble Average theory can be found in [134].

Three essential parts are composed of the formulation in the Eulerian approach: the
derivation of field equations, constitutive equations, and interfacial conditions. The field
equations state the conservation principles for, e.g. the mass, momentum, and energy.
Constitutive equations close the equation system by taking into account the structure of
the flow field and material properties by experiment correlations. It’s noted that both of the
Mixture approach and the Averaging approach are not closed and methods of closure, or the
constitutive equations for the interaction terms, are required to put the equations of motion
into a form suitable for application. Because of its close relation to measuring techniques,

the Averaging method is most widely used in the multiphase flows.

98

7 .2 Hydrodynamics of ensembles of particles

7 .2.1 Theory of ensemble averaging

As mentioned in the previous section, there are several kinds of averaging methods, i.e.,
time average, volume average, and ensemble average can be applied in the average ap-
proach to solve for the multiphase flow. Comparing with the other averaging methods,
ensemble averaging method has some advantages and is widely used in the current analysis
of multiphase flow [23,134, 135, 139, 140]. First, the data acquired by time and/or volume
averaging can be easily used as the ”sample” of the ensemble. Second, ensemble averaging
does not require that a control volume contain a large number of particles in any given
realization. Third, ensemble averaging is easily implemented. Forth, the ensemble average
allows for that all realizations are only approximations of the ideal. Detail information on
the ensemble average is given by Drew and Passman [135]. The definition of ensemble

average is
f(x, t) = ff(x, t;,u)dm(y) (7.1)
E
where dm(-) is the density for the measure (probability) on the set of all precesses 6. Some

results can be applied to the ensemble average in order to average the equations of motion:

(1) Reynolds rules

Qfi+Qh=mfi+QZ 03

(2) Treating generalized function

f. s

 

 

)vd dt- _ -6tf 6“" t)f(x ,t)dvdt (7.3)
Q

99

and
f ¢(x,t)Vf(x,t)dvdt = — f V¢(x, t) f(x, t)dvdt (7.4)
Q Q

in which 05 is a test function belongs to (I), Q ,a compact set in space and time to

support (I) E (I).

(3) Interface Delta function and Topological Equation

an
— = - VX -
6n nk k (7 5)

This is the interface Delta function where Xk is the characteristic function as

1 if phase k occupies x
Xk(x, t) = (7-6)

0 otherwise
The Topological equation is

6X k

 

 

‘5;- + Vi ' VXk = 0 (7.7)
(4) Gauss and Leibniz rules
Xka = ‘7ka - fVXk
= V(ka) - fkiVXk (7.8)

This is called the Gauss rule in which fki is the value of the function f evaluated on

the component k side of the interface.

Similarly, the Leibnitz rule is

X (if — _anf _ 4939‘.
at _ at at
__ 3m 3X1.

100

7 .2.2 Averaged balanced equations

A two-fluid model is used in the balanced of mass, and momentum equations can be ob-
tained by taking the product of the balanced equations with the phase indicator, Xk, ma-
nipulating using the product rule for differentiation, and then performing the averaging

process. Mass balanced equation can be written as

 

ka —— 6Xk
—— V-X = — -VX .
6t + kpv p( at +v k) (710)

By using the topological equation 7.7, the right-hand side can be reduced to

 

Fk = [.0(V - Vi)] ' VXk (7.11)

This is the interfacial source of mass due to the phase change. If (v — vi) - n = 0, then

I‘k = O. The averaged density and averaged velocity of phase k can be defined by

akin. = m (7.12)

 

akfikvk = kav x (7.13)

Substituting (7.11 - 7.13) into (7.10), the averaged balanced of mass for the phase k can be

obtained by

001.51.
at

 

+ v . 8,8ka = r, (7.14)
Multiplying the equation of balance of momentum by Xk and taking average to it, the

averaged momentum equation for phase k can be obtained

6va
at

 

 

 

 

+ V - kavv = v . XkT + kag + pv[(v - v,-) . vxk] — T - VXk (7.15)

in which T is the stress tensor, and g is the body force, e. g., gravity and magnesium. Defin-

ing the averaged stress by

0ka = 3031". (7.16)

101

the Reynolds stress by

 

ak'rfe = —W, (7.17)
the interfacial velocity by
v73. = pv(v — vi) - VXk (7.18)
and the interfacial force by
Mk = 5773?; (7.19)

finally the equation of balance of momentum can be reduced to

(90/14.)ka

at + v - 8,,kaka = v 6km + T?) + akpkg + Mk + v2.1“), (7 .20)

7 .3 Equations of motion and orientation for a dilute sus-

pension

To distinguish each constituent of two-phase flows, subscripts ”f” and ”s” are used to rep—
resent the continuous phase (fluid) and the solid phase (particle) respectively. If the con-

centration of the solid phase is a, the concentration of the fluid phase will be
ale-as=l—ar (7.21)

The ”—” over the variables means ensemble averaged quantities. For convenience, the
overline may be taken off to yield simpler expressions. Assuming I‘d/f = 0, the conversa-

tion equations for mass of the constituents are

 

6t!
—6t +V-arvs = 0 (7-22)
6(1—0)
at +V-(1—a)vf = 0 (7.23)

102

The corresponding equations of momentum for the dilute suspension are of the form

Ups (9'5"; + V3 ' VVS) = V ' GT5 + apsg + MS (7.24)
an
(1—a)pf(E-+vf-va)=V-(l—a)Tf+apfg—Ms (7.25)

In the above, a is the concentration of the solid phase. p s, and p f are the partial densities
of the solid and fluid phase respectively; the mass weighted averaged of solid and fluid
velocity are V, and vf respectively; the phase interaction force per unit volume is denoted
by M s and the mass-weighted stress for the solid and fluid phase are T S and T s respectively.

For dilute suspensions, the interaction between the particles is negligible. According to

(6.10) the orientation of the particles can be written as

g.=—W-a—a-WT+4(S-a+a-S—ZS:A) (7.26)

7 .4 Stress model

It is noted that the averaged balanced equations are not closed. In order to make the above
equations solvable, constitutive relations must be obtained for the phase interaction force
and the phase stresses. There are many distinctly different modes. Hwang and Shen [141,
142] provided a derivation of the solid phase stress using a a control volume/control surface
approach. Consider a rrrixture of an incompressible Newtonian fluid and rigid particles with

a uniform size. The fluid and solid phase stress may be decomposed, respectively, as
Tf = -pr+T’f (7.27)

and

T, = —p,1 + T; (7.28)

103

where — p f and — p s are the phase pressures, I is the unit tensor, and T; and T} the devi-
atoric parts of the phase stresses. The deviatoric parts of the phases can be decomposed,
respectively, as

T’f = T, + T, ' (7.29)

and

in which TV is the fluid viscous stress, T, the fluid turbulence stress (or Reynolds stress), Tc
the collision/contact stress, Tk the kinetic stress (equivalent to the solid turbulence stress)
and Tp the particle-presence stress resulting from the hydrodynamic forces acting on the
particles.

Consider dilute suspensions of rigid particles in an incompressible Newtonian flows at

a low Reynolds number. Due to the low Reynolds number, the fluid is laminar. Hence
T, = 0 ' (7.31)

For dilute suspensions, the concentration of the particle approaches zero so that there is no
particle collisions, consequently

TC = 0 (7.32)

Furthermore, driven by a laminar fluid motion and a stationary body force, particles do not

fluctuate. In the absence of particle collision induced random motion, this implies
Tk = 0 (7.33)
For the Newtonian fluid considered here there shear stress is well defined as
T, = gm]. + (va)T] (7.34)

104

The physical interpretation for T p was first given by Batchelor [21,143] using a volume

averaging concept. The resulting stress has the following form

1 f
T " = — Z-knkr'dA- f 6 2' r'dV] (7.35)
pl] VO[ A0 I J V0 ’6 1k ]

where V0 and A0 are the volume and the surface of a single particle, zik the hydrodynam-
ically induced local stress at dA on the surface of a particle or at dV in side a particle, nk
is the kth component of a unit outward normal on the particle’s surface and r J- is the jth
component of the position vector of the infinitesimals dA or dV. The averaged particle

stress for uniform suspensions of ellipsoids with neglecting the inertia force is written as
Tp = kA : S (7.36)

in which A is the forth order orientation tensor of the particles, k a factor depending on the

fiber length and the fiber concentration. The expression for k is

7! vL3
k— ”f

— mflé‘) (7-37)

 

, and e = [1n 2L/d]‘1. The

1
where vL3 the volume fraction of the particles, f (8) = 1 1 58

concentration ratio of the particle in suspension is
_ 2 _ 3 2
a — 7rd Lv/4 — 7rvL /(4ap) (7.38)

in which 0,, = L/d is the aspect ratio of the particle. It can be seen that the factor k depends
more on the concentration ratio of the particles [107].
An alternative formulation of the solid phase stress is given by Hwang and Shen [141]

based on the concept of utilizing a control surface and considering stress as the force per

105

unit area on such a surface. The resulting formulation of this stress is identical to Batche-
lor’s [21]. Because the stress components Tk and Tc are strictly modelled from the momen-
tum transfer rate across a control surface, it is different from volume-averaging of internal
solid stress. In order to be consistent with the concept of deriving Tk and Tc, a control
volume/control surface approach is adopted to derive the particle-presence stress Tp by

Hwang and Shen [141]. The form of particle-presence stress is
n
T ~ = — 2- -dA— 62' -dV 7.
PU “(j/40 lknkrj 170 k rkr} ) (39)

I
: —[f ziknkrjdA—f akzikrjdV] (7-40)
V0 A0 V0

which is identical to Batchelor’s [21] result for a slow flow of a dilute uniform fluid-solid

mixture. For this flow, the particle-presence pressure is the form

10:12 370 —f02” f2” psin¢d¢d0 (7.41)

7 .5 Interfacial force

In [142], the derivation of the phase interaction term M 3 is provided based on the same
concept of control volume/control surface approach used in deriving Tp. The form of M S

is

(I

where h is the hydrodynamic force, acting on a single particle, V0 the volume of a single
particle. For a dilute mixture, his approximated by hydrodynamic forces of a single particle

in an infinite fluid flow. With the additional assumption of low Reynolds number of the

106

particle and fluid, h is composed of several contribution as the follows:

in which f5 is the Stokes drag acting on a particle, fa the additional forces including the
added mass effect, the Basset force [75] and the Saffman [144] force due to the fluid inertia.

Substituting (7.43) into (7.42), the phase interaction is given by
a

The analytical hydrodynamic force on a single particle suspending in unbounded creep-
ing flows with a constant velocity gradient is derived in Chapter 2. Induced by the linear
shear flow the particle may rotate and translate inside the flow. The hydrodynamic force

described in the rotating coordinate system is
f; = K’(v} — V's) (7.45)

where v} and v; are the velocities of the surrounding fluid and particle respectively, K’ is
a resistance tensor. If a spheroid suspending in linear shear flow with no slip boundary at

the interface, K’ is given as

6 —5 + 28 0 0
K’ = —§7rr0p 0 —5 + 28 0 (7-46)
0 0 —5 + 8

where 3 is the deformation coeflicient defined in Chapter 2. The hydrodynamic force de-

scribed in the fixed coordinate system has the form

fr

K(vf — vs) (7-47)

= RT - K’ -R(vf - vs) (7.48)

107

in which RT is the transformation matrix between the rotating coordinate system and the

fixed coordinate system, defined in (2.78). It is known that the orientation of the particle is

defined by
sin 6 sin (b
p = - sin 6 cos (15 (7-49)
cos 6
Hence, the dyadic of pp is
a = PP
sin2 6 sin2 (b — cos (1) sin2 6 sin (13 cos 6 sin 6 sin 45
= - cos (1) sin2 6 sin (12 cos2 ¢ sin2 6 — cos 6 cos ¢ sin 6 (7.50)
cos 6 sin 6 sin (15 — cos 6 cos ()5 sin 6 cos2 6

Substituting (2.78) into (7.47) and comparing with (7.50), then
6
K = —§7rro/J [(-5 + 28)] — 8a] (7.51)

If the influence of addition force acting on spheroids is not considered for a general transient

flow, the phase interfacial force can be written as the follows:
a

where V0 = gurgfl - 8) for a spheroid. From (7.52), it can been seen that the interaction
force on the particles depends on the velocity difference between the solid phase and the
fluid phase, the averaged orientation of the particles, the viscous stress term, and the particle
stress term. As discussed in Chapter 2, the variable 8 in (7.51) is the deformation of a
particle from a sphere. The K matrix is accurate to 0(82). Hence, the model for the

interfacial force on the particles shown by (7.52) is valid for a small number of 8.

108

7 .6 Summery

A numerical model is developed to describe solid-fluid two phase flows using a continuum
approach. A so-called Eulerian-Eulerian technique is adopted to deal with the motion of
the spherical particles and Newtonian fluid. Based on the moments of the distribution
function, the evolution of the second moment of the orientation tensor is used to govern the
orientation of particles statistically. The concept of control volume/control surface method
is used to develop closure models for the stresses and interfacial force on the particles. The

model for the interfacial force is valid for small deformation of the particles from spheres.

109

CHAPTER 8

SIMULATION OF THE FLOW

INDUCED FIBER ORIENTATION

AND MIGRATION USING A FINITE

ELEMENT METHOD

8.1 Governing equations for 2-dimension problems

8.1.1 Basic assumptions

In this chapter, a simple 2-dimensional problem is investigated by using the finite element
method to solve the governing equations. According to the governing equations of the
dilute suspension system introduced in Chapter 7, some further assumptions are taken into

account to build the equations of motion for the 2-dimensional problems.

110

Velocities of the fluid and the particles in the x direction are assumed to be zeros. Par-
ticles rotates in the y-z plane induced by the surrounding flow fields. Due to their weak
contributions the components of the orientation state a, i. g. a xy and a xz and the compo-
nents of the particle stress term Tpxx, Tpxy, and Tpxz (are assumed to be zeros. Even though
the normal stress component T p x x is not exact zeros for the FSQ closure model, it is negli-
gible compared with the effect of the magnitudes of the other components, i.g. pry, prz,
and szz. According to the normalization condition of the orientation of a particle, the
component of a xx can be determined by a xx = 1 — ayy — an. Therefore for 2-dimensional
problems, both the velocity field of fluids and particles have two independent entries v fy,
vfz and vsy, vsz; due to the symmetric property ayz = azy and the normalization condition
a xx = 1 — ayy — azz, the average orientation tensor a has three independent entries, i.g.
ayy, azz, and ayz; only three independent entries pry, szz, and prz are evaluated for the
particle stress term Tp. Including the concentration ratio of the particles and the pressure of
the fluid, there are totally 12 unknowns for the 2-dimensional suspension system, namely
Vsy, Vsz’ ny’ sz’ Tim” szz’ prz’ “yr, azz’ “th a, and Pf:

The fourth order orientation tensor A in (7.26) and (7.36) requires a closure model for
it. The fully symmetry model [102, 103, 145] developed at Michigan State University is
applied in the following simulations. This closure model retains all the six symmetry and
projection properties of the fourth order tensor. Different forms of this model have been
discussed in the work of Manda] [145]. Herein, only linear and quadratic terms in the fully
symmetric model are considered (the so-called FSQ model) in implementing the equations.
The coefficient C2 in the FSQ model depends on the third invariant 111;, of the anisotropic

part of the averaged orientation tensor [111,112]. It has been noticed that when C 2 = 0.37,

111

the closure model can provide pretty good results for simple shear flows as well as other
different flow field, i.g. uniaxial flow [145]. For simplicity a constant value of C 2 = 0.37 is
selected in the FSQ model.

The computed microstructure can be interpreted with the help of eigenvalues and eigen-
vectors of a. The eigenvalues of a has the property of 1 / 3 S /l S 1.0, in which .lmax is
the maximum eigenvalue in the domain. Furthermore, if .1 —- 1 (i.e. Am)" 2 0.0), fibers are
aligned in the direction of corresponding eigen vector. On the other hand, if 11mm, 2 1 / 3
(i.e. Am)" 2 1/3), fibers are orientated randomly in all directions. Setting emax as the
normalized eigen vector associated with the maximum eigenvalue, the nricrostructure has

1

2
been interpreted by plotting (4mm — 3) / 5. This will result in a vector of zero length for a

random orientation state and unit length for a uniaxial alignment state.

8.1.2 Governing equations

The simplified governing equations for the migration and alignment of ellipsoidal particles

are presented below.

(1) Continuity equation of solids

Vector form:

0’
E+V'(O’VS)=O (81)

Component form:

 

 

a
aa+a( Vsy+avsz) ( 6a Ba):0 (8.2)

5 6y 6z v9.5; + VHS—z-
(2) Continuity equation of fluids

112

Vector form:

6(1—a) +

 

6:

Component form:

6t 6y

(3) Momentum equations of solids

Vector form:

6v
—§9/-+(1—ar)(—fy

v.[(1— a)vf] = 0 (8.3)

(9sz 60 60
+ 7)— (V1779; + W232?) — 0 (8.4)

6v 0
psa (ii—ts + v5 - Vvs) — V—OK(vf — vs) — aV - {—pfl + 2prS] —psag = 0 (8.5)

9apf

(91/3,- (9V3,-
psa(— (9t +V Vsj axj ——'+) [m] [(—5 + 28)5ij - 861,7] (ij - vsj) (8.6)

6v,

6va

’a—[_ pf‘Sij+ ”f(b—fj‘l'a—xn— psagi=0
l

y-Component form:

 

6vsy at.”
p50 (7! +1030 Vsy— 0y +V Vsz az

[ 9apf

—)

(8.7)

90;:
+ m] [(-—5 + 28) — Sayy] (ny "' Vsy) + ("W—1:83] {—anzj (va — v52)

avfy avfy
”El—”’1‘ ”(6 ““67

z-Component form:

6vsz 6vsz
psa’ TI +950 VsyE'+Vsz—62:—

9apf 9011f
[10%(1- t:.~)][_8a”z](vfy — v”) +[10r2(1 —

-a-[ (sf—mi- :
062:

62

)

6v 6v
_a_3_[# ( fy fz
6z

‘6.— “(97117483 = 0

(8.8)

8)][(-5 + 28) — 8azz] (sz — vsz)

6v 6v
‘Pf +—+#f( f+ —f—Z)J-pragz=0

62 6z

(4) Momentum equations of fluids
Vector form:
6v f
pf(1-a)(-—a-t—+Vf-VVf)-(l -Za)V-[—pfl+2flVS] (8.9)

+ {—pr + ZpVS] - Va + gamvf — v,) — v - (an) —pf(1 - a)g = 0

an _if 6Vf, 6v_f,- (9ij
pf(1- a)— +v ij— 6x] --(1 20)£— -pf6ij +,uf a—-xj + a—x (8.10)
r

(”(3% (3311') 6a [x 90w

+ _ _ _—
5x j 6x: 10r(2)(1 — 8)

[ea-m. ]1<>1<>

6 6a
-05x—ijij - TpijaTj "pf(1 - a)8.° = 0

y-Component form:

 

6vf 6vf
paw-a)( fy—)+pf(l (l)(vfy ayy+VfZ azy) (8.11)
6va +6va _ a"fy +6va
“1 za’gl’pf ”flay _ll‘(1’2")6_z ”(62 6)]
' 6v 6v 6v 6v 6
Wt 412—: 1 (6+2: .4) .—:

9 9
_ LJR— —5—+28) 80M] _vs_y) [10,-2:1 f ][—80)z] (sz Vsz)
a 6 6 60

(106(1— e)
a
"9 5T1)” (Ya—szrz pryg'y‘ — przaj ‘Pf (1 ‘ “My = O

 

 

 

z-Component form:

avf —Z)avf (8.12)

pf(l—a)(: __f_:;Z)+pf(l—a)(vfyfz 6y +sz 6z
6"_fz +_6V_z_fy)]_ (1_

 

av—fz 6sz

-(1—2a)— ,uf— Za)— —pf+ pf 6—z + 6z
anz 6Vfayy 60' avf: +avfz (9C!

[+"f(a_y+az)a_y +"[Pflf‘fLi—zz‘i 6_z)j6z

_ 9041f _ _ _ 9apf _ _ _
[10%(1_8)][ “921(ny V5?) [10%(1_8)][( 5+28) 8azz](sz V32)

6 6 6a 6a

114

(5) Particle stress

Tensor form:
TP 2 kA : S (8.13)
yy-Component form:

1 2 6v

62

 

—%[—1 +C2 +2C2(1—ayy —azz)

“W F ‘9va 6va
-3; £10 (30 - 65C2 + 94C2ayy) 797 + (5 — 4OC2 + 35C2ayy) 797
_%z

35

 

_ avfy
(5 5C2 + 35cm),y + 4C2azz) :97-
6V
:1.
6z

_3_—aw fz
3—5— (5— 5C2+2502ay.—, locza:)(aav +—)

5)’
6ny+ 6sz

zz-Component form:

1 2 fy
szz —7[—1 +C2 +2C2(1—ayy—azz) — (8.15)
1 2 anz

 

(1ny all—Lu 312
a 6v
"3% L5 (1 - 8C2 + 14cm), + 7C2au) 65—31 (8.16)
6va
+ (30 — 65C2 + 35cmy + 94C2azz) Tz—
_3_f__yz (9va 6va
35 (5— 5C2 -10C2ayy + ZSCzaZZ)(— 62 + 6y —)
3% 3" _"_f2
= O

115

yz-Component form:

1 ,

TPYZ

"35
_ 1
35

1%

=0

35.
1.
_351
1,

 

 

 

 

_1+ C2 + 2C2 (1 — a)... — azz)2 62—?

_ —1+C2 +2C2(1—ayy—azz)2 9%

_ (5 — 5C2 — 8C2ayy)(a:;:y + 6:52)]11”
(5 502 + 35c2ayy - 802%) (6%? + 6:52)](13
(2— 9C2 + 24C2ayy + 3cm.) (6%) :4ng ayz(_

1 6v fz
+77- (2 - 9C2 + 3C20yy + 24C2aZZ) (WH ayz

(6) Orientation equations

Vector form:

6a

a:

+vf-Va+W-a+a-W—A(S'a+a-S—ZS:a)=0

yy-Component form:

+9}
35

T5

——4

35

2,1 (5 — 5C2 + 35cm), + 4c2azz)
+

f

 

(8.17)
‘9ny + 19.11:)
62 6y
(8.18)
(8.19)

(—5 — 65C2 + 94C2ayy) £39 + (5 — 8C2 + 7C2ayy) % lay).

6ny
72%—

+ (10 — 10C2 + 70cm), — 10C2azz)

ém _ 10%.)

6va “a
‘82—

—5 [7 + ,1 (—1 — 6C2 + 30C2ayy — 12C2azz)] 6%.?

1 +5 [7 + 11(1+ 6C2 - 30C2ayy +12C2“zz)]§v3¥

116

1ayz

 

zz-Component form:

aazz 661?; 6082
8.2
_at +(ny—6y +VfZ az ) ( 0)
2/1 , 2 . 6ny 6sz_
.1.3 h_1+C2+2C2(1—ayy—azz) ]( 6y +7 62)

g; ’ (—10 +10C2 + 1062a”) 6% + (—5 + 5C2 - 4C20yy) avg? 1“»

F'

6v
+2 (5 — 8C2 + 14C2ayy + 7C2azz) 46; a
35 6v fz zz
( + (—5 — 65C; + 35C2ayy + 94C2azz) Tz

r 6v av
lelCzayz (—10?§Z + 773%)
1

‘3—5* +5 {—7 + 1(1 + 6C2 + 12C2ayy - 30C20zz)] :31 my

 

 

 

1 +5 [7 + 1(1 + 6C2 +12C2ayy — 30C20zz)] 6%“)? 1

= 0
yz-Component form:

6ayz ( 6ayz (MW) (8 21)

797- + vfy—ay- + VfZ—a—Z— .
2,1 2 6ny 6va
+§§[—1+C2+2C2(1—ayy—azz) 1(32—4'797
, a
1 [_35 + ,1 (15 + 20C2 + 3202a”) g?-
75‘ + [35 + 1(15 + 20C2 + 32C2a ) 6va a”
1 y’ 797

6
_i I [_35 + 1(15 + 20C2 + —140C2ayy + 32C20zz)] 7:)? a

22

 

 

 

a
70 1 + [35 + ,1 (15 + 20C2 - 140C2ayy + 32C2azz) -§¥
r 6ny 6va 5 ‘
31 3663“” (7&— + 76—y
a
+3? + (—5 - 3OC2 + 80C2ayy + 10026122972;fl ”W
6sz
( + (—5 — 30C2 + loczayy + 80C2azz) 33- J

 

117

8.1.3 Boundary conditions

The boundary conditions imposed on the geometry are

(1) Dirichlet boundary conditions:

Vf/S=fu on 1‘“ (8.22)

(2) Neumann boundary conditions:
t = (—pI + [1 f [(va) + (vafl) - n on 1“,. (8.23)

where n is the unit normal to the boundary and 1“,, and I", are Dirichlet boundary and

Neumann boundary and shown in Figure 8.2 [146].

8.2 Mixed finite element model

8.2.1 Weak form

The finite element method is used to solve this problem numerically. The starting point
to develop the finite element models of (8.3)—(8.21) is their weak statements. The weak
forms of (8.3)-(8.21) over an element (2" can be obtained by a three-step procedure. These
steps are briefly reviewed here. First we multiply the differential equations with different
wight functions, and integrate over the element. To distribute differentiation equally among
all variables such that the finite element approximation functions satisfy the continuity
requirement it is necessary to take integration by parts in the second step. The third step
consists in expressing the boundary integral terms as functions of known quantities. The

weak form development is shown below.

118

(1) Weak form of the momentum equation of solids

6v__s_,- va—VSi
90 [1f .
+ L} st {[m] [(:5 + 28)6ij - Eaij] (ij — vsj)_ psagi}dQ

6sta 3W ‘9ij f _
+L{ axj [pf6ij+[1f(axj +—67i- d0 erat,dl"—O

(2) Weak form of the momentum equation of fluids

 

w 1 avf’ avf‘ do 825
1;; vf{pf( 'a)(— at +vaa—xj')} ( )
+1.1 [mm +r.)1}dQ-f..w<l-2a>ridr

P (9Vfi +6ij 60;
+vaf{ -pf5ij+/Jf(-a-—xj+ an —)J]a—xj}dfl

( 9apf
_. Qwvf{km][(—5+28)6ij—Saij] (vfj—vsj)}d§2

6 (90
— vaf {aaTij-j - Tpij‘aTj ”Pf“ — (1)8i}dQ = 0

 

 

(3) Weak form of the particle stress
f pr {Tp — kA : S}dQ = o (8.26)
Q
(4) Weak form of the average orientation state

fWa{%+W-a+a-WT-/l(S-a+a-S—2S:A)}dQ=O (8.27)
(2

(5) Weak form of the continuity equation of solids

f Wa ((9—0 + V - avs)d§2 = 0 (8.28)
Q at

119

 

 

 

 

 

 

U)
U)

 

 

 

 

 

(a) (b)

Figure 8.1. Quadrilateral elements used for the finite element model. (a) A nine-node bi-
quadratic element is used for the shape function of velocities. (b) A four-node continuous-

bilinear element is used for the shape function of the pressure of fluids.

(6) Weak form of the continuity equation of fluids

L—wpf[i(%93+v-(1-a)vf dQ=0 (8.29)

where W are the weight functions and the superscripts vs, vf, Tp, a, a, and pf denotes the
weighting functions for the velocity of solids, velocity of fluids, particle stress, orientation

tensor, concentration ratio, and pressure of fluids respectively.

8.2.2 Finite element model

Since Galerkin method is applied in finite element models, the same interpolation functions
as the weight functions (isoparametric) are used to approximate the dependent variables vs),

vfi, Tp,-j, aij, a, and p f. Suppose the dependent variables are approximated by expansions

of the form
”I . .
v,,- z 2 wgsvfw. = wfsvsi (8.30)
i=1

120

vf, z 2 ijvk. = wffvf, (8.31)
i=1
”3
~ ,- ,- _ T ..
Tp ,j ~ 2 WTPTPU. ‘. WTTPU (8.32)
j=1
N4 .
aij x Z Wgaf = wZaij (8.33)
i=1
N5 .
a x Z Wéaj = wl‘a, (8.34)
i=1
N6 .
Pf 3 Z “’5pr = ngfi (835)
i=1

Lagrangian type of polynomials are used for the interpolation functions. In order to
prevent an overconstrained system of discrete equations, the interpolation functions for
pressure should be at least one order lower than that used to velocities field to satisfy
the LBB (Ladyzehskaya, Babuske, Brezzi) conditions [147]. For two-dimensional flows
nine-node rectangular element shown in Figure 8.1(a). The velocity component, and other
variables, i. g. particle stress tensor, and the orientation state tensor are approximated by bi-
quadratic Lagrangian functions. These functions are expressed in terms of the normalized

coordinates s, t for the element, which vary from —1 to 1, given as the following

( W

(s2 - sxr2 - t)
(s?- + s)(t2 - z)
(s2 + s)(t2 + z)
(s2 — s)(t2 + t)
4 2(1 — s2)(:2 — t) ) (8.36)
2(s2 + s)(1 — 9)
2(1 - 32x:2 + :)
2(s2 — s)(l — :2)
2(1- 52)(1— :2) J

 

 

121

A 4-node continuous-bilinear element shown in Figure 8.1(b) is used to approximate the
pressure of fluids. The bilinear interpolation functions are defined as

'(1-s)(1—t)l
(1+s)(1'—z)
( (1+s)(1+t)
(32-S)(t2+t)
((l-s)(1+t))

Substituting (8.30-8.35) to (8.24-8.29), we can get the matrix form of the weak from of the

V

(8.37)

#I"

 

 

governing equations

Vsz Vsz
ny ny
sz sz
pry pry
6 T T
M5- P“ +K P” =F (8.38)
I prz prz
ayy “yy
azz azz
ayz ayz
(I (I
1 Pf J 1 Pf J

 

 

 

 

in which M is the mass matrix, K the stiffness matrix, and F the force vector. Their explicit
forms are shown in Appendix A. Eq. (8.38) can be rewritten into a more symbolic format

as

MU + KU = F (8.39)

where

T
U = ( v5); vSZ ny sz pry szz prz ayy azz ayz a pf ) (8'40)

122

In (8.40) there are 12 unknowns and 12 equations. Suitable boundary conditions and initial
conditions are needed to solve this equation. The general form of (8.40) is nonlinear and
time-dependent. Therefore, the first order backward difference scheme is applied with a
relaxation factor of 0.5 and a Picard iterative method is adopted to obtain the solution. A
finite element code is developed to predict the flow induced orientation and migration of

suspensions in complex geometry by using Matlab7.0.

8.3 Simulation of a plane Poiseuille flow

Consider the slow flow with particle suspension between two long parallel plates at rest
shown in Figure 8.2(a). This flow is driven by a pressure gradient in the axial direction.
This kind of flow is often called a plane Poiseuille flow. When the length of the plate is
very large compared to both the width and the distance between the plates, it is a case of
a plane flow. 2H and 2L denote, respectively, the distance between and the length of the
plates [see Figure 8.2]. At the inlet both of the velocity profiles for the fluids and particle
V0 have been specified as a parabolic function of z. Particles are ejected randomly at the
inlet with the constant concentration ratio of 0.01 (a semidilute concentration if L/d = 50).
The parameters associated with glycerin for the material of fluids and sand for the material
of particles have been used in this problem, i.e. density of fluids p f = 126Okg/m3, density
of solids p s = 2500kg/m3, dynamic viscosity of the fluid [if = 1.5Ns/m2. Due to the
axial-symmetry in this problem, it suffices to model only half of the domain. BL, BR, BB,
and BU shown in Figure 8.2(b) represent the left, right, bottom, and the upper boundary of

the half domain respectively. The boundary conditions for the half domain are set as the

123

Computational domain

 

fl:

 

 

III/IIIIIIIIIIIIII [III/111111111]

 

L
(a)

 

 

 

Quadratic element Linear element

 

(b)

Figure 8.2. Domain and mesh for a plane Poiseuille flow with particle suspensions. (a)
Geometry, computational domain, and (b) the finite element mesh used for the analysis of

the slow flow with particle suspensions between parallel plates.

124

following:

AtBL:
WU
sz
1/3 0
a: O 1/3
0 O
AtBR:
AtBU:
AtBB:

Uzpf
6v

0
0 (random orientation)
1 / 3
a=0m

 

any = (9:5),— :0
6y 6y
ny—V5y=0

vfzzvsz =0

6W2 &w

 

82 dz

125

(8.41)
(8.42)
(8.43)

(8.44)

(8.45)

(8.46)

(8.47)

(8.48)

(8.49)
(8.50)
(8.51)

(8.52)

(8.53)

(8.54)

9:0 $53

In the mixed finite element model, it is necessary to specify the pressure at least at one
node. In the present case, the node at (y, z) = (L, O) is specified to have zero pressure. All
the fibers in the computational domain have been set initially random . At the initial time,
the velocity of the fluids and solids are set to be the same parabolic function of 2 as the
boundary condition at the inlet (BL).

In this problem the average Reynolds number is specified as 36 at the inlet and defor-
mation coefficient a = 0.2. The computational domain is meshed by 12 x 6 = 96 nine-node
quadratic elements for the velocity variable, and 12 x 6 = 96 four-node bilinear elements
for the pressure. Fiber orientation states at different time are shown in Figure 8.3. Dot
points at the inlet and near the center line indicate that the fibers are randomly oriented at
these regions. Away from the center line fibers become oriented faster and rotate inside
the fluids (i.e. xlmax is higher near the wall than the center region). The particle stresses
due to the presence of fibers depend on the orientation tensor a of the particles and the
strain rate of the fluids. The contour of main eigenvalue of the particle stress and its cor-
responding eigenvector are shown in Figure 8.4. It can been seen that the main eigenvalue
of the particle stress Tp is higher near the wall than the center region. Uniform suspension
at the initial time is assumed for this problem. At different times, the concentration ratio
a of the particles still keeps uniform shown in Figure 8.5. No migration is found for this
case. The pressure fields are shown in the Figure 8.6 at different times. The pressure field

is not same at any downstream cross section. The pressure is higher near the inlet than that

126

 

 

(a) t = 5.0
04 v .. Mx . '.'
N 0.2 ‘ ~ .
O r r r r L . L r i r . 1 l r
o 0.5 1 Y 1.5 2
(b) t = 10.0

 

 

04

 

 

.ZrLOBS

0.5
0.475
0.45
0.425
0.4
0.375
0.35

0.5
0.475
0.45
0.425
0 4

i 0.375

0.35

0.5
0.475
0.45
0.425

0.375

Figure 8.3. Contour plots of the principal eigenvalues Amax of the orientation te'nsor super-

posed with corresponding eigenvecotors for the problem of spheroids suspended in a plane

Poiseuille flow. The results are shown for three different times.

127

(a) t = 5.0 Tpmax

 

 

 

 

I
y--

I / I I I I /
/

/
‘ .r
’ ‘1
I / / I f / / I 0
4'1 Ll I l

 

 

 

 

 

0 g .r I. 2. 2 ,. .. ., , - -
o 0.5 1 v 1.5 2
(b)t=10.0 Tpmax

 

 

 

 

 

 

Tp max

0.0006
0.0005
0.0004
0.0003
0.0002
0.0001

 

 

 

 

 

 

Figure 8.4. Contour plots of the principal eigenvalues Tpmax of the particle stress super-
posed with corresponding main eigenvecotors for the problem of spheroids suspended in a

plane Poiseuille flow.

128

(a) t = 5.0 , a
0.01
0.01
0.01

: . 0.01

. 0.01

0.01

0.01

0.4

N 0.2

 

 

 

 

 

 

O

0.5 1 y 1.5 2

=1 .
(DH 00 0.01

0.01
0.01
0.01
0.01
0.01
0.01

0.4

N 0.2

 

 

 

 

 

 

O

0.5 r y 1.5 _2

(c) t = 15.0 a

0.01
0.01
0.01
0.01
0.01
0.01
0.01

0.4

 

I
N 0.2

 

 

 

 

 

 

Figure 8.5. Contour plots of concentration of the particles a for the problem of spheroids

suspended in a plane Poiseuille flow

129

0.4

No.2

 

 

0.4 _-

N 0.2

I I V I

 

 

0.4 -

 

N 0.2

 

 

00‘

 

Figure 8.6. Contour plots of the fluid pressure p f for the problem of spheroids suspended

in a plane Poiseuille flow

130

near the outlet. The specified parabolic velocity profile at the inlet is changed a little in the
downstream due to the particle stress, the fluid has some behaviors of non-Newtonian flow

to some extents.

131

CHAPTER 9

SUMMARY AND CONCLUSIONS

In the first part of this dissertation, various factors affecting the hydrodynamics of a single
particle suspended in a viscous fluid are studied. These factors include the presence of slip
on the particle surface, the influence of flow fields, non—Newtonian viscosity, and the pres-
ence of inertia forces. The drag force and rotary motion of a single particle are analytically
computed to study the effects of all these phenomena. In the second part, a closure model
for the orientation tensor of nearly arbitrary shape particles is developed and a framework
is proposed to estimate the alignment and spatial migration of spheroidal particles. The
findings associated with these studies are presented below.

The dynamics of a rigid particle shaped as a slightly deformed sphere surface in creep-
ing flows is studied with consideration of slip on the particle surface. Analytical expres-
sions are obtained for the hydrodynamic force and torque exerted by the fluid on a deformed
sphere using an asymptotic method wherein the normalized amplitude of the deviation from
sphericity is assumed to be a small parameter. The Stokes’ resistance calculated by this

method is validated by comparing with existing solutions of the limiting cases of no slip

132

and perfect slip. The analytical results for the axial and equatorial drag and torque on a
slightly deformed spheroid reproduce previously reported results for three limiting cases:
the perfect slip case, the no-slip case, and the case with an aspect ratio of unity (sphere).
This new theory has thus the potential to account for the presence of slip in multiphase
flows. In addition, the equations describing the motion of a deformed sphere with a slip
surface induced by a simple shear flow are also derived and solved. The motion of the
deformed sphere is shown to differ significantly from the no-slip case for low values of a
dimensionless parameter that incorporates the coefficient of sliding friction. The period of
the motion of a deformed sphere is longer, and for cases where the coefficient of sliding
friction is low, the spheroid rotates to a fixed angle and reaches a steady orientation state.

Analytical expressions for the drag force on a slightly deformed sphere suspended in
quadratic and cubic flows are derived by assuming that there is no slip on the interface
between the particle and the surrounding fluid. For a slightly deformed sphere, no rotary
motion is induced by the quadratic flow while periodic motion is induced by a cubic flow
with no slip. Comparison with the motion of a deformed sphere in a linear shear flow
reveals that the period of the motion of particle in a cubic flow is much longer if the same
coefficients as for the cubic flow are used.

Consideration of inertial forces on the uniform motion of a particle can be achieved for a
viscous flow by using Burgess’ general solution for Oseen flows. No-slip and slip boundary
conditions are considered on the interface between the particle and the fluid respectively.
Two kinds of geometry of the particle,. a sphere and a slightly deformed sphere, are studied.
Four cases are calculated respectively according to different boundary conditions on the

interface and the shape of the particle. They are: (l) the motion of a sphere with nonslip,

133

(2) the motion of a sphere with slip, (3) the motion of a deformed sphere with no slip,
and (4) the motion of a deformed sphere with slip. Both of the solutions for case (2), case
(3) and case (4) can be reduced to the solution of case ( 1) by letting the slip coeflicient
,B and the deformation coeflicient 8 equal to zero. The error between the solution of case
(4) when the slip coeflicient goes to co and the solution of case (3) is negligible when the
velocity and the diameter of the particle are small. The boundary condition at infinity is
well satisfied and the boundary condition at the interface is approximately satisfied if the
length dimension of the particle or the velocity of the particle is small.

It is well recognized that numerous fluids cannot be described by a Newtonian con-
stitutive model and the influence of a non-Newtonian viscosity is studied by allowing the
viscosity to vary with the shear rate. A Power-Law model is used to predict the viscosity
of the fluid. The no-slip boundary condition is also applied on the interface. It is found that
the non-Newtonian flow has much influence to the motion of a deformed sphere.

The consideration of particles of arbitrary shape as led to the development of a new
closure model to complete the description of the motion of ensembles of rigid particles of
complex shapes. Each particle is non-axisymmetric and its orientation is described with
a second order tensor < R >. An evolution equation for the second moment of the dis-
tribution function, which forms a fourth order tensor < R >, is used in order to obtain
the average orientation of the particles in homogeneous flows. As suggested by Rallison
(1978), the rotation operator is used to predict the microstructure of a suspension of rigid
, non-axisymmetric particles as a state variable rather than the Euler angles. This disserta-
tion has identified a closure for the 4th moment of the orientation distribution function in

terms of the 2nd moment that satisfies all six-fold symmetry and projection properties of

134

the exact the 4th moment.

In the last part of this work, models describing solid-liquid two-phase flows are de-
veloped using a continuum approach. A so-called Eulerian-Eulerian technique is adopted
to deal with the motion and migration of the particles and the fluid. Based on the mo-
ments of the distribution function, the evolution of the second moment of the orientation
tensor is used to govern the orientation of particles statistically. The concept of control vol-
ume/control surface method is used to develop closure models for the stresses and interfa-
cial force. The Fully Symmetric Quadratic model, developed at Michigan State University,
is applied to close the problem associated with computing the orientation tensor. A finite
element code is deve10ped to simulate the alignment and migration of dilute suspensions of
spheroids in a flowing liquid. Simulations results for flow between two parallel plates show
that at the inlet and near the center line the fibers are randomly oriented while away from
the center line fibers become oriented faster and rotate inside the fluids and no migration is

found for the plane Poiseuille flow.

135

APPENDICES

A. Tensor notation used in this dissertation

Dot product

Dot product of two vectors (a, b) is as the follows

3
a - b = Zaibi = a1b1+ azbz + a3b3 (A-l)

i=1

Dot product of two second-order tensor (A, B) is as follows

3 3 3
A . B = Z Z eiel [Z A,- 1231-1] (A-2)

A111311+ A1213214" A131331 A111912 + A12322 + A131332 A11313 + A121923 + A131933
= A211911 + A22321 + A231331 A211912 + A221922 + A231332 A211313 + A221323 + A23B33
A31311 + A32321 + A331331 A311912 + A321322 +4331932 A311313 + A321323 +A33B33

Double dot product

Double dot product of two second-order tensor (A, B) is as follows

3 3
A:B= 22.40.81,- (A-3)

i=1 j=1
= A111911+A12321 +A131331

+421 321 + A221322 + A231932

+4313 13 + A321323 + A33333

Double dot product of two fourth-order tensor (A, B) is as follows

3 3 3 3 3 3
A B = Z Z Z Z ......,.,[Z Z Arm...) («x-4)
=1

i=la=116= j=lfi=l

136

Dyadic product

Dyadic product of two vectors (x, y) is as follows

 

3 3
a®b = Zinyje,eJ (A-5)
i=1 j=l
albr 01192 6111?3
'-‘ azbr a2b2 a2b3
a3b1 (13192 (13193
Dyadic product of two second-order tensors (A, B) is as follows
3 3 3 3
A®B=ZZZZAiaBjfieieaejefl= (A-6)
i a j ,3
A11311 A111312 A111313 A121311 A121912 A121313 A131911 A131312 A13313
A111921 A111322 A111923 412321 A121322 A121323 A131321 A131922 A131923
A111331 A111332 A111333 A121331 A121332 A12333 A131331 A131332 413333
A213“ A21312 A211313 A22311 A22312 A22313 A23311 A231912 A23313
421321 421322 A21323 A221321 A221922 A221323 A231921 A231322 A231923
A211931 A211332 421333 A221931 A221932 A221333 A231931 A23332 A231933
A311311 A31312 A311913 A321311 A321312 A32313 A331311 A331312 A331313
A31321 A31322 A311323 A321921 A32322 A32323 A331321 A331922 433323
A311931 A31332 A311933 A321331 A321932 A321933 A331331 A331932 A331333

Dyadic product of two fourth-order tensors (A, B) is as follows

3 3 3 3 3 3 3 3
A®B=ZZZZZZZZAI'QIBBWMC,’eaej‘efiekeyeled (A-7)
r a j B k y l 6

137

 

B. Matrix form of the weak form of the governing equa-

tions
MassmatrixM
M11 o o 0 0 0 o 0 o
0 M22 0 o 0 o o 0 o
o 0 M33 o 0 0 0 0 0
0 o o M4'4OOO 0 0
o 0 0 0 0 0 o o 0
M: o 0 o 0 000 0 0
o o 0 0 o o o o o
0 o o o 000M833 0
o 0 o o o o o 0 M9"9
0 0 0 0 000 0 o M
o 0 0 o 0 o o o 0
(o 0 0 o 000 0 0

 

Entries of the mass matrix
M” = jg; stp,(w2,‘a)W,Tsdn
M22 = f9 stps(W;a)W;,l;dQ
M4.4 = jg; wvfpfu —w§a)wffdn
M83 = f WaWEdQ
9
M9’9 = f WaWEdQ
Q
M10,10 =f wawg‘dn
Q

Mll’” = jg; Wan/Eda

138

OOOOOOOOO

fl
9

CO

OOOOOOOOO

 

OOOOOOOOOOOO

(B-l)

(B-Z)

03-3)

03-4)

03-5)

03-6)

(B-7)

(B-8)

(B-9)

 

 

 

 

 

 

M12” zLWawg‘dQ (B-lO)
Stiffness matrix K
1 K“ K13 [(13 KL4 0 0 q 0 0 0 0 0
K11 K12 K23 K14 o 0 0 o 0 o 0
K3,] K3’2 [(3.3 K3,4 K3,5 O K3,7 0 O 0 K3,] 1
K4,1 K42 K433 K4,4 O K4,6 K47 0 O 0 [(4.11
O 0 K5,3 K5,4 K5’5 0 0 K5,8 K53 [(5.10 0
K_ o 0 K63 K94 0 K86 0 K68 K68 Km 0
— O 0 K7 ,3 K7,4 O 0 O K7,8 K7,9 K7’ 10 0
0 0 K8’3 K8’4 0 O 0 K8’8 K83 K8’10 0
o 0 K93 K” o 0 0 K93 K99 K930 0
O O K10,3 K10,4 O 0 0 K108 K103 K10,10 0
K1” K112 o 0 o o 0 0 0 0 Kll'”
\ O O K12,3 K12,4 O 0 O 0 0 0 K12’11
(B-ll)
Entries of the stiffness matrix
awT .awT
Km : vasp,a[(w3;vsy)—ay‘£+(Wfiv5z)—5z"—S)dn (B-12)
9(WTa)p
+ wT—"— (—5+28)—8(WTa ) (—wT)dQ
fa ”mafia—ed “ W] ”s
T
9(Wa’ a)p
K"2 =f WT—— .9 Wu.) (—wT dQ B-13
Q vs10r(2)(l—s)[( a )2] vs) ( )
K"3 I W Mk-smm-qw‘h )](—WT)dn (13-14)
52 V510r3(1-8) a .V) vf
'6st T (6W3 W] 6W3;
+ (W a)+ W a 2 d!)
Is 6y a "‘ 6y ) " 6y
' T ‘ T
6st T aWVf aWVf
+ — W + W d9
jg) aZ ( 00) vs[ dz p 62:

 

 

139

K1,12
[(2.12
K3,]2

K4,12

OOOOOOO

 

6st

 

Kl’4 = +faz
Q

 

(WT a>+ st[

10rg(1 -

T
K2,] : -Lstm[3 (W3 aw)“-

10rg(1 - )

ans
=+f stps( (“W 0) (stvsy)— (9y

VS lorgu - 8)

K2:3 = +f
Q

 

 

T 6Wv f
rL

K“ = + f vamp—s+2.)-.(wgau)](w3f)dn
Q

’ aWT

Vf

W1 62]

dz ]

1- w.»

rorgu — a)

aw... aw;
—6y- —(WT a) + st (3)—0]]

 

f.
f.

K212 = f
n

T

"67 —(Wa a)+w,,s(——6Z a

 

 

Lg” (WT )+ Wag.

aWT

7:11——

+ f stm [3(W3ayz)
Q 8).

+(WTsvsz)—

+15; W mk—S+28)—3(Wgazz)](—W3;)d§2

63:5 (W30) + st {-a—av-ng-an [+:1[

9<W§ )
+LWVSW£L18WTW2H1

 

 

+211 [—

 

T
6Wf

 

Wff] dQ

awT
az

 

- W..-(— W g-fwr f3 st<Wpf)dr

I‘U

K“ = + I W [-M] [(—5 + 28) — £(Wgayy)] (—w3;)dn
Q

10%(1 — e)

140

 

T
awvf

5y

W”) d9 vsy

 

 

d0

d9

i)dfl vsz

(B—15)

(B-16)

(B-17)

(B-18)

(B-19)

(B-20)

(B-21)

 

 

 

 

 

9(WTcr)p
K3’2=+ f W -——"-— -e WTa . -WT dQ 13-22
a ”f[ 10%(1-s)][ ( a ”W ”3) ( )
3 T T aWT T 6W3}
K' = Lvapf-(l-Waa) (vafy)-E‘i-+ (W fvfz)—— az d!) (B 23)
'aw ( WT ‘ BWT
f -—&y"—f(1—2W;,Ta)—2va{aayaaT pr gym
a_ J.
_ . T
anf ' T (WE ‘ 6va
jg;_ az (1-2W00)—2vaT az aj‘pf aZ d9

 

 

 

 

 

 

 

9(WTa)p T
jg; va [— mhe -5+—28) 5(W (1y),)](—WT f)dfl (B-24)
an T (9W3
K“ = fa 6Zf(1—2W§a)— 2va[;v" 0]],uf ayfdn

 

fW _ 9(Wga)fl
52 ”f largo—s)

f aWTp
K35 =+fQva -WTa

 

 

 

(
K3’7 = +Lvaa

aWT
6y

K3’11 : vaf|:-(W (WTPpry)——

6WT

paWTf
+1: vaflf —az

aWT
[(3.12 2 f
a

ay"—f-(1-2wuz—§) 2va
1L va(1—2W§a)

WW]

 

T
Wpde

141

dQ
K y

_(W;p Tm?)—

6W;

 

][ S(WTayZ)]( vaym

(B-25)
(B-26)

6W3

6y
aWT

V___f_
Ty vfy

d9

d0 (B-27)

6W;
6y

    

62:

[LWT (1)] (—Wpr)dQ (B-28)

9(Wga)p
K4,]z+1;W”f[—T0_rg(_1-_5][_8(WTG’Z)](—W3;)dfl (B-29)

Ta
1(4’2=+fQ va[--9—(XV——)—][(—5+2e)— 8(WTazz)](-Ws)dfl (B-30)

 

 

 

 

 

 

 

 

 

 

 

103(1
BWT T 8WT
,3 _ Vf T 6We Vf
K4 _ jg; ay (l—2Waa—) 2va[—6£a]pf[ az d!) (3.31)
9(Wga)p
+ W, —— - WT
jg; f[10r3(1—s)[8(“ayz)](w"f)
T _6WTf
- W 1—2W dr
.[I‘i Vf( aa)l1f—az
4 T T aWIf T 6W3
K4’ = vafpfo—Waa))vv——(vafy)w+(wfvfz) az d9 (B 32)
T
6va T [(9—wa ] anf
1—2Waa 2—W dQ
fm 6y ( "f ”f 6y
T
6W va T 6W; , anf
fm 6z (1 —2Wa :—az2wvf[ a)2pf az dQ

 

 

9(WT )
m mag-«Mn e< w

 

 

 

 

     

     

K46 ' T awTPT
. yr]; vaL(—Wa,a a) az dn (13-33)
K47 P T aWTP1
=+ In W (—Wa a) ay d0 (13-34)
aWT aWT
K4»ll = f va[— (WTPprZ)— ay —(W}przz) 701m (13-35)
6W5 anf aWT
-vaf(W (Tfpf) ——dQ+f vapr[a—z sz] 620
awff 6WTf aw;
+ vafpf 6y vfz+ az __va a),

 

142

   

if
I‘U

 

 

 

 

 

 

 

 

 

aWTf(

 

1—2W3a) — 2va[—

(WT
az

)]<- wzfm

awvf(1-2W;Ta)Wder— f3 W(1vf —2W TWa)T pde‘

 

 

 

 

 

 

 

 

 

 

. awT
W 2
3=-kf 5” —1+C2+2(1—W;Tayy-W§ayyazz) C2] avfdsz
Q L
W ' 6W3
Ki4 = —k a 7“ —1 +C2+2(1— WT aayy— WT ayyazz) )ZCZ] 62 d0
= f stwgpdsz
Q
T 6W3}
st (30—65C2 +94C2Wa aW)Ty-vfy T
K58 - ——k W d9
9 35 6wa “
+(5— 40C2+3SC2W Ta”) 62 __va
V awT W
5- 5C+35CWHa+4CWa ——v
K53 _ —k WW ( 2 2 W 2 a) fy >WaTdQ
Q 35 BWTf
1 +10(1—C2+7C2WTayy—C2WaTazz) az _va J
awff anTf
_— + _—
K5,10:_kf3st{ az vfy 6y vfz who
9 35 T awff aWTf
L+C2VVaayz -14—6y—ny+20'—éZ—'sz J
. . aWT
W 2 vf
K6'3——k ”3 —1 C +21—WT —WT C d9
(2 7 1 + 2 ( aayy aayyazz) 2‘ 6y
_ awT
W 2 1
K64 = —k 5” -1+ C2 + 2(1— W312” — WgTayyau) C2 avfdg
Q 1 . Z
K65: f WVSWT
Q
(9wa 1
T v
K68 st 10(1—C2—C2Waayy)—6y vfy T
’ = —k 35 6WT >WadQ
Q
L +-(5 5C2+4C2W draw) (92 —sz 1

143

 

(B-36)

(B-37)

(B-38)

(B-39)

(B-40)

(B-41)

(13-42)

(B-43)

03-44)

03-45)

(B-46)

K4,12 = L

+£Uwvf(1—2W§a)wadT- f3 W(vf1-2WT)WT pde‘

 

 

 

 

 

 

 

 

 

 

6WT aWT
ava(1 MW; )— 2va[——az— a]](- wwgf)

 

 

 

 

 

 

 

 

 

 

 

 

5 3 Qst_ T T 2 6W3}
W 6%,:
= I); WVSWTde
T 6W3}
5 8 ms (30 — 65C2 + 94C2Wa a”) __ay vfy T
K ’ = _k T Wa d9
Q 35 6W f
+—(5 4OC2 + 35C2W Ta”) 62 _va
r (5 5C 35C WT 4C WT )aW‘T ‘
— 2 + 2 0);); + 2 Ga ny
K5’9 — —k st< a a 6y]. >ngfl
Q 35 an
10 1— 7C WT C Wa
\ + ( C2 + 2(1-a”; 2 Tazz)— az —sz J
T T
W awva + anfv
K5 1°=-kf—"i 62 fy 6y fz >WaTdQ
Q 35 T T
T awvf awvf
+C2Wa ayz -14 6y ny + 206—2sz J
, 1 (WT
W 2
63 VS T T Vf
K - "'k Q 7 L—l + C2 + 2 (1 - W(1 ayy — Wa ayyazz) C2‘ 6)) ‘19
K6’4 st ' T T 2 ' 6W3}
2 —k 9 5 L—1+C2 +2(1- Waayy — Waayyazz) C2‘ 62 d9
= f WWW} d9
Q P
aWTf 1
10 1— C2 —C2W;Tayy #va
K63 = —k W” ( ) W T >ngQ
Q 35 6W f
+-(5 5C2 + 4C2W Ta”) 62 _va J

 

143

(B-36)

(B-37)

(B-38)

(B-39)

(B-40)

(B-41)

(B-42)

(B-43)

<B-44)

(B-45)

(B-46)

aWT

 

 

 

1 ”f
_ 1- 8C2 +14C2WTayy + 7C2WTazz _ny
[{6’9 = —k st 7( a a ) 6y 6WT Wg‘dQ
Q l —
+55- (30 — 65C2 + 35C2Wgayy + 94C2W3au) agfvfz
' (B-47)
' (5 — 5C2 - 10C2Wgayy + ZSCZWgaZZ) ‘
aWT BWT
3W vaf . + _ivf
K6,10 = _k vs 4 62 ) 6y 2 *WTdQ (BT48)
Q 35 T
awvaf
+C2WT ayz 20— _ 146W
K 6y

 

 

T

73 st T Czawzvf

K' =—k 35 —1+C2+2(1—Wa ayy— WT “My dQ (B-49)
Q

 

 

W 6153f
K” = —k f ”[ 1+C2 +2(1— WT a)”, WT (wawazzf Czay] d9 (1350)
Q

 

 

 

35

= jg; WVSWTPdQ (13-51)

7 8 st T awff 6W3} T
K ’ = -k 9 35 (5 - 5C2 — 8C2Wa a”) —a-Z—ny + —ay—sz Wa d9 (B-SZ)

7 9 st T T 6W3, aWTf T
K » = —k a 35 (5 -- 5C2 — 35C2Wa ayy — 8C2Wa au) a—zvfy + By _va Wa d!)
(353)
f 1 aWTf ‘
T V
7 (2 — 9C2 + 24C2Wa ayy + 3C2W3azz) '57va
awT

1
K7,10 = —k Wm +7 (2 — 9C2 + 24Czwgazz + 3C2W3ayy) 7%va >W3‘dn
T T
54 T awvf anf
+§§C2 (Wa ayz) {-E—ny + W172

 

 

k

 

 

(B-54)

2,1 2 3” f

8,3 _ T T V
K _fnwa{? [—1+C2+2C2(1—Waayy—Waazz)1 }dQ (3-55)

83 2’1 T T 6|le
K. = 9W0 —,-7—[—1+C2+2C2(1—Waayy—Waazz)2 62 an (B-56)

144

 

 

 

 

 

 

 

 

 

 

 

aWT aWT
8,8 _ T a T a
K _ L; W, (vavsy) __By + (vavsz) az ]dfl
( T )awff
+ f Wall a Q Wldn (13-57)
9 35 T anf
+ (5 — 8C2 + 7C2Wa ayy) '—a—Z—sz
T T 6W3"
K89 = + I W 2 (—5 — 5C2 + 35C2Wa a” + 4C2Wa an) 6y vfy WTdQ
Q 035 6W3} a
+ (10 — 10C2 + 70C2W3ayy — 10C2W3azz) 79?va
(B-58)
' T 6W3} 6W3} ‘
1211C2 (Wa ayz) 7—6y—vfy — 10—6z—vfz
1 . aWT
8,10 _ _ _ T
K ‘ L W“ 35 ‘ —5 [7 + a (—1 — 6C2 + 30C2W3ayy - 1.2C2W;Tazz)] ({va ’W0 m
an
+5 [7 + a (1 + 6C2 — 3002W3ayy + 12C2WaT (13)] vfy
k
(8-59)
9 3 24 T T 2 6W3].
K ’ .—. fawn 7[-1+C2 +2C2(1— Waayy — Waau)] 0y dQ (B-60)
aWT
2/1 2 f
9,4 __ T T V
K .. 1) Wu {—7— [—1 + C2 + 2C2(1- Waayy — Waazz)] az }dQ (B-61)
( T )6W}
-1() +10C2 +10C2W ayy —ny
K9,8 = + f Wag/l 0 ‘2y WIdQ (B-62)
Q 35 T (9va
(WT aWT
9,9 _ T a T a _
K _ L W, (vavsy) —6y + (vavsz) _az d9 (B 63)
T T 6W3
+ I W 2 ,1 (5 — 8C2 + 14C2Wa ayy + 7C2Wa au) vfy WTdQ
Q “35 aw} ‘1

 

+ (—5 - 65C2 + 35C2w3‘ayy + 94C2WaTazz) 62 sz

145

f

___v_fvyf+7___ vf

aWT aWT 1
6y az fz

121C2 (Wgayz) {—10—

9,10 _ _ _ T
K ‘ I, W“ 35 ‘ +5 [—7 + 1(1 + 6C2 + 12C2W3 ayy — 30C2WaT (13)] ’ Wa ‘19

 

foz
(B 64)

 

 

[ +5 [7 + x1 (1 + 6C2 + 12C2W3ayy — 30C2WaTazz)]a:;

10,3 24 T T 6W3}
K = 0W0 33[—1+(:2 +2C2(1— Waayy — Wa azz)2] 62 d9 (B-65)

 

10,4 24 T T 2aWTf
K = 9 W0 5 [—1 + C2 + 2C2(1— Wa ayy — Waazz)2 0y d9 (B-66)

 

 

awa
T
10,3 _ _ .1. {—35 + 1(15 + 20C2 + 32C2Wa aw)] az fvfy T
K _ W, 70 0W3} Wa d9
+ [35 + ,1(15 + 20C2 + 32C2Wa1-ayy)]— ayf"
(B-67)
1 1 T 116—‘15»
_35+,115+20C +32C W a
1 2 2 a yy
K109 = — f Wa7—O 3y? WIdQ (B-68)

Vf

+ [35 + 1(15 + 20C2 + 32C2WaTayy)] az vfy

 

 

 

 

aWT aWT
10,10 _ T
K _ jg; Wa (vavsy) 6y +(W vasz) 62 “[8152 (B-69)
aWT ‘
(—5 — 30C2 + 8OC2WTayy + 10C2W3‘azz) -6—yvf-vfy
T T
321 6W vf 8Wv f T
+f Wa— +36C2Wa T-—-ayz[ aZ -——"ny + —6y—sz >Wa d9
6W3}
+ (—5 — 3OC2 + 10(,*2W;“ayy + 80C2W} au) _371» fz
J

 

 

K11,1__ W WT (1’1ng _
— Q a( a a) 6y (B 70)
(WT

K11,2=_ fa wa(W}a)—azl’£d§2 (13-71)
aWT aWT
K1131: [9 W0, (st1)”)?y +(W,T_,vsz) a “[8192 (13-72)

 

146

T
K123 = —f W f(1— WTa)%dQ
Q p a (9Z

K12'4z—f W (1-WTa)%dQ
Q pf a 0y

33W aWT
[(12,112 fQWpf [(Wffoy) 3y (WT VfZ) $1“)

Force vector
F=(f1f2f3f400000000)T
‘ Entries of the force vector
f1: + f9 w,,p,(wga)g,m
f2 = + f9 wvsps(W§a)gzdn

f3 : +vafpf(1— WEa)gde

f1 = +IQWprf(1— Wga)gzd9

147

(B-73)

(B-74)

(B-75)

(B-76)

(B-77)

(B-78)

(B-79)

(B—80)

BIBLIOGRAPHY

[1] GB. Jeffery. The motion of ellipsoidal particles immersed in a viscous fluid. Proc.
Royal Soc. Math. London, 102:161—179, 1922.

[2] H. Brenner. The stokes resistance of an slightly deformed sphere. Chem. Engng.
Sci, 19:519—539, 1964.

[3] R. W. Burgess. The uniform motion of a sphere through a viscous liquid. American
Journal of Mathematics, 38:81—96, 1916.

[4] J. Happel and H. Brenner. Low Reynold number hydrodynamics with special appli-
cations to particulate media. Prentice-Hall,N. J ., 1965.

[5] S. Kim and S. J . Karrila. Microhydrodynamics: principles and selected applications.
Bulterworlt—Heineman, 1991.

[6] H. Lamb. Hydrodynamics. Dover, New York, 1945.

[7] B. J. Trevelyan and S. G. Mason. Force-free particles. Journal of Colloid Science,
6:354—370, 1951.

[8] F. P. Bretherton. The motion of rigid particles in a shear flow at low reynolds number.
Journal of Fluid Mechanics, 14:284—304, 1962.

[9] D. Broday, M.Fichman, M.Shapiro, and C.Gutfinger. Motion of spheroidal particles
in vertal shear flows. Phys. Fluids, 1:86—100, 1999.

[10] D. J. Jeffery. The calculations for teh low reynolds number resistance functions for
two unequal spheres. Phys. Fluids, 4:16—29, 1992.

[11] H. Keh and S. Chen. Low-reynolds-number hydrodynamic interactions in a sus-
pension of spherical particles with slip surfaces. Chemical Engineering Science,
52:1789—1805, 1997.

[12] E. D. Wetzel and C. L. TuckerHI. Droplet deformation in dispersions with unequal
viscosities and zero interfacial tension. Journal of Fluid Mechanics, 426:199—228,
2001.

[13] H. Niazmand and M. Renksizbulut. Surface effects on transient three-dimensional

flows around rotating spheres at moderte reynolds number. Computers and Fluids,
32:1405—1433, 2003.

148

[14] H. Brenner. The stokes resistance of an arbitrary particle. Chemical engineering
science, 27:1—25, 1963.

[15] H. Brenner and M. E. O’Neill. On the stokes resistance of multiparticle systems in
a linear shear field. Chemical engineering science, 27:1421—1439, 1972.

[16] G. K. Batchelor. Brownian diffusion of particles with hydrodynamic interaction.
Journal of fluid mechanics, 74:1—29, 1976.

[17] G. K. Batchelor. Sedimentation in a dilute polydisperse system of interacting
spheres. part 1. general theory. Journal of fluid mechanics, 119:379—408, 1982.

[18] D. J. Jeffery and Y. Onishi. Calculation of the resistance and mobility functions for
two unequal rigid spheres in low-reynolds-number flow. Journal of fluid mechanics,
139:261—290, 1984.

[19] E. J. Hinch and L. G. Leal. The effect of brownian motion on the rheological proper-

ties of a suspension of non-spherical particles. Journal of Fluid Mechanics, 52:683—
712, 1972.

[20] M. Doi and S. F. Edwards. The theory of polymer dynamics. Clarendon Press,
Oxford, 1987.

[21] G. K. Batchelor. The stress system in a suspension of force-free particles. Journal
of Fluid Mechanics, 41:545—570. 1970.

[22] E. J. Hinch and L. G. Leal. Constitutive equations in suspension mechanics, part 1.
general formulation. Journal of Fluid Mechanics, 71:481—495, 1975.

[23] D. Z. Zhang and A. Prosperetti. Momentum and energy equations for disperse two-
phase flows and their closure for dilute suspensions. International Journal of Multi-
phase Flow, 23:425—453, 1997.

[24] S. M. Dinh and R. C. Armstrong. A rheological equation of state for semiconcen-
trated fiber suspensions. Journal of Rheology, 28:207—227. 1984.

[25] W. R. Blakeney. The viscosity of suspensions of straight rigid rods. J. colloid
Interface Science, 22:324—330, 1966.

[26] R. O. Machmeyer and C. T. Hill. Trans. Soc. Rheol, 21:183—195, 1977.
[27] J. N. Miles, N. K. Murty, and G. F. Molden. Polym. Eng. Sci., 18:271—, 1981.

[28] T. Loimer, A. Nir, and R. Semiat. Shear-induced corrugation of free interfaces in
concentrated suspensions. Journal of Non-Newtonian Fluid Mechanics, 102:115—
134, 2002.

149

[29] G. K. Batchelor. The stress generated in a non-dilute suspension of elongated parti-
cles by pure straining motion. Journal of Fluid Mechanics, 46:813—829, 1971.

[30] J. D. Goddard. Tensile behavior of power-law fluids containing oriented slender
fibers. Journal of Rheology, 22:615—622, 1975.

[31] A. G. Gibson and S. Toll. Mechanics] of the squeeze flow of planar fiber suspensions.
Journal of Non-Newtonian Fluid Mechanics, 8211—24, 1999.

[32] H. Ramkissoon and St. Augustine. Slip flow past an approximate spheroid. Acta
Mechanica, 123:227—233, 1997.

[33] S. Deo and S. Datta. Slip flow past a prolate spheroid. Indian Journal of Pure
Applied Mathematics, 33:903—909, 2002.

[34] E. Bonaccurso, M. Kappl, and H.-J. Butt. Hydrodynamic force measurements:
Boundary slip of water on hydrophilic surfaces and electrokinetic effects. Phys.
Rev. Lett., 88:076103, 2002.

[35] V. S. J. Craig, C. Neto, and D. R. M. Williams. Shear-dependent boundary slip in an
aqueous newtonian liquid. Phys. Rev. Lett., 87:054504, 2001.

[36] R. Pit, H. Hervet, and L. Léger. Direct experimental evidence of slip at hexade-
cane/solid interfaces. Physical Review Letters, 85:980, 2000.

[37] PG. de Gennes. On fluid/wall slippage. Langmuir, 18:3413—3414, 2002.

[38] S. Lichter, A.Roxin, and S. Mandre. Mechanisms for liquid slip at solid surfaces.
Physical Review Letters, 93:086001—1—086001-4, 2004.

[39] M. Vincent, E. Devilers, and J. F. Agassant. Fiber orientation in injection moulding
of reinforced therrnoplastics. J. Non-Newtonian Fluid Mech., 73:317—326. 1997.

[40] P. A. Thompson and SM. Troian. A general boundary condition for liquid flowat
solid surfaces. Nature (London), 389:360—362, 1997.

[41] S. Granick, Y. Zhu, and H. Lee. Slippery questions about complex fluids flowing
past solids. Nature Mater, 2:221, 2003.

[42] J. R. A. Pearson and C. J. S. Petrie. Polymer Systems: Deformation and Flow.
Macmillian, London, 1968.

[43] S. Richardson. On the no-slip boundary condition. J. Fluid Mech., 59:707—719,
1973.

[44] L. M. Denn. Issues in viscoelastic fluid mechanics. Annu. Rev. Fluid Mech., 22:13—
34, 1990.

150

[45] N. V. Priezjev, A. A. Darhuber, and S. M. Troian. Slip behavior in liquid films on
surfaces of patterned wettability: Comparison between continuum and molecular
dynamics simulations. PHYSICAL REVIEW E, 71:0416081—041608-1 1, 2005.

[46] P. Ball. How to keep dry of water. Nature (London), 423:25, 2003.

[47] A. Karlsson. Molecular engineering: Networks of nanotubes and containers. Nature
(London), 409:150, 2001.

[48] E. Schnell. Slippage of water over nonwettable surfaces. J. Appl. Phys, 27:1149,
1956.

[49] O. I. Vinogradova. Slippage of water over hydrophobic surfaces. Int. J. Min. Process,
56:31, 1999.

[50] K. Watanabe, Y. Udagawa, and H. Udagawa. Drag reduction of newtonian fluid in a
circular pipe with a highly water-repellent wall. J. Fluid Mech., 381 :225—238, 1999.

[51] C. Cottin-Bizonne, J. L. Barrat, L. Bocquet, and E. Charlaix. Low-friction flows of
liquid at nanopattemed interfaces. Nat. Mater, 2:237—240, 2003.

[52] S. Richardson. On the no-slip boundary condition. J. Fluid Mech., 59:707—719,
1973.

[53] L. M. Hocking. A moving fluid interface on a rough surface. J. Fluid Mech., 76:801—
817, 1976.

[54] K. M. Jansons. Determination of the macroscopic (partial) slip boundary condition
for a viscous flow over a randomly rough surface with a perfect slip microscopic
boundary condition. Phys. Fluids, 31:15, 1988.

[55] D. Einzel, P. Panzer, and M. Liu. Boundary condition for fluid flow: Curved or rough
surfaces. Phys. Rev. Lett., 64:2269, 1990.

[56] M. J. Miksis and S. H. Davis. Slip over rough and coated surfaces. J. Fluid Mech,
273:125, 1994.

[57] L. Bocquet and J. L. Barrat. Hydrodynamic boundary conditions, correlation func-
tions and kubo. relations for confined fluids. Phys. Rev. E, 49:3079, 1994.

[58] H. A. Barnes. A review of the slip (wall depletion) of polymer solutions,emulsions
and particle suspensions in viscometers. J. Non-Newtonian Fluid Mech., 56:221—
251, 1995.

[59] D. Andrienko, B. aneg, and O. I. Vinogradova. slip as a result of a prewetting
transition. J. Chem. Phys, 119:13106—13112, 2003.

151

[60] M. M. Denn. Extrusion instabilities and. wall slip. Annu. Rev. Fluid Mech., 33:265-
287, 2001.

[61] K. B. Migler, H. Hervet, and L. Leger. Slip transition of a polymer melt under
shear-stress. Phys. Rev. Lett., 70:287—290, 1993.

[62] V. Mhetar and L. A. Archer. Slip in entangled polymer melts. 1. general. features.
Macromolecules, 31 :8607—8616, 1998.

[63] C. H. Choi, K. J. A. Westin, and K. S. Breuer. Apparent slip flow in hydrophilic and
hydrophobic. microchannels. Phys. Fluids, 1522897, 2003.

[64] P. A. Thompson and S. M. Troian. A general boundary condition for liquid flow at
solid surfaces. Nature (Londond), 389:360—362, 1997.

[65] N. V. Priezjev and S. M. Troian. Molecular origin and dynamic behaviour of slip in
sheared. polymer films. Phys. Rev. Lett., 92:018302, 2004.

[66] R. G. Horn, 0. I. Vinogradova, M. E. Mackay, and N. Phan-Thien. Hydrodynamic
slippage inferred from thin film drainage measurements in a solution of nonadsorb-
ing polymer. J. Chem. Phys, 112:6424—6433, 2000.

[67] Y. Zhu and S. Granick. Rate-dependent slip of newtonian fluids at smooth surfaces.
Phys. Rev. Lett., 87:096105, 2001.

[68] O. I. Vinogradova, N. F. Bunkin, N. V. Churaev, O. A. Kiseleva, A. V. Lobeyev, and
B. W. Ninham. Submicrocavity structure of water between hydrophobic and by-
drophilic walls as revealed by optical cavitation. J. Colloid Interface Sci., 1732443-
447, 1995.

[69] H. Brenner. The stokes resistance of an arbitrary particle ii an extension. Chem.
Engng. Sci., 19:599—629, 19643.

[70] H. Brenner. The stokes resistance of an arbitrary particle iii shear fields. Chem.
Engng. Sci., 19:631—651, 1964b.

[71] H. Brenner. Coupling between the translational and rotational brownian motions of
rigid particles of arbitrary shape: Ii general theory. J. Colloid Interface Sci., 28:407—
436, 1967.

[72] I. L. Claeys and J .F. Brady. Suspensions of prolage spheroids in stokes flow. part
1. dynamics of a finite number of particles in an unbounded fluid. J. Fluid Mech.,
251:411—442, 1993.

[73] S. G. Advani and C. L. Tucker H1. The use of tensors to describe and predict fiber
orientatioin in short fiber composites. Journal of Rheology, 31:751—784, 1987.

152

[74] TM. Macrobert. Spherical Harmonics. Dover, New York, 1947.
[75] AB Basset. A Treatise on Hydrodynamics. Dover, New York, 1961.

[76] G. G. 11 Lipscomb and M. M. Denn. The flow of the fiber suspensions in complex
geometries. Journal of Non-Newtonian Fluid Mechanics, 26:297—325, 1988.

[77] R. B. Bird, R. C. Armstrong, and O. Hassager. Dynamics of polymeric liquids,
Volume 2, Fluid Mechanics. John Wiley and Sons, 1987.

[78] P. K. Kundu. Fluid mechanics. Academic Press, Inc, 1990.

[79] T. E. Garstang. The flow of viscous liquid past spinning bodies. Proceeedings of
the Royal Society of London, Series A, Containing Papers of a Mathematical and
Physical Character, 1422491—508, 1933.

[80] S. Goldstein. On the two-dimensional steady flwo of a viscous fluid behind a solid
body ii. Proceeedings of the Royal Society of London, A, 142:563-537, 1933.

[81] L. Bairstow, B. M. Cave, and E. D. Lang. The resistance of a cylinder moving in a
viscous fluid. Philosophical Transactions of the Royal Society of London. Series A,
223:383—432, 1923.

[82] P. V. Southwell and H. B. Squire. A modification of oseen’s approximate equation
for the motion in two dimensions of a viscous incompressible fluid. Philosophical
Transactions of the Royal Society of London. Series A, 232:27—64, 1933.

[83] P. A. Davies, D. L. Boyer, H. J. S. Fernando, and X. Zhang. On the periodic motion
of a circular cylinder through a linearly stratified fluid. Philosophical Transactions:
Physical Sciences and Engineering, 346:353—386, 1994.

[84] A. Berry and L. M. Swain. On the steady motion of a cylinder through infinite
viscous fluid. Proceeedings of the Royal Society of London, A, 102:766—, 1923.

[85] D. Meksyn. Solution of oseen’s equations for an inclined elliptic cylinder in a vis-
cous fluid. Proceeedings of the Royal Society of London, A, 162:232—251, 1937.

[86] G. J. Richards. On the motion of elliptic cylinder through a viscous fluid. Phy-
ilosiphical Transaction of the Royal Society of London, A, 233:279—301, 1934.

[87] N. Piercy and H. F. Winny. The skin friction of flat plates to oseen’s approximation.
Proceeedings of the Royal Society of London, A, 140:543—561, 1933.

[88] K. K. Kabanemi, J. F. Hetu, and A. G. Rejon. A 3-d coupled solution for the flow
and fiber orientation in injection molding fo fiber-filled systems. Rheology and Fluid
Mechanics of Nonlinear Materials, 217:207—227, 1996.

153

[89] R. S. Bay and C. L. Tucker HI. Fiber orientation in simple injection moldings, part
1: Theory and numerical methods. Polymer Composites, 13:317—321, 1992.

[90] R. S. Bay and C. L. Tucker 111. Fiber orientation in simple injection moldings, part
2: Experimental results. Polymer Composites, 13:322—342, 1992.

[91] H. H. de Frahan, V. Verleye, F. Dupret, and M. J. Crochet. Numerical prediction of
fiber orientation in injection molding. Polym. Eng. Sci., 32:254—266, 1992.

[92] M. Gupta and K. K. Wang. Fiber orientation and mechanical properties of shor-
fiber-reinforced injection-molded composites: simulated and experimentatl results.
Polymer Composites, 147:367—382, 1993.

[93] E. J. Hinch and L. G. Leal. Time-dependent shear flows of a suspension of particles
with weak brownian rotations. Journal of Fluid Mechanics, 57:753—797, 1973.

[94] J. Ko and J. R. Youn. Prediction of fiber orientation in the thickness plane during
flow molding of short fiber composites. Ploym. Campos, 16:114—124, 1995.

[95] F. Folgar and C. L. Tucker III. Orientation behavior of fibers in concentrated sus-
pensions. J. Reinf Plast. Compos, 321095—1122, 1995.

[96] A. C. Papanastasiou and A. N. Alexandrou. Isothermal extrusion of non-dilute fiber
suspensions. J. Non-Newtonian Fluid Mech., 25:313—328, 1987.

[97] G. Ausias, J. F. Agassant, and M. Vincent. Flow and fiber orientation calculations in
reinforced therplastic extruded tubes. Int. Polymer Processing, 9:51—59, 1994.

[98] E. J. Hinch and L. G. Leal. Constitutive equations in suspension mechanics, part 2.
approximate forms for a suspension of rigid particles affected by borwnian rotations.
Journal of Fluid Mechanics, 76: 187—208, 1976.

[99] G. L. Hand. A theory of anisotropic fluids. Journal of Fluid Mechanics, 13:33-46,
1962.

[100] S. G. Advani and G. L. Tucker. Closure apporximations for three-dimensional struc-
ture tensors. Journal of Rheology, 34:367—386, 1990.

[101] J. S. Cintra and C. L. Tucker. Orthotropic closure approximation for flow-induced
fiber orientation. Journal of Rheology, 39:6—, 1995.

[102] S. M. Parks and C. A. Petty. Prediction of Low-Order Orientation Statistics for
F low-Induced Alignment of Fibers and Platelets. paper presentation, Fundamental
Research in Fluid Mechanics: Particulate and Multiphase Flow 1, AIChE Annual
Meeting, Dallas, TX, October 31- November 5, 1999a.

154

[103] S. M. Parks, C. A. Petty, and M. Shafer. F low—Induced Alignment of Fibers in the Ab—
sence of Fiber-Fiber Interactions. paper presentation, Symposium on Suspensions,
APS/DFD, New Orleans, LA, November 21-23, 1999.

[104] A. Acrivos and C. Taylor. The stokes flow past an arbitrary particle. Chem. Engng.
Sci., 19:445, 1964.

[105] H. J. Workman. Concerning the orientation distribution function of rigid particles
in a suspension which is undergoing simple shear flow. J. Colloid and Interface
Science, 29:664—669, 1969.

[106] J. M. Rallison. The effect of brownian rotations in a dilute suspension of rigid parti-
cles of arbitrary shape. J. Fluid Mech., 84:237—263, 1978.

[107] R. G. Larson. Consititutive equations for polymer melts and solutions. Butterworth,
1988.

[108] W. C. Jackson, A. D. Advani, and C. L. Tucker III. Predicting the orientation of short
fibers in thin compression moldings. Journal of Composite Materials, 20:539—557,
1986.

[109] M. C. Altan, S. G. Advani, S. I. Guceri, and R. B. Pipes. On the description of the
orientation states for fiber suspensions in homogenous flows. Journal of Rheology,
33:1129—1155, 1989.

[110] J. Betten. Integrity basis for a second-order and a fourth-order tensor. Intemat. J.
Math, 5:87—96, 1982. i

[111] Y. Kim, C. T. Nguyen, S. M. Parks, D. Mandal, A. Bnard, and C. A. Petty. Mi-
crostructure of Liquid Crystalline Polymers Induced by Simple Shear. Symposium
on Manipulation of Nanophases by External Fields,Annua1 AIChE Meeting, Indi-
anapolis Convention Center, Indianapolis, 2002, November 3 - 8.

[112] Y-C Kim, H. Kini, D. Mandal, A. Bnard, and C. Petty. Microstructure of Liq-
uid Crystalline Polymers Induced by Simple Shear. Symposium on Suspensions,

56th Annual Meeting, American Physical Society: Division of Fluid Dynamics,East
Rutherford, NJ, 2003, November 23-25.

[113] D. Leighton and A. Acrivos. The shear-induced migration of particles in concen-
trated suspensions. Journal of Fluid Mechanics, 181:415-439, 1987.

[114] J. F. Brady and G. Bossis. The rheology of concentrated suspensions of spheres in
simple shear flow by numerical simulation. Journal of Fluid Mechanics, 155: 105—
129, 1985.

155

[1 15] J. F. Brady and G. Bossis. Stokesian dynamics. Annual Review of Fluid Mechanics,
20:111—157, 1988.

[116] D. Chen and M. Doi. Simulation of aggregating colloids in shear flow. Journal of
Chemical Physics, 91 :2656—2663, 1995.

[117] N. Phan-Thien and S. Kim. Microstructures in Elastic Media: Principles and Com-
putational Methods. Oxford, New York, 1994.

[118] M. S. Ingber, A. A. Mammoli, and L. A. Mondy. A parallel obundary elitent formu-

lation for determining effective properties of heterogeneous media. J. Numer. Meth.
Eng, 37:3905—3919, 1994.

[119] G. Segre and A. Silberberg. Radial poiseuille flow of suspension. Nature, 189:209,
1961.

[120] G. Segre and A. Silberberg. Behavior of macroscopic rigid spheres in poiseuille
flow. part 1. Journal of Fluid Mechanics, 14:385—413, 1962.

[121] J. H. Aubert and M. Tirrell. Macromolecules in nonhomogeneous velocity gradient
fields. J. Chem. Phys, 72:2694—2701, 1980a.

[122] J. H. Aubert, S. Prager, and M.Tirre11. Macromolecules in non-homogenous. veloc-
ity gradient field. part ii. J. Chem. Phys, 73:4103—4112, 1980.

[123] R. C. Armstrong G. Sekhon and M. S. Jhon. The origin of polymer migration in a
nonhomogeneous flow field. Journal of Polymer Science: Polymer Physics Edition,
20:947—952, 1982.

[124] C. Crowe, M. Sommerfeld, and Y. Tsuji. Multiphase flows with droplets and parti-
cles. CRC Press, Boca Raton, FL, 1998.

[125] S. S. Sadhal, P. S. Ayyaswamy, and J. N. Chung. Transport phenomena with drops
and bubbles. Springer, New York, 1997.

[126] J. F. Davidson, R. Clift, and D. Harrison. F luidization. Academic Press, Orlando,
1985.

[127] L. S. Fan and C. Zhu. Principles of gas—solidfiows. Cambridge University Press,
Cambridge, 1998.

[128] R. D. Marcus, L. S. Leung, G. E. Klinzing, and F. Rizk. Pneumatic conveying of
solids. Chapman and Hall, London, 1990.

[129] M. Ungarish. Hydrodynamics of suspensions. Springer, New York, 1993.

156

[130] M. Massoudi. Constitutive relations for the interaction force in multicomponent par-
ticulate flows. International Journal of Non-Linear Mechanics, 38:313—336, 2003.

[131] T. B. Anderson and R. Jackson. Fluid mechanics] description of fluidized beds. Ind.
Engng. Chem. Fundam, 6:527, 1967.

[132] R. Jackson. Locally averaged equations of motion for a mixture of identical spheri-
cal particles and a newtonian fluid. Chemical Engineering Science, 52:2457—2469,
1997.

[133] D. Drew. Mathematical modeling fo two-phase flow. Annual Review of Fluid Me-
chanics, 15:261—291, 1983.

[134] D. D. Joseph and T. S. Lundgren. Ensemble averaged and mixture theory equations
for incompressible fluid-particle suspensions. International Journal of Multiphase
flow, 1:35—42, 1990.

[135] D. A. Drew and S. L. Passman. Theory of multicomponent fluids. Apringer-Verlag,
New York, 1998.

[136] M. Ishii. Thenno-fluid dynamic theroy fo two—phase flow. Eyrolles, Paris, France,
1975.

[137] D. Lakehal and C. Narayanan. Numerical analysis of the continum formulation
for the initial evoluiton of mixing layers with particles. International Journal of
Multiphase Flow, 29:927—941, 2003.

[138] D. Z. Zhang and A. Prosperetti. Average euqations for inviscid disperse two-phase
flow. Journal of Fluid Mechanics, 267:185-219, 1994.

[139] T. Lendgren. Slow flow through stationalry random beds and suspensions fo sheres.
Journal of Fluid Mechanics, 51:273—299, 1972.

[140] D. A. Drew. A turbulent dispersion model for particles or bubbles. Journal of
Engineering Mathematics, 41:259—274, 2001.

[141] G. J. Hwang and H. H. Shen. Modeling the solid phase stress in a fluid-solid mixture.
Int. J. Multiphase Flow, 15:257—268, 1989.

[142] G. J. Hwang and H. H. Shen. Modeling the phase interaction in the momentum
equations fo a fluid-solid mixture. Int. J. Multiphase Flow, 1:45-57, 1991.

[143] G. K. Batchelor and J. T. Green. The determination of the bulk stress in a suspension
of spherical particle to order c2. Journal of Fluid Mechanics, 56:401-472, 1972.

[144] P. G. Saffman. The lift on a small sphere in a slow shear flow. Journal of Fluid
Mechanics, 22:385—400. 1965.

157

[145] D. K. Mandal. Similation of flow-induced fiber orientation with a new closure model
using the finite element method. PH.D dissertation, Michigan State University, 2004.

[146] J. N. Reddy and D. K. Gartling. The finite element method in heat transfer and fluid
dynamics. CRC press, 2001.

[147] F. Brezzi and K. J. Bathe. A discourse on the stability conditions for mixed finite

element formulations. Computer Methods in Applied mechanics and Engineering,
82, 1990.

158

   

n]‘u]‘]j]“n]][][I]