Dynamics of Colloidal and Polymeric Suspensions via a Fluctuating Lattice-Boltzmann Model

Material Information

Dynamics of Colloidal and Polymeric Suspensions via a Fluctuating Lattice-Boltzmann Model
USTA, O. BERK ( Author, Primary )
Copyright Date:


Subjects / Keywords:
Center of mass ( jstor )
Diffusion coefficient ( jstor )
Fluids ( jstor )
Hydrodynamics ( jstor )
Momentum ( jstor )
Monomers ( jstor )
Peclet number ( jstor )
Polymers ( jstor )
Simulations ( jstor )
Velocity ( jstor )

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright O. Berk Usta. Permission granted to University of Florida to digitize and display this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Embargo Date:
Resource Identifier:
660034147 ( OCLC )


This item has the following downloads:

usta_o ( .pdf )








































































































Full Text







Copyright 2007


0. Berk USTA

To my family, friends and that somebody I do not know or recognize yet, and .,_1:one else

who will appreciate this dissertation.

We're all lonely for something we don't know we're lonely for. How else to explain the

curious feeling that goes around feeling like missing somebody we've never even met?

-David Foster Wallace, Loneliness


I would like to thank my advisors Jason E. Butler and Anthony J.C. Ladd for their

professionalism, tremendous knowledge, patience, support and help during the last five

years. I have disagreed with them on many occasions and still do. However we all agree

that we can disagree which I found one of the most valuable things during my doctoral

degree. However, I do wish my disagreements to have been on more solid grounds. I am

yet to completely disprove anything they have told me except for small issues. It was a

great pleasure and also a distinct privilege to have worked with two such quality people

so different yet so similar. I am grateful to Dr. James Dufty not only for his presence on

my committee but also for his lectures, kindness, sincerity, wisdom and great example of

how a scientist should be. I would also like to thank my other committee members Dr.

Richard Dickinson, and Dr. Sergei Obukhov and Dr. C'!5 ,i:-Wong Park for their time and

advice. I am also grateful to my undergraduate advisor Dr. Tiirkan Haliloglu who still

provides a her valuable support and friendship. I would now like to turn to my colleagues

Dr. Piotr Szymczak, Dr. Jonghoon Lee and Dr. Byoungjin Chun for being members of

the original lunch gang. As Jonghoon put it in his thesis \\ discussed subjects without

boui. 1 ,i Piotr has probably been one of the most influential person in my life besides

my family, my friends and my advisors. I wish one d4v I can be half as good as he is in

being able to reduce complex problems into simple terms, among all the other things he

does so well. Jonghoon and Byoungjin have both been great inspiration to me by their

calmness, dedication and hardwork. I will never forget neither their personal nor their

professional help and support. I would also like to thank Jonathan M. Bricker for many

nights we have spent together in the lab sharing not only science but our lives and Quyen

Nguyen for helping me both personally and professionally when I first joined the group. I

would also like to acknowledge Gaurav Misra who, within the limited time we knew each

other, impressed me greatly with his focus and love for science. I would also like to thank

Joontaek Park, Hyun-Ok Park, Phillip D. Cobb, Murali Rangareai and Dazhi Yu.

I guess the most important people in shaping my life and me as a person would be

my family and my greatest gratitude goes to all of them, my mom, Olcay Usta, my dad.

M. Nafiz Usta, and my sister, Burcu Usta Uslu. Besides our frequent conversations and

their visits and continuing financial support they, too, have provided me what I value most

dearly; to disagree. I have not become what they probably have envisioned for me nor

what I have envisioned for my self many years ago. However, I do feel their support and

love with me regardless of this and I hope they will continue to do so even if I do fail every

once in a while.

The presence of my friends Kerem Uguz, Ozgiir Ozen, Derya Gtilsen, Monica James,

Milorad Dmjoiilii was invaluable to me both emotionally and intellectually. Dr. Anuj

C'! i ,,!1i i was another big support whom I can talk to about anything and has enriched

my life and visions. Although, we have gotten to know each other late I would like to

acknowledge here, Heng Zhu, Chi-Chuing Li, Vicky Huang, Darren McDuff, Colin Sturm,

Young Sun, Young-Seok Kim, Chia-Yi C'!, i, and Andrew Kuo. I would also like to thank

my friends here who have not been here with me physically during my degree but were

with me emotionally. Just knowing I could talk to them every once in a while through the

speaker or the keyboard was a big relief and really kept me going. To name a few, I should

start with my best friend Aykut Arg6ntil. He's been the one that I talk to most when

I have problems, the true definition of "' .-I when I needed to. Emek Seyrek is a close

second, whom I have known for a long long time and had so much fun with over those

many years. G6kge Yilancioglu, Delal Dink, Aysa Akad, Hande Tiirgak, Aysegiil Ayok,

Mehmet Mercang6z, Cem Salahifar, Cem Dinger, and Asim Okur are a few others to name

here. I apologize to those whom I left out but love so dearly.


ACKNOW LEDGMENTS ................................. 4

LIST OF FIGURES .. .. .. ... .. .. .. .. ... .. .. .. ... .. .. .. 8

A B ST R A C T . . . . . . . . . . 12



1.1 Introduction . . . . . . . . .14
1.2 Fluid Motion and Hydrodynamic Equations ................. 15
1.3 Hydrodynamic Interactions ........................... 17
1.4 Numerical Methods to Solve Hydrodynamic Equations .......... 19

SU SPEN SIO N S . . . . . . . . . 24

2.1 Boltzm ann Equation .............................. 24
2.2 Lattice Boltzmann Method ........................... 26
2.2.1 An Overview of Lattice-Boltzmann Method ............. 26
2.2.2 Lattice Boltzmann Equation ...................... 28
2.3 Fluctuations . . . . . . . . 30
2.4 Solid Particles . . . . . . . . 31
2.4.1 Finite-Size Spherical Particles . . . . . 31
2.4.2 Point Particles . . . . . . . 33
2.4.3 Particle M otion . . . . . . . 36
2.5 Verification of the Fluctuations . . . . . . 37

GEOMETRIES . . . . . . . . . 43

3.1 Introduction . . . . . . . . 43
3.2 R results . . . . . . . . 45
3.2.1 Fluctuation-Dissipation . . . . . . 45
3.2.2 Hydrodynamic Interactions . . . . . 50
3.2.3 Diffusion of an Isolated Polymer . . . . . 52
3.2.4 Polymer Diffusion in a C ii, . . . . . 57
3.3 Conclusions . . . . . . . . 59


4.1 Introduction . . . . . . . . 62
4.2 M ethod . . . . . . . . . 63
4.3 Results and Discussion . . . . . . . 65
4.4 Conclusion . . . . . . . . . 73


5.1 Introduction . . . . . . . . 75
5.2 M ethods . . . . . . . . . 76
5.3 R results . . . . . . . . . 76
5.4 C conclusion . . . . . . . . . 85

6 CONCLUSIONS .. .. .. ... .. .. .. ... .. .. .. .. ... .. .. 86


A DISCRETE LANGEVIN EQUATION ....................... 88


T H E O RY . . . . . . . . . 93

T H E O RY . . . . . . . . . 94

REFERENCES ....................................... 96

BIOGRAPHICAL SKETCH ................................ 102

Figure page

2-1 Location of the boundary nodes for a circular object of radius 2.5 lattice spacings.
The velocities along links cutting the boundary surface are indicated by arrows.
The locations of the boundary nodes are shown by solid squares and the lattice
nodes by solid circles . . . . . . . . 32

2-2 Momentum change correlations on a plane wall. AP(t) is the total momentum
exchange between the fluid and the particle at any instant t. Results from simulations
are compared with a theoretical expression derived from continuum theory. The
theoretical expression is higher than the simulation data at the first time-step,
error accumulates for a couple of time-steps and then the difference between
these two curves turn into a constant offset. The first point for this correlation
is calculated correctly (l1 ) from a discrete theory starting with the fluctuations
in population densities (section D) . . . . . . 39

2-3 Momentum exchange correlations on an infinitely massive sphere. AP(t) is the
total momentum exchange between the fluid and the particle at any instant t.
Results from simulations are compared with a theoretical expression derived
from continuum theory . . . . . . . . 40

2-4 Verification of Onsager's linear response theory. The decay of the velocity correlations
is compared with the velocity decay of a particle with an initial given velocity.
The two curves lead to the same diffusion coefficient upon integration. . 41

2-5 Decay of the velocity autocorrelation function for a single particle with mass
factor of 1 & 10 and a radius of 5.4, normalized by the equipartition value U2
kT/mB. The two different curves correspond to different integration schemes.
These two curves eventually lead to the same diffusion coefficient. . . 42

3-1 Relation between the effective friction coefficient and the input friction o (Eq. 2-21).
The effective friction coefficient has been obtained from the drag force on a single
monomer (circles) and the diffusion of the monomer in a thermally fluctuating
fluid (squares). The slope of the line is unity, indicating a constant offset between
-1 and o1-. The unit cell size L 20Ax. . . . . ... 46

3-2 Offset between -1 and io1 as a function of system volume. The variation in
the dimensionless mobility g (Eq. 3- 1) is shown for a cubic cell of length, L.
Results are shown for input hydrodynamic radii ao = 0.32 (circles) and ao =
0.15 (squares). The monomer friction 0o = ',.,/ . . . . ... 48

3-3 Mean kinetic energy of a single monomer as a function of the dimensionless frictional
coupling, CAt/m. The calculated kinetic energy for implicit (circles) and explicit
(squares) updates is compared with the model described in Appendix A (lines). 49

3-4 Velocity due to a point force. The velocity component parallel to the direction
of the force, vx(r, 0, 0), is plotted as a function of distance from the force. Results
are shown for the Oseen flow field (dashed), the Stokes flow field in a periodic
cell of length 100Ax (solid line) and lattice-Boltzmann simulations (circles) in
the sam e size cell . . . . . . . . 51

3-5 Comparison of the parallel component of the two particle mobility tensor, pl2.
Results obtained by dissipative (Eq. 3-10, circles) and fluctuating (Eq. 3-11,
squares) simulations are compared with an independent solution of the periodic
Stokes equations (Eq. 3-9, triangles). . . . . . . 52

3-6 Comparison of the perpendicular component of the two particle mobility tensor,
pf2. Results obtained by dissipative (Eq. 3-10, circles) and fluctuating (Eq. 3-11,
squares) simulations are compared with an independent solution of the periodic
Stokes equations (Eq. 3-9, triangles). . . . . . . 53

3-7 End-to-end distance and radius of gyration in a periodic unit cell. The scaling
of Re (upper points) and Rg (lower points) is shown for different monomer separations,
b, and temperatures, T. The dashed lines are the theoretical scalings for an excluded
volume chain, with exponent v = 0.59. . . . . . 54

3-8 Center of mass diffusion coefficient as a function of chain length in different size
cells. The diffusion coefficients are normalized by the monomer diffusion coefficient,
Dm = T/l, and are shown for periodic systems with lengths L ~ 3.5R, and L
~ 7R,. In all cases the mean separation between the monomers b Ax and the
temperature T = 0.1. Results are shown for the raw diffusion coefficient (open
symbols), as well as the corrected values (filled symbols). . . ..... 55

3-9 Center of mass diffusion coefficient, with finite-size corrections, for different values
of b and T. The length of the periodic cell L ~ 7R9 in each case. The statistical
uncertainties in the diffusion coefficient are smaller than the size of the symbols,
except for the longest chain, N = 1024, where they are shown explicitly . 56

3-10 Center of mass diffusion coefficients parallel (Dll) and perpendicular (D') to
the channel walls as a function of channel width H; the radius of gyration of
the polymer (N = 128, b = 1) is approximately 7Ax. The confined diffusion
coefficients, normalized by the monomer diffusivity, are shown for a square cell,
L, with various ratios of L to H . . . . . . 58

3-11 Center of mass diffusion coefficient for a confined polymer relative to the unconfined
diffusivity, DH /D. The solid line indicates the Faxen correction [7] for weak confinement,
assuming the polymer lies in the center of the channel. The dashed line is the
theoretical (R,/H)-2/3 scaling [46]. The aspect ratio of the cell, L/H, was set
to 4 in each case . . . . . . . . 60

4-1 Center of mass distribution for pressure driven flow as a function of Peclet number
and channel width. Results are for N 16. The distributions are normalized
such that o HRg P(y/R,)d(y/R,) 1. Due to the symmetry of the system, the
distribution is plotted up to 0.5H. The bounding wall is at the origin. . 66

4-2 Center of mass distribution for uniform shear flow as a function of Peclet number
and channel width. Results are shown for N 16. The distributions are normalized
such that H/Rg P(y/R)d(y/R) . . . . . ..... 67

4-3 Components of the radius of gyration along the flow, (x2), and confinement,
(y2), axes. Results for chains of length N 16 are shown at two different Peclet
numbers, normalized by R ((x2) + y2) + (2)05. . . . ..... 68

4-4 Center of mass distribution for different levels of discretization of the polymer
in a tightly confined channel, H/R, 3. The distributions are normalized such
that fo P(y/R )d(y/R ) 1. . . . . . . ... 69

4-5 Average axial velocity of the center of mass of polymers normalized by the average
fluid flow. The results are presented with respect to the polymer Peclet number
for different channel sizes . . . . . . . 70

4-6 Hydrodynamic dispersion coefficients of a polymer chain in Poiseuille flow. Numerical
simulations are shown as filled symbols, along with Taylor-Aris theory (TA) for
tracer particles (Ld = 0, dashed line) and for finite sized particles with a fixed
depletion l v.r (Ld = Rg = H/8, solid line). Theoretical calculations, accounting
for the depletion l-1 vr thickness determined from simulation, are shown as open
symbols. The standard error in the simulation data corresponds to the size of
the sym bols . . . . . . . . . 72

4-7 Comparison of dispersion coefficients with different inertia. The highest particle
Reynolds number for kT =0.1 is almost 3.4, whereas the channel Reynolds numbers
are up to 30 for the same simulations. . . . . . 73

4-8 Comparison of center of mass distributions at 3 different temperatures at Pe =
300. The three different temperatures correspond to different Reynolds numbers.
The highest particle Reynolds number at kT =0.1 is around 3.4. . . 74

5-1 Center of mass distribution of a single polymer in a channel H/Rg = 8 under
the application of an external body force. Half of the distribution is shown, with
the boundary at y/H = 0 and the center of the channel at y/H = 0.5. .. . 77

5-2 Illustration of different migration mechanisms. In each figure the solid lines represent
the forces while the dashed lines represent the velocities generated by these forces.
The velocities can be calculated using an image system for a point source near
a planar boundary. A) The lift due to shear: the tension in the polymer near
a solid boundary generates a net velocity .1-iv-, from the boundary. B) Rotation
due to the external field: the velocity field due to the external force on each bead
results in a rotation about the center of mass of the polymer. C) Drift of a rotated
polymer: two beads oriented at an angle to the external forces drift due to the
hydrodynamic interactions between the beads. . . . . 79

5-3 a) Scaling of the distribution function with respect to the channel width under
constant polymer Peclet number, Pe. b) Scaling of the distribution function
with respect to polymer Peclet number using different levels of chain discretization.
The boundary is at y/H = 0 and the center of the channel is at y/H = 0.5. 81

5-4 Center of mass distributions for concurrent application of an external body force
and pressure-driven flow. The solid curve shows the level of migration under
the pressure-driven flow only. The flow Peclet number in all three cases is Pef =
12.5. The boundary is at y/H = 0 and the center of the channel is at y/H = 0.5. 82

5-5 Center of mass distributions for countercurrent application of an external body
force and pressure-driven flow. The solid curve shows the level of migration under
the pressure-driven flow only. The flow Peclet number in all cases is Pe = 12.5.
The boundary is at y/H = 0 and the center of the channel is at y/H = 0.5. 83

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy



0. Berk USTA

May 2007

C!i ,i': Jason E. Butler
CoC' ki': Anthony J. Ladd
Major Department: CI, i ii' Engineering

We investigate the dynamic and static properties of polymeric suspensions via

a numerical simulation method. This method couples Newtonian mechanics of solid

phase suspended particles with a fluctuating lattice-Boltzmann fluid. This results in the

solution of the Navier-Stokes equation, including thermal fluctuations and hydrodynamic

interactions for the solid particles. The coupling can be achieved either through boundary

conditions for a finite sized particle or through a Stokes-like frictional coupling for point

particles. The frictional coupling procedure is a very efficient approach due to the use of

point particles which take up much smaller volumes. This in turn results in simulations

of far fewer lattice nodes, while retaining essential hydrodynamic features for dynamics of

polymers. This approach has enabled us to simulate polymer chains with more than 1000

segments, including hydrodynamic interactions. This discretization is nearly an order of

magnitude higher than what is feasible using traditional Brownian dynamics methods.

The fluctuating lattice-Boltzmann method is presented for both colloidal and

polymeric suspensions. A verification of the fluctuations for the case of colloidal finite

sized particles is presented as a test case. However, the main theme of this thesis is the

simulations of dilute polymer solutions.

The hydrodynamic interactions p1 i, an important role in these systems and the

effect of these interactions is the focus of this study. In unbounded geometries, the

effects lead to important differences in the scaling of the polymer properties, such as

the radius of gyration and diffusion coefficient, with the linear length of the polymer. In

bounded geometries the hydrodynamic interactions can result in a transverse migration.

In unidirectional shearing flow or external fields, this leads to a nonuniform concentration

of the polymer across the lateral axis. In both cases although the underlying phenomenon

is the hydrodynamic interactions, the motion created by these two driving forces are

distinctly different.

The magnitude and direction of the migration can be manipulated using a combination

of hydrodynamic and external fields. The migration observed in this study also affects

macroscopic transport properties of the polymer, such as dispersion and segregation

velocities, due to changes in the center of mass distribution. The migration under all flow

conditions depends on the dimensionless ratio between the channel size and the polymer

size. Therefore, separation of polymers based on chain length is possible as long as the

shear rate across the channel varies. Our results for combined flows, when compared with

recent experiments, also indicate that hydrodynamic interactions, usually neglected in the

electric field driven motion of polymers, may only be partially screened.


1.1 Introduction

Colloidal and polymeric suspensions are found in a great v ,ii I v of natural and

artificial materials and processes. These range from the foods we consume every w to

biological systems. Therefore understanding both the static and dynamic properties of

such systems is imperative. Such an understanding enables us to successfully use these

materials and create new ones. An understanding of natural processes involving such

systems can also lead to mimicking these in artificial devices.

Colloidal and polymeric systems are often classified as "complex fluids" owing to the

nonlinear relationship between the applied forces and resulting dynamics. This complexity

arises from the wide range of time and length scales that govern their motion. The motion

of these systems depend on a variety of different factors such as particle interactions

and large scale collective motion related to structural relaxations. The simplest class of

interactions are the direct inter-particle interactions such as electrostatic interactions and

hard-sphere collisions which also exist in simple fluids. In addition to these, sub-micron

solute particles in a colloidal suspension are coupled in their motion through the fluid.

This coupling is called the hydrodynamic interaction [1] and constitutes the focus of this

work. These hydrodynamic interactions are long ranged and ii ,niv-1 body, therefore they

create a great deal of complexity for theoretical calculations.

In this work we investigate the static and dynamic transport properties of hard-sphere

colloidal suspensions and dilute polymer solutions with a numerical simulation method.

This method combines Newtonian dynamics of solid particles with a fluctuating

lattice-Boltzmann fluid [2-4]. Such an approach leads to the solution of N ',.i. -Stokes

equations, with suspended particles, on large time and length scales. This method is

closely related to and is a descendent of the lattice-gas cellular automata for colloidal

suspensions [5].

In C'!i ipter 1, I will briefly discuss the hydrodynamic equations governing the motion

of the suspending fluid (1.2) along with a discussion of the lowest order description of

hydrodynamic interactions between suspended solid particles (1.3). I will also present a

brief history and summary of numerical techniques on the subject.

In C'!i Ipter 2, I will present the fluctuating lattice-Boltzmann method. This includes

two different approaches for the solid particles coupled to the lattice-Boltzmann fluid for

hard-sphere colloidal and polymer suspensions. I will then present static and dynamic

problems (C'!i ipter 3) involving dilute polymer solutions. This is followed by a numerical

investigation of the migration of polymer chains under unidirectional shearing flows and

external fields in C'!i Ipters 4 and 5.

1.2 Fluid Motion and Hydrodynamic Equations

Suspensions are two phase materials usually consisting of a solid phase and a

suspending fluid. Therefore it is necessary to understand the motion of both phases,

the fluid and the solid, and also their coupling. In this section I briefly describe the mass

and momentum conservation equations in such suspensions which are better known as

the continuity and Navier-Stokes equations. I then describe the coupling between the

two phases through boundaries and the hydrodynamic interactions that arise due to this

coupling. In this study we assume that the suspending fluid is a Newtonian liquid.

The macroscopic equations describing the motion of a fluid is based on the continuum

hypothesis and a mass and momentum balance on a finite volume element. This volume

is small enough for the properties to be local but large enough to contain a high number

of molecules such that fluctuations arising from from the different properties of these

molecules have no effect on the observed average [6]. In physical terms this requires the

ratio of mean free path (A) of a molecule to the linear dimension (L) of the volume to be

small such that A/L < 1. This ratio is known as the Knudsen number.

The conservation of mass applied to this volume element gives the continuity

Op ( )
=- --V.(pv), (1-1)

where p is the local density of the fluid and v is the local mass average fluid velocity. In

this study we will assume that the fluid is incompressible, in which case Eq. 1-1 reduces


V v 0. (1-2)

The conservation of linear momentum for this finite volume element can be obtained

by applying Newton's laws. Newton's laws may be interpreted as stating that the external

forces exerted by the the surroundings on a stationary fluid element are equal to the time

rate at which which momentum is being created within the element. The external forces

are contact forces exerted by the fluid stresses, acting over the surface of the element, and

the volume body forces such as gravity. The rate of creation of momentum is equal to the

accumulation plus the rate of out-flux of momentum through its surface. [7]. Therefore the

momentum balance equation reads,

a + V (pvv) V H + pF, (1-3)

where H is the stress tensor and F is the external force. For a Newtonian fluid the stress

tensor is given by,

H -pl + K V v + 2 (VV)'. (1-4)

Here p is the hydrostatic pressure at rest, I is the unit tensor, p is the viscosity, and

(VV)8 is the rate of deformation tensor [7]. The bulk viscosity, K, can be ignored for a low

density monatomic gas [7] and an incompressible fluid. Substitution of this information

and assuming incompressibility in Eq. 1-3 gives the Navier-Stokes equation:

p (O + v Vv -Vp + /V2v + pF. (1-5)

When the viscous terms, pV2v, are much larger than the inertial terms, pv Vv, the

latter can be neglected. Formally the ratio of inertial to viscous forces is known as the

Reynolds number, Re = where L is a characteristic macroscopic length. For small

enough Reynolds numbers (Re < 1), the inertial terms can be neglected, and with the

continuity equation this gives rise to the creeping motion or Stokes equations of fluid flow


pV2v = Vp (1-6)

V v 0. (1-7)

The Stokes equations constitutes a great simplification over Navier-Stokes equations for

the solution of fluid motion, where suitable. In this study we will only consider creeping

motion for the fluid unless otherwise noted.

1.3 Hydrodynamic Interactions

The motion of a particle suspended in a fluid is coupled to the motions of the

fluid through momentum transfer between the two phases. In turn this coupling also

gives rise to an interaction between suspended particles that is mediated through the

suspending fluid. In this section I will briefly discuss the coupling of the two phases and

the hydrodynamic interactions.

Solid particles are rigid bodies which do not deform, in contrast to fluids. The rigidity

constraint dictates that under pure translation, every point on the particles has the same

velocity. The particle may also undergo rotation when different points in the particle have

different velocities. The velocity of any point U(r) is given by [8],

U(r)- U + Qx (r- R), (1-8)

where U is the center of mass translational velocity, Qf is rotational velocity, and r R is

the position vector of a point relative to the center of mass vector of the particle.

No-slip (stick) boundary conditions apply at the solid-fluid surface when a rigid

particle is immersed in a fluid. The fluid velocity field on the surface of a particle under

this condition, following Eq. 1-8, is

v(r) -U + x (r- R) for r R = a, (1-9)

where a is the particle radius [9]. If the relative velocity (U v(r)) of the solid and fluid

phases is small, the effects of the fluid can be described in terms of forces and torques

exerted on the solid particles [9]. The forces, F, and torques, T, can be expressed in terms

of the particle velocities, U, and f2, through configuration dependent friction tensors (,
Fi -- U + Q] (1-10)
[13 Uj i (1310)
Ti I j + Ci Qj]. (1-11)

Here the superscripts T and R refer to translational and rotational components of the

friction tensors and the summations extend over all particles in the suspension. These

equations can be inverted to give mobility matrices p,
Ui [ TFj + F Tj] (1-12)
i -- [i Fj + Tj]. (1-13)

The mobility matrices are often termed "the hydrodynamic interaction" tensors since

they represent the interactions between suspended particles through the fluid. These

hydrodynamic interactions can be calculated using the method of induced forces [10]

where a distribution of point forces on the surface of each sphere is assumed and this

distribution is then expanded about the center of each sphere in a series of spherical

harmonics [5]. These interactions are long-ranged, not pairwise additive, and often

impossible to calculate analytically except for special cases such as pairs of suspended

particles. The difficulty arises due to the divergence of the integral from the effect of

the induced forces which create a velocity field decaying as r-1 and from the n iiv-1 body

nature of the interactions.

When the particles in the suspension are small enough such that they are much

smaller than the average inter-particle distance, these particles can be represented by point

particles. On these scales, the rotational degrees of freedom and microscopic details of

solid fluid interaction such as the stick boundary conditions can be neglected. This allows

for a simple frictional coupling between the solid particles and the fluid [11]. This is a

great simplification and such conditions can usually be assumed for a polymer chain that

consists of hydrodynamic centers (beads) connected by springs. The interaction of such

particles can be described locally in terms of the relative velocity between the fluid and

the solid particle and the collective influence of particles becomes additive in the creeping

flow limit. Each particle in the suspension induces a velocity in the fluid that follows

v(r) H f (1-14)

where H is the Oseen tensor. This tensor is given by

H I + ]. (1-15)
7STrr r2

1.4 Numerical Methods to Solve Hydrodynamic Equations

Numerical techniques, over the past few decades, have proven to be invaluable tools

to study the dynamics of colloidal and polymeric suspensions. Their use has led to the

discovery of important physical phenomena such as long time tails [12] and difference

in the shear-rate dependence of the hydrodynamic and Brownian components of stress

[13]. Computer simulations can resolve problems which are not accessible by laboratory

experiments. Simulations allow for investigation of ideal systems which help identify subtle

individual effects, one of which is hydrodynamic interactions; the focus of this thesis. In

this section, I list and briefly review a few important simulation methods for colloidal and

polymeric suspensions

Molecular dynamics methods are the most simple, straightforward, and versatile

simulation method for a wide variety of different problems including colloidal and

polymeric suspensions. If the rotational motion of particles in a system is negligible,

Newton's equations of motion are applicable for each particle i:

d2 r
m 2 Fi, (1-16)

where m is the mass of the particle, r is the position vector, and F is the force on the

particle. These equations of motion are then integrated for each particle in the system.

Hydrodynamic interactions arise naturally due to the collision events between the

solid and fluid particles. Transport properties such as viscosity and diffusivity can be

calculated by using time-correlation formalism [14]. Although quite straightforward, such

an algorithm is computationally very ineffective when one is interested in the large scale

dynamics of colloidal particles which are bigger than the surrounding fluid particles. The

simulation of smaller fluid particles takes up most of the computational time, therefore a

very limited time scale can be simulated with such an approach.

The size separation of suspended colloidal particles and the fluid also results in a time

scale separation for the motion of the two phases. One can exploit these scale separations

and describe the fluid as a continuum which results in a huge reduction of degrees of

freedom. Stokesian dynamics and Brownian dynamics are two such methods where

hydrodynamic interactions are assumed to be instantaneous due to the above mentioned

time-scale separation.

Stokesian dynamics is an accurate simulation technique which uses a reformulation of

the Stokes equations (Eq. 1-6) as an integral equation for the fluid velocity field generated

by a prescribed force density find [2]. The force density is localized to the particle surfaces

and for spherical particles one needs to expand the force density and resulting velocity

field in spherical harmonics. This leads to a linear system of equations relating moments

of the fluid velocity on the particle surfaces to moments of the induced force fid. The

exact solution of this problem requires an infinite hierarchy of multipoles, but Stokesian

dynamics methods usually truncate this at L = 1 order which includes torques and stresses

on a particle. In addition to this exact solution of the far field interaction between distant

particles, Stokesian dynamics also incorporates lubrication forces when particles are in

near contact. This interaction is incorporated in an additive manner although in general

hydrodynamic interactions are not additive.

Although Stokesian dynamics has been a successful simulation method it has some

limitations besides the computational cost which usually scales as O(N3); this could be

improved by fast-multipole methods to O(N2). The suspending fluid is assumed to be

Newtonian and only Stokes flows of (Re = 0) can be simulated. The method is difficult

to extend to non-spherical particles or to suspensions with bounded geometries which

requires recalculation of the Green's function for each different geometry.

Brownian dynamics simulations are based on a stochastic reformulation of the

equations of motion of solid particles. They are particularly well suited to studying the

structure and rheology of complex fluids in hydrodynamic flows and other nonequilibrium

situations. Explicit solvent molecules are represented by a stochastic force since there

is a large separation in time scales between the rapid motion of solvent molecules and

the long-time scales for the motion of polymers or colloidal particles. Therefore one can

assume, neglecting inertia, that all forces on a particle are balanced at each instant:

FH + F +FNH 0, (1-17)

where H denotes hydrodynamic forces, B Brownian forces, and NH non-hydrodynamic

forces such as external and potential forces. The hydrodynamic forces can be written as a

drag between the particle and fluid phase based on the Stokes friction:


where ( is the friction coefficient and ur(ri) is the unperturbed fluid velocity at the

position ri of the particle. Equation 1-18 can be rewritten, with the substitution for the

hydrodynamic term, which yields a Langevin type equation

dr 1(19
= 4(ri) + (F + FNH ()),19)

where FB(t) is a stochastic force and

< FB >- 0 and

< F (t)Ff(t') > 2kTU6iy6(t- t').

This equation can be integrated forward in time to create trajectories of molecules. In

simpler implementations of Brownian dynamics, hydrodynamic interactions are usually

neglected. However, it is possible to incorporate hydrodynamics through the use of an

interaction tensor which represents the correlated forces on each particle due to motion

of all other particles in the suspension. The interaction tensor is configuration-dependent

and the calculation for the Brownian contribution requires decomposition of this tensor.

This calculation is computationally intensive and, similar to Stokesian Dynamics, usually

scales as the square or cube of the number of particles in the system, which puts Brownian

Dynamics methods at a disadvantage for large systems. The interaction tensor needs to be

modified for each different geometry which becomes cumbersome and is inflexible.

In contrast to Stokesian dynamics and Brownian dynamics methods described above,

time dependent methods, such as finite difference solutions of Navier-Stokes equations,

lattice-gas and lattice Boltzmann methods, rely on an explicit solution of the fluid motion.

A solid phase can then be coupled to fluid through appropriate boundary conditions.

This leads to a temporal evolution of hydrodynamic interactions unlike the assumption

of instantaneous interactions in the traditional methods. This natural evolution of

hydrodynamics avoids the need for calculation of complicated covariance matrices for the

correlated motion of the particles. As a result the computational effort usually depends

on the fluid volume rather than the number of particles. The temporal evolution not only

allows for a reduction in computational cost but also allows for probing of shorter time

scales than is possible with Brownian dynamics, which only probes beyond diffusive time


In this study we have chosen to study a coupling of a lattice-Boltzmann fluid with

different descriptions of the solid phase. The method will be described in Chapter 2.

Verification tests along with physical problems will be presented in Ch'! pters 3, 4, and 5.


The lattice-Boltzmann method is a basis for simulating solid-fluid suspensions. The

method uses a discretized Boltzmann equation for the evolution of the fluid phase, and

the solid particles are coupled to the fluid through appropriate boundary conditions and

integration of Newtonian dynamics [4]. The most important features of this method

are the natural evolution of hydrodynamic interactions and a linear scaling of the

computational cost with respect to the number of particles. However, it is important

to note that the computational cost in fact scales directly with the fluid volume.

This C!i Ipter includes a brief review of the basic concepts of the lattice-Boltzmann

method, starting with a discrete Boltzmann equation and description of the fields and the

lattice. The incorporation of fluctuations and the coupling of solid particles to the fluid

are also discussed. Readers interested in a more detailed discussion should refer to papers

by Ladd and his group [2, 4, 15 17].

2.1 Boltzmann Equation

The Boltzmann equation is one of the greatest achievements of theoretical physics

and a milestone in kinetic theory [18]. In this section I will start by describing some key

points of this equation. The probability density or distribution function, f(x, p, ), gives

the probability of finding a molecule around position x with momentum p at time t. This

is the key object of the kinetic theory. The Boltzmann equation, derived in 1872, describes

the evolution of this distribution function, f, in terms of microdynamic interactions. This

equation reads as follows.

Dtf [Ot + P- + F FOp]f(x,p,t) 12, (2-1)

where the left hand side represents the streaming motion of the molecules along

the trajectories associated with the force field F and C12 represents the effect of

intermolecular(two-body) collisions taking molecules in and out of the streaming trajectory


The collision operator involves the two body distribution function f12 expressing the

probability of finding a molecule 1 around x, with velocity vi and a second molecule 2

around x2 with velocity v2, both at time t. The difficulty in writing f12 is that it calls into

pl i, the three-body distribution function f123 which in turn depends on f1234 and so on

according to the BBGKY(Bogoliubov, Born, Green, Kirkwood, Yvon) hierarchy [19]. In

order to achieve a closure in his equation Boltzmann made a few stringent assumptions

on the nature of the system (2- 1); a dilute gas of point-like, structureless molecules

interacting via a short-ri,.- two-body potential.

Under such conditions intermolecular interactions can be described solely, in terms

of localized binary collisions, with molecules spending most of their lifespan in free

trajectories (in the absence of external fields) perfectly unaware of each other. Within this

approximation we can split the collision term into gain and loss components:

12 G L (fl'2' f12)go-(g, Q)dQdp2 (2-2)

corresponding to direct and inverse collisions taking molecules in and out of the volume

element dpl. The link to underlying molecular dynamics is channeled into the differential

cross section a expressing the number of molecules with relative speed g = gQ around the

solid angle Q2. This is followed by the famous closure approximation made by Boltzmann:

f12 f2 (2-3)

The closure relies on the idea of molecular chaos (or Stosszahlansatz), which is fairly

plausible for a dilute gas with short range interactions in which molecules spend most of

the time traveling in free space only to meet occasionally for very-short lived interactions.

We should also note that the molecules are strongly correlated after the collision by virtue

of the conserved quantities, mass, momentum and energy.

The probability for two molecules that met at time t 7 to meet again at time t with

the same moment pi and P2 is exponentially small with T. More precisely this scales

like e-/ i1* where Tint is the duration of a collision event. In Boltzmann's theory Tint (this

scales as s/v, where v is a typical particle speed and s is the effective molecular diameter)

is negligibly small [18]. In a liquid the situation is completely different due to much higher

density. In this case the particles are in constant interaction with one another. Violations

of Boltzmann's molecular chaos can occur in relatively dilute fluids due to nonlinear

fluid-particle interactions. One important example is the long-time tails first detected by

Alder and Wainwright [12] .

To summarize all the above information, we write the Boltzmann equation in its final

form [18].

[Lt + P- + F p] f(x, p, t) = (fl'2' f12)g Q)dQdp2. (2-4)

The left hand side is similar to Newtonian single-particle dynamics, while the right hand

side describes intermolecular interactions under the molecular chaos approximation. The

most important feature of the equation is that it represents the local conservation laws

which are of great use in deriving the continuum result, the N i,. r-Stokes equations.

2.2 Lattice Boltzmann Method

2.2.1 An Overview of Lattice-Boltzmann Method

Historically Lattice Boltzmann is the descendent of the Lattice Gas Automaton

(LGA) where the discrete positions and the velocities of the single particles are tracked on

the lattice while satisfying the conservation equations locally. Lattice gas methods proved

to be quite simple and useful for hydrodynamics, though had some drawbacks such as the

statistical noise which requires long averaging times.

The lattice Boltzmann methods, though they can be independently derived from the

Boltzmann equation, evolved from the idea of fixing the drawbacks of the LGA approach.

Lattice Boltzmann follows the evolution of the distribution function for the particles

rather than tracking individual fluid particles, which results in pre-averaging in momentum

space. This approach eliminates the statistical noise together with the thermal noise. In

order to simulate fluctuating phenomena such as the Brownian motion of solid particles,

we introduce the thermal fluctuations consistent with fluctuating hydrodynamics [20]. This

is achieved by introducing a discrete random stress tensor in agreement with a discrete

fluctuation-dissipation theorem.

The computational method I present in this Ch'! pter is based on a fluctuating

lattice-Boltzmann model of the fluid phase, [3, 4, 21] which includes thermally driven

fluctuations in the fluid velocity field via random fluctuations added to the stress

tensor. The key difference with Brownian dynamics is that these random fluctuations

are local in space and time, avoiding any factorization of the covariance matrix. Instead,

initially uncorrelated fluctuations develop in space and time by viscous diffusion of

momentum, giving rise to hydrodynamic correlations at sufficiently large spatial and

temporal scales. Numerical simulations have shown that the dissipative and fluctuating

hydrodynamic interactions between finite-size particles can be reproduced quantitatively

[4]. Computationally, the lattice-Boltzmann method is inherently and straightforwardly

parallelizable, enabling large simulations to be spread across a cluster of processors with a

relatively low network bandwidth.

The lattice-Boltzmann method has been used for a wide variety of complex flows,

including suspensions of colloidal [3, 22] and non-colloidal [23-25] particles, and porous

media [26-28]. Since the fluid equations are solved on a grid, external boundaries and

imposed flow fields can be included at no additional computational cost. However,

suspended solid particles occupy a substantial fluid volume, typically of the order of 100

grid points, and interact with the fluid through boundary nodes distributed on the particle

surface. Since the computation time of the method is proportional to the fluid volume, the

suspension model is unnecessarily expensive for polymer solutions, where the monomers

can be modeled as point forces coupled to the fluid by a friction coefficient, Ahlrichs

and Diinweg integrated this frictional coupling into the fluctuating lattice-Boltzmann

model and investigated the static and dynamical scaling laws of individual polymers and

semi-dilute solutions [11, 29].

2.2.2 Lattice Boltzmann Equation

The fundamental quantity in the lattice-Boltzmann model is the discretized

one-particle velocity distribution function ni(r, t), which describes the mass density of

particles with velocity ci at a particular node of the lattice r at a discrete time t. The

hydrodynamic fields, mass density p, momentum density j = pu, and momentum flux H,

are moments of this velocity distribution:

P EYi j E i nici, H ncici. (2-5)

The time evolution of ni(r, t) is described by a discrete analogue of the Boltzmann

equation, [30]

ni(r + ciAt, t + At) = ni(r, t) + Ai [n(r, t)], (2-6)

where Ai is the change in ni due to instantaneous collisions at the lattice nodes and At

is the time step. Somewhat surprisingly, this simple evolution equation is second-order

accurate in space and time. The numerical diffusion that usually accompanies a low-order

grid-based method is eliminated by the relationship between the eigenvalues of the

linearized collision operator and the fluid viscosity.

The standard 19 velocity model comprises stationary particles and the 18 velocities

corresponding to the [100] and [110] directions of a simple cubic lattice. The population

density associated with each velocity has a weight ac' that describes the fraction of

particles with velocity ci in a system at rest; these weights depend only on the speed ci

and are normalized so that L ac = 1. The optimum choice of weights for this model is

ao 1 a = 1 (2-7)
3' 18' 36 (27)

In these simulations we use a 3-parameter collision operator, allowing for separate

relaxation of the 5 shear modes, 1 bulk mode, and 9 kinetic modes. The post-collision

distribution nt = nt + A, is written as

nt (= ac/ C 2 (2-8)

where the sound speed c = Ax/v3At, Ax is the spacing between neighboring fluid
nodes and At is the timestep. In suspensions of moving particles, the non-linear puu
term in Eq. 2-8 should be retained, since it maintains Galilean invariance and prevents an
artificial cross-stream drift [21]. Empirically, we have found that the domain of validity of
Stokes flow is considerably larger with the non-linear term in place.
The non-equilibrium momentum flux fIneq neqC C. relaxes due to collisions at
the lattice nodes,
Fneq,* (1 + A) f + -(1 + A,)(fI" : 1)1, (2-9)
where Fnueq I I eq, IPq p + puu is the equilibrium momentum flux and Ineq
indicates the traceless part of fleq. The parameters A and A, are eigenvalues of the
linearized collision operator and are related to the fluid shear and bulk viscosities:

p pAt + _- AAt ( + (2-10)

The factor of 1/2 corrects for numerical diffusion, so that viscous momentum diffuses at
the expected speed for the given viscosity.
In the presence of an externally imposed force density f, for example a pressure
gradient or a gravitational field, the time evolution of the lattice-Boltzmann model
includes an additional contribution f/(r, t),

rt(r + cAt, t + At) -r (r, t) + f,(r, t). (2- 11)

This forcing term can also be expanded in a power series in the velocity,

f c (uf + fu) : (cc, cc 1)
fi arc + 2cc At. (2-12)
cS 2 SI

More accurate solutions to the velocity field are obtained if it includes a portion of the

momentum added to each node, [21]

j' pu' = nic, + f At. (2-13)

The macrodynamical behavior arising from the lattice-Boltzmann equation can be

found from a multi-scale cin i ,'-i- [30] using an expansion parameter c, defined as the

ratio of the lattice spacing to a characteristic macroscopic length; the hydrodynamic limit

corresponds to c
the Navier-Stokes equations with corrections that are of the order u2 and c2. Thus, at

sufficiently low Mach numbers, the method is second-order accurate in space with relative

errors proportional to Ax2. It is also second-order accurate in time if the viscosities are

defined according to Eq. 2-10.

2.3 Fluctuations

The lattice-Boltzmann model can be extended to include thermal fluctuations,

which leads to Brownian motion of colloidal particles and polymers. The fluctuating

lattice-Boltzmann model [3, 4] incorporates a random component into the momentum flux

during the collision process (c.f. Eq. 2-9):

neq, (1 + A)eq+q + }(1 + )(fleq 1)1 +If
ni ( + (R 1, (2-14)

where R, is a random variable with zero mean and unit variance, and R is a symmetric

matrix of random variables of zero mean. The off-diagonal elements of R have a variance

of 1, while the diagonal elements have a variance of 2. In this particular implementation,

6 random numbers are required to generate the components of the symmetric matrix R,

together with the constraint that R is traceless. Random numbers spanning a finite space

are advantageous since the change in any one step is then limited. We use integer random

numbers [-2, -1, 0, 1, 2] with weights [1/12, 1/6, 1/2, 1/6, 1/12] to obtain a unit variance
and second-order accuracy in time.
An analogy with fluctuating hydrodynamics [20] sI-- -1- that

( 2T 1/2 ( 2T2, 1/2
AXAt, AXAt, (2-15)

where T is the thermal energy in units of moAx2/At2 and mo 0 pAx3/36 is taken to be
the rest mass of an individual population density. However, detailed calculations show that
discrete lattice artifacts again modify the results for continuous fluids [4, 21]. The correct
coefficients, determined by relating the decay of long wavelength stress fluctuations to the
viscosity of the lattice-Boltzmann fluid are [4, 21]

2TA2 1/2 2TA 1/2 ( 16)

There is an additional multiplicative factor of A-2 (c.f. Eq. 2-15), but for A = A, = -1,
an exact correspondence with the fluctuation-dissipation relation for continuous systems
is obtained. However the discrete version (Eq. 2-16) can also be applied to systems where
the viscous eigenvalues are not equal to -1 [4]. This may be advantageous in problems
where a time-scale separation is required between hydrodynamic and kinetic modes.
2.4 Solid Particles
2.4.1 Finite-Size Spherical Particles
In our implementation of the Lattice-Boltzmann model, the solid particles are defined
by a surface which cuts some of the links between lattice nodes (Fig. 2-1). Fluid particles
propagating along the links interact with the solid surface at boundary nodes placed
halfway along the links. This provides a discrete representation of the particle surface
which becomes more precise as the particle dimension to grid size increases. A detailed
explanation of the treatment of solid-fluid boundaries can be found in references [2, 4, 16].
Boundary Conditions.

Figure 2-1. Location of the boundary nodes for a circular object of radius 2.5 lattice
spacings. The velocities along links cutting the boundary surface are indicated
by arrows. The locations of the boundary nodes are shown by solid squares
and the lattice nodes by solid circles.

A moving boundary condition, without any interior fluid is implemented as follows.

We take the set of fluid nodes just outside the particle surface, and for each node all the

velocities Cb such that r + cbAt lies inside the particle surface. An example of a set of

boundary nodes is shown in Fig. 2-1. Each of the corresponding population densities

is then updated according to a simple rule which takes into account the motion of the

particle surface. [4]

n(r, t + At) (r, t) (2-17)
where n (r, t) is the post-collision distribution at (r,t) in the direction cb and c -cb.

The local velocity of the particle surface,

Ub =U + QX (rb R) (2-18)

is determined from the particle center-of-mass velocity U, angular velocity fl, and

center-of- mass R; rb r + oCbAt is the location of the boundary node. We use the

mean density po in Eq.2-17 instead of the local density p(r, t) since it simplifies the update

procedure. The differences between the two methods are small, of the same order (pu3)

as the error terms in the lattice-Boltzmann model. Test calculations show that even large

variations in fluid density (up to 10' .) have an effect on the force less than 1 part in 104.

As a result of the boundary node updates, the momentum is locally exchanged

between the fluid and the solid particle. This results in the conservation of the total

momentum both locally and globally. The forces exerted at the boundary nodes can be

calculated from the momentum transferred in Eq.2-17,

f(r, t + At) 2n (r, t) 2abpOUb (2-19)
2 At cS
The particle forces and torques are then obtained by summing f(rb) and rb x f(rb)

over all the boundary nodes associated with a particular particle. It can be shown

analytically that the force on a planar wall in a linear shear flow is exact and examples of

lattice-Boltzmann simulations of hydrodynamic interactions are given in Ref. [16].

2.4.2 Point Particles

A polymer chain is represented by a bead-spring model with an equilibrium bond

length b between neighboring beads. The spring stiffness, K = 30T/b2, was chosen so that

the fluctuations in bond length, ((r b)2)/2 b/ /30, were small in comparison with the

radius of gyration of the chain, R.. A hard-sphere excluded-volume interaction between

the monomers was included, with collision diameter b/2. Other excluded volume and

intra-chain potentials could be used.

The equations of motion for the monomers are written in inertial form,

d2X, ViIJXN) + F, (2 20)

where KP is the total potential energy and F{ is the hydrodynamic force exerted on bead

i by the fluid. The Brownian dynamics algorithm is derived from Eq. 2-20 by neglecting

the inertial term and integrating the random component of the hydrodynamic force over

the duration of the time step [31]. In this thesis evidence will be presented to support

our contention that it is more efficient to integrate Eq. 2-20 directly. The large time-scale

separation between the dynamics of the polymer and the individual monomers allows

time for the hydrodynamic interactions to reach a quasi-steady state, without imposing

this condition at each and every time step. We will show that both inertial and diffusive

simulations can use similar time steps, of the order of the monomer diffusion time. This

allows for orders of magnitude savings in computation time by keeping the interactions

local in space. We note that the lattice-Boltzmann model is just one possible framework

for such calculations, which could also be carried out within a fluctuating finite-difference

or finite-element simulation of the fluid.

The hydrodynamic interactions between the beads are determined by coupling them

individually to the fluctuating lattice-Boltzmann model, from which the mesoscopic

equations for a fluid with thermal fluctuations can be derived [21]. The coupling between

the beads and the fluid is derived from a frictional force based on the difference in velocity

between the bead, U, and the surrounding fluid, u(r, t), [11]

Ff -o [U(t) u(X, t)] + Fr. (2-21)

The random force F' is introduced to balance the additional dissipation caused by using a

frictional coupling instead of a no-slip boundary condition on the bead surfaces [11]. Since

the fluid satisfies its own fluctuation-dissipation relation, F' has a local covariance matrix

(Fr(t)Fr(t')) 2ToJ6(t t')1. (2-22)

Hydrodynamic interactions between the beads are transmitted through the fluid via

correlated fluctuations in the fluid velocity field which develop over the inertial time scale,

pr2/2, where r is the separation between beads.
The point-particle approach leads to an Oseen interaction between individual

monomers as will be shown numerically in Sec. 3.2.2. The time scale for hydrodynamic

interactions to develop, pR21/T, must be much shorter than the diffusion time for the

polymer, RJ/D. In other words, viscous momentum must diffuse considerably faster than

the timescale for appreciable changes in polymer configuration. This scale separation is

characterized by the Schmidt number Sc = ]/pD, which for the temperatures and bond

lengths used in the simulations is in the range 3N' 60N'. This scale separation is larger

and more easily controlled than in Dissipative Particle Dynamics simulations [32] and the

effect of variations in Schmidt number is usually small. The Euler method was used to

integrate the equations of motion of the monomers with a time step that was typically

1/100 of the monomer diffusion time. When necessary, the thermodynamic forces were

integrated with a smaller time step than the hydrodynamic forces to maintain stability. A

hard-sphere molecular dynamics code was used at each step to detect collisions between

pairs of monomers.

Since the monomers move continuously over the grid, while the velocity field is

evaluated at a discrete set of points, we employ a trilinear interpolation to evaluate

u(X, t), using the fluid velocities at the nodes of the cube surrounding the monomer

[11]. After calculating the weights, ,, for each of the nodes, the mass density, p, and the

momentum density pu are interpolated to determine the velocity field at the bead location


u(X, t) = (2-23)

To conserve momentum, the accumulated force exerted by the bead on the fluid

is distributed to the surrounding nodes with the same weight function. Trilinear

interpolation is computationally efficient and works satisfactorily as demonstrated

by numerical tests. However, there are discontinuities in the spatial derivative of the

velocity field at the cell boundaries, which can be greatly reduced by using second order

interpolation. Nevertheless, the increased accuracy does not clearly justify the additional

computational cost of employing a higher order interpolation. Consequently, the results

in this thesis use linear interpolation, though higher order interpolation schemes will be

considered in future work.

2.4.3 Particle Motion

Three different schemes for the translational velocity update of the particles is

performed with an adjustable parameter a as follows. The equation for the velocity

update is

U(t + At) U(t) + F At,
U M)-t


where U is the particle velocity and F is the force on the particle caused by the stresses in

the fluid which can be calculated as follows:

Fimp F0 U(t + At)

Fep Fo U(t)



where ( is a high frequency resistance, and subscripts "imp" and "exp" stand for implicit

and explicit update schemes respectively.

Thus for the implicit update the equation becomes

U(t + At)(1 + i-At) U(t) + Fo At
mn m
U(t + At)(1 + -At) U(t) + U(t) + [ FAt U(t)At]
n rn rn r



We easily recognize the term in brackets as -XPAt and we can write the following set of

equations for the update scheme:

U(t + At) U(t) + AU


AU (a + )_ Fexp At,
m m
where a is our variable parameter and a = 1 for implicit and a = 0 for explicit updal

scheme and our midpoint method uses a = 0.5.

Further information on the calculation of high frequency resistance, rotational

velocity update, and discussion of stability with these schemes can be found in [2, 16].



These equations coupled with the forces generated at the solid-fluid surface due to the

fluctuating stresses, results in a discrete Langevin-like equation for the update of the

particle velocity.

2.5 Verification of the Fluctuations

The fluctuating Lattice-Boltzmann method has been previously tested with some

success, but has failed to give the correct equipartition temperature for the single

Brownian particle [4]. Nevertheless in this same work it was shown that a re-scaling of

the velocity autocorrelation function shows a good agreement with the velocity decay of a

particle set in motion with an impulsive kick, according to Onsager's "Linear Response

Theory" [15].

The current study extends the results of these papers [4, 15] and that of [2] with the

improvements in the algorithm such as better handling of the boundary nodes, addition of

a stationary mode to the 18 velocity model and an additional parameter (A,) for the decay

of the normal stress modes in an effort to study the self-diffusion in colloidal suspensions.

This section consists of the tests we have performed to verify our simulation method.

The lattice-Boltzmann method used in this study is an isothermal model and for the

fluctuating scheme and the temperature is characterized by the stress fluctuations. We

have calculated these fluctuations over a node, over a plane, and over the whole periodic

box. The average stress fluctuations are local in space and time with a variance given by

< >f 2( + 66 ) + 2 (2-3 )

The fluctuations we calculated all corresponded to the right input temperature, according

to equation 2-31 upon the use of different stress relaxation parameters (A and Av) and also

in the presence of solid particles of considerable volume fraction.

In the current model, the only interaction between a stationary planar wall and

a fluid is momentum exchange between the two arising from the fluctuations in the

fluid. A physically similar problem is an oscillating planar wall in a stationary fluid,

since the dissipation mechanism is the same in both problems. Time correlations of the

momentum exchanges provide useful information to calculate transport properties and

these correlations can be calculated either discretely or using a continuum approach.

Since the shear rate is time dependent in both cases, the frequency dependent friction

kernel is identical in both problems. Therefore we have chosen to calculate this for the

oscillating planar wall which is easier to tackle analytically. This information was then

used to calculate the momentum change correlations between the fluid and the wall for

the fluctuating fluid case (Appendix C). A discrete theory can also be constructed based

simply on the information from lattice Boltzmann theory to calculate the instantaneous

correlations (Appendix D). Figure 2-2 illustrates the agreement of the simulation results

and the theoretical expression. The theoretical expression is slightly in error at the first

few time steps and this error turns into a constant offset after the first few steps, implying

a similar first and second derivative for the two curves which determine the transport


A similar and practical problem for the simulation of colloidal suspensions is the

interaction of an infinitely massive, immobile, sphere with the fluctuating fluid. Here the

different friction kernel is different due to the differences in geometry. Using the friction

kernel [33] the momentum exchange correlations between the fluid and the sphere can be

analytically calculated (Appendix B). The agreement between the theoretical expression

and simulation values are given in Figure 2-3. Arguments similar to the case of the planar

wall for the transport coefficients also apply here. These two examples with different

geometries show that the method is capable of producing the correct transport coefficients.

The velocity auto-correlation function of a particle reveals transport properties such

as the diffusion coefficient and characteristic time scales for the relaxation of the particles.

According to the linear response theory of Onsager, the velocity relaxation of a single

particle which is given an impulsive kick is equivalent to the relaxation of the velocity

Figure 2-2.

0 200 400 600 800 1000
Momentum change correlations on a plane wall. AP(t) is the total momentum
exchange between the fluid and the particle at any instant t. Results
from simulations are compared with a theoretical expression derived from
continuum theory. The theoretical expression is higher than the simulation
data at the first time-step, error accumulates for a couple of time-steps and
then the difference between these two curves turn into a constant offset. The
first point for this correlation is calculated correctly(t1 .) from a discrete
theory starting with the fluctuations in population densities (section D).

fluctuations of a Brownian particle going through rn ii: random kicks since the means of

relaxation is the dissipation of the energy in the surrounding fluid in both cases. Figure

2-4 demonstrates this for our simulation technique and is yet another mean of verification


As illustrated in Figures 2-4 and 2-5, a single Brownian particle never reaches its

equipartition velocity fluctuations due to the artifacts of the discrete Langevin like

equation utilized in our simulations. Upon integration of this discrete equation, using

a constant friction term, to calculate the velocity fluctuations and the autocorrelation


6000 Theory /

5000 -

A 4000 -

V 3000- 1000

2000 600 -

1000 200-

0 1 0 10 20 30 40 50
0 200 400 600 800 1000
Figure 2-3. Momentum exchange correlations on an infinitely massive sphere. AP(t) is the
total momentum exchange between the fluid and the particle at any instant t.
Results from simulations are compared with a theoretical expression derived
from continuum theory.

it is found that the the velocity fluctuations differ from the equipartition value by
a factor of 2 The F sign corresponds to the explicit and implicit integration

schemes. Thus the ratio of the Brownian relaxation time step to the simulation step

p1- i, an important role for the velocity fluctuations. Figure 2-5 shows calculations of the

velocity autocorrelation function for different values of 2 and illustrates the effect of the

artifact. Yet our simplifying assumption of a constant friction term does not apply to

the simulations and the theory is only qualitatively correct. Nevertheless, the diffusion

coefficient calculated by integrating the velocity autocorrelation function according to the

Green-Kubo relationship is free of this artifact and gives the correct Stokes-Einstein value

(equation 2 32) as verified by simulations. These calculations are summarized in Appendix

Figure 2-4. Verification of Onsager's linear response theory. The decay of the velocity
correlations is compared with the velocity decay of a particle with an initial
given velocity. The two curves lead to the same diffusion coefficient upon

A. These preliminary tests show that the fluctuating lattice-Boltzmann scheme is capable

of producing the correct transport coefficients for spherical Brownian particles of

D = BT (2-32)

This concludes our test for colloidal finite sized particles. Further tests on the

hydrodynamic interactions and scaling behavior of point particles will be discussed in

C'!h pter 3.

mass factor =10

S- implicit -

t/(a vV)

Figure 2-5.

Decay of the velocity autocorrelation function for a single particle with mass
factor of 1 & 10 and a radius of 5.4, normalized by the equipartition value
U2 = kT/mB. The two different curves correspond to different integration
schemes. These two curves eventually lead to the same diffusion coefficient.

mass factor = 1




t (a v)


The numerical method, presented in Ch ipter 2, to simulate the dynamics of polymer

solutions in confined geometries has been implemented and tested. The method combines

a fluctuating lattice-Boltzmann model of the solvent [34] with a point-particle model of the

polymer chains. A friction term couples the monomers to the fluid [11], providing both the

hydrodynamic interactions between the monomers and the correlated random forces. The

coupled equations for particles and fluid are solved on an inertial time scale, which proves

to be surprisingly simple and efficient, avoiding the costly linear algebra associated with

Brownian dynamics. Complex confined geometries can be represented by a straightforward

mapping of the boundary surfaces onto a regular three-dimensional grid.

The hydrodynamic interactions between monomers are shown to compare well with

solutions of the Stokes equations down to distances of the order of the grid spacing.

Numerical results are presented for the radius of gyration, end-to-end distance, and

diffusion coefficient of an isolated polymer chain, ranging from 16 to 1024 monomers in

length. The simulations are in excellent agreement with renormalization group calculations

for an excluded volume chain. We show that hydrodynamic interactions in large polymers

can be systematically coarse-grained to substantially reduce the computational cost

of the simulation. Finally, we examine the effects of confinement and flow on the

polymer distribution and diffusion constant in a narrow channel. Our results support

the qualitative conclusions of recent Brownian dynamics simulations of confined polymers

[35, 36].

3.1 Introduction

Numerical simulations of the dynamics of polymer solutions are computationally

demanding because of the wide range of length scales and time scales in the system.

Hydrodynamic interactions are necessary to obtain even qualitatively correct scaling

laws, but Brownian dynamics, the standard simulation technique for polymer solutions,

then becomes computationally expensive. The time to calculate a single step for a chain

of N segments scales as N3 when hydrodynamic interactions are included, with the

computational cost dominated by the factorization of the mobility matrix. To calculate

average properties, the simulations must be performed over many Zimm relaxation

times, tz ~ qb3N3V/T; here r is the solvent viscosity, b is the segment length, T is the

temperature, and the excluded-volume exponent v 0.6. The number of time steps

required for comparable statistical accuracy in the diffusivity or other average property

is therefore proportional to N3". Thus, the overall computation time for a dynamical

simulation scales as N48, although an approximate factorization of the mobility matrix

[36, 37] reduces this to N3.8, but with possible loss of positive definiteness [36]. Still,

the maximum chain length that can be studied by Brownian dynamics is of the order

of 100 monomers [35, 36, 38-40], which is not ..v ,i,- sufficient to determine accurate

theological properties of flexible polymers [41]. In a semi-dilute solution of N, chains, the

computational cost of Brownian dynamics increases further, in proportion to N/ or N3

depending on the method of factorization.

In addition to the unfortunate scaling of the computational cost, Brownian dynamics

is not easily extended to include external boundaries and flow fields. No-slip boundaries

require a complicated Green function even for channel flows [36], while more involved

geometries can only be handled by finite-element or boundary-element methods. Brownian

dynamics has no straightforward way of including any external flow field beyond a simple

linear shear flow. More complex flows require a coupled calculation of the evolution of

the fluid velocity field with the dynamics of the particle chains. Even the great advantage

of Brownian dynamics, which is the capability to solve directly for the dynamics on the

diffusive time scale, is negated by the very small time step that is necessary to integrate

the stochastic equations, typically lO-stz for moderate length chains (N ~ 100). We will

see that comparable time steps can be used in inertial simulations.

In this C'!i ipter we report extensive tests of the lattice-Boltzmann method, presented

in C'!i ipter 2, for individual monomers as well as applications to the diffusion of polymer

chains in periodic and confined geometries. We show that the Oseen limit of the

dissipative and fluctuating hydrodynamic interactions between two monomers is obtained

for separations larger than Ax, the grid spacing of the lattice-Boltzmann model. We also

show that the radius of gyration, end-to-end distance, and diffusion coefficient of a flexible

excluded volume chain agree quantitatively with renormalization group calculations [42]

for chains up to 1024 monomers. A coarse-graining of the hydrodynamic interactions

has been investigated, which -, r-.- -I the possibility of a considerable speed up in the

computation for large chains. Despite the extra inertial time scale, our simulations use

time steps that are comparable to Brownian dynamics, and therefore achieve a more

favorable scaling of the computational time with chain length. We briefly investigate the

effects of confinement and flow on the polymer distribution in a two-dimensional channel.

3.2 Results

In this section we summarize the results of our simulations of polymer chains in

periodic and confined geometries. We begin with detailed tests of the numerical method,

validating the code for both dissipative and fluctuating hydrodynamic interactions.

Next we calculate the diffusion coefficient for single chains in a periodic cell, using well

established finite-size corrections to estimate the diffusion coefficient at infinite dilution

for chains up to 1024 monomers. Finally, we calculate the effects of confinement on the

diffusion coefficient, and investigate how the polymer distribution is modified by an

imposed flow field.

3.2.1 Fluctuation-Dissipation

The fluctuation-dissipation theorem implies that the diffusion coefficient of an isolated

monomer D, = T/o, where the temperature, T, is related to the variance of the random

stresses in the fluid (Eq. 2-16) and o is the monomer friction. However, as Ahlrichs and

Diinweg[ll] originally indicated, the diffusion coefficient of an isolated point particle is not







Figure 3-1.

O Dissipative
* Fluctuating

0 0.2 0.4 0.6

Relation between the effective friction coefficient and the input friction o
(Eq. 2-21). The effective friction coefficient has been obtained from the drag
force on a single monomer (circles) and the diffusion of the monomer in a
thermally fluctuating fluid (squares). The slope of the line is unity, indicating
a constant offset between -1 and $o1. The unit cell size L 20Ax.

precisely proportional to the inverse of the friction coefficient, but is offset by a constant

value that is independent of the input friction 0o,

Dm 1
T -o A'Ax

(3 1)

We therefore attempted to verify the Stokes-Einstein relation for a single monomer

by calculating the velocity of a monomer under an external body force as well as the

diffusivity of the monomer due to thermal fluctuations in the fluid. Fig. 3-1 summarizes

the most important finding, namely that the friction-coefficient determined from the

mean-square displacement of the monomer at finite temperatures, -1 = Dm/T, matches

the friction coefficient determined from the drag force, Fd = -U. Figure 3-1 also shows

that there is a constant offset between the effective friction i, measured from the drag

force or diffusion coefficient, and the input friction coefficient o (Eq. 2-21). In essence,

there is a renormalization of the friction coefficient by the flow field induced by the moving

particle. A point force F applied at the origin of a continuum fluid creates a velocity field,

1 + rr/r2
u(r) 1- + r/r- F, (3-2)

that is singular at the origin. However, the discrete nature of the lattice-Boltzmann

velocity field leads instead to a finite result,

u(0) Fg (3-3)

where g is a numerical constant. A steady-state force balance between the imposed force F

and the particle velocity,

F -Fd o [U u(0)] =- U, (3-4)

leads directly to Eq. 3-1. Further tests verify that the parameter g is independent of

solvent viscosity, confirming the relation given in Eq. 3-1. However, g does depend weakly

on system size as shown in Fig. 3-2, but not on the input friction, o. The offset, g, is

easily calculated for any cell geometry, using a measurement of the fluid velocity at the

origin of a point force together with Eq. 3-3.

When the effective size of the monomer is small compared to the grid spacing (weak

friction), qAx/l > 1 and the effect of the lattice renormalization is small. In this case

the diffusivity is controlled by the input particle radius, Dm/T -+ 1/C.,/.,,, The opposite

limit, where the particle is large compared to the grid spacing, violates the underlying

assumption of the point-particle model. However, even here the diffE' -iv i., is finite but

controlled by the grid spacing, Dm/T g/qlAx. In these numerical simulations the

lattice contribution to the diffusivity is of the order of 21 '. The fluctuation-dissipation

theorem for a discrete model does not necessarily imply equipartition of energy, as is the

case for a continuous system described by the classical Langevin equation. The data in



Figure 3-2.


0 ao=0.32
a =0.15

0.021 1 1 1 7
0 20 40 60 80
Offset between -1 and iO1 as a function of system volume. The variation in
the dimensionless mobility g (Eq. 3-1) is shown for a cubic cell of length, L.
Results are shown for input hydrodynamic radii ao 0 0.32 (circles) and ao 0
0.15 (squares). The monomer friction 0o 6= 6/ao.

Fig. 3-3 shows that the mean kinetic energy of a single monomer in a thermally fluctuating
fluid varies between 1 and 10W'. from equipartition. The results can be understood by
analyzing the Langevin equation for a single particle with constant friction,

mU= -U+F ,


where the random force, F8, has a variance (F(t)F(t')) = 2T(6(t t'). Solving this
equation leads to the well known results m (U2) = T and D = T/l.
A numerical integration of this equation can be accomplished by discretization in

U(t + At) U(t) QU(t) + AU(t),



- Implicit
-- Explicit

4- r


0 0.05 0.1 0.15 0.2

Figure 3-3. Mean kinetic energy of a single monomer as a function of the dimensionless
frictional coupling, CAt/m. The calculated kinetic energy for implicit (circles)
and explicit (squares) updates is compared with the model described in
Appendix A (lines).


U(t + At) U(t) QU(t + At) + AU(t),


where Q = (At/m and ((AU)2) 2T//m. Equation 3-6 corresponds to a first-order

explicit integration of Eq. 3-5, while Eq. 3-7 corresponds to an implicit integration.

In Appendix A, we show that this choice of random force, derived directly from the

continuous case, still satisfies the Stokes-Einstein relation D = T/l, regardless of the

integration method, but classical equipartition of energy is modified,

mKU2 li /2'

(3 -8)

explicit integration leads to the minus sign in the denominator of Eq. 3-8, while implicit

integration leads to the plus sign. Figure 3-3 shows that this simple model accounts



V 1



for the kinetic energy measured in the simulations, although here the fluctuations in

particle velocity are largely driven by the fluid. The long-time dynamics of the polymer

are controlled by the diffusion of the monomers, and the absence of equipartition at short

times will not affect these results. We note that equipartition can be recovered by using a

weaker frictional coupling or a larger particle mass, but this is not necessary for accurate

results on the diffusive time scale. Therefore the Brownian relaxation time, Q 1, does not

have to be large.

3.2.2 Hydrodynamic Interactions

The computational method described in Chapter 2 provides an accurate solution for

the flow-field around a point force, as can be seen in Fig. 3-4 where a section of the fluid

velocity field is shown for a periodic cell of length L = 1OOAx. This example is typical of

the method, which quantitatively describes the flow field on distances larger than a single

grid spacing, Ax. There is a significant difference at large distances between the Oseen

field for an infinite system (dashed line) and Stokes flow in a periodic cell (solid line),

which was calculated from the Fourier representation, [43]

1 eikr(1 kk/k2)
u(r) V F. (3-9)

The hydrodynamic interactions between a pair of monomers have been calculated

from both the dissipative response to an external force and the correlations between

particle velocities in a fluctuating fluid. The pair mobility pi, is defined in terms of the

velocity induced on one particle by a force on the other,

U, pi Fl (3-10)

and also from the autocorrelation function of the particle velocity in a fluctuating fluid,

Ai (U<(t)U (0)) dt. (3-11)
T Jo

0.01- U 1u -^
01 U -. Oseen
0.03- 0 LBE

I1 0.02-

% X 0.01-

S\ 0-
1 2 3 4 5

0.00 10 20,
1 0 2 0 3 0t "4 0 "

Figure 3-4. Velocity due to a point force. The velocity component parallel to the direction
of the force, vx(r, 0, 0), is plotted as a function of distance from the force.
Results are shown for the Oseen flow field (dashed), the Stokes flow field in a
periodic cell of length 1OOAx (solid line) and lattice-Boltzmann simulations
(circles) in the same size cell.

The linearity of Stokes flow allows us to simplify the analysis by keeping the particles

fixed in place during the course of the simulation, even though the particle velocities are


The single-particle mobility was found to be unaffected by the presence of a

neighboring particle, as would be expected for point forces. The pair mobility, which

can be decomposed into components parallel (pl) and perpendicular (p') to the

separation vector, is shown in Figs. 3-5 and 3-6. The results for dissipative and fluctuating

calculations of the mobility agree with one another and with an exact Stokes flow

calculation in the same periodic geometry, Eq. 3-9, down to separations of the order

of Ax. The discrepancies largely reflect errors in interpolating the flow field to off-lattice

A 0 Dissipative

D Fluctuating
0.5- A Stokes


0.3- O

0.2 6


0- I I
0 2 4 6 8
Figure 3-5. Comparison of the parallel component of the two particle mobility tensor,
/12. Results obtained by dissipative (Eq. 3-10, circles) and fluctuating
(Eq. 3-11, squares) simulations are compared with an independent solution
of the periodic Stokes equations (Eq. 3-9, triangles).

locations; the results shown in Fig. 3-4 were evaluated at the lattice nodes and are much
closer to the Stokes flow result at small distances. Improved interpolation should lead to

more accurate values of the mobility of closely spaced monomers.

3.2.3 Diffusion of an Isolated Polymer

We have simulated the Brownian motion of individual polymers over a wide range of

chain lengths, from N = 16 monomers to N = 1024. As demonstrated in the previous

sections, the accuracy of the hydrodynamic interactions depends on the separation
between .,-1i ,ent monomers, b, in comparison to the grid spacing of the lattice-Boltzmann

model, Ax. Consequently, three different values of b were used to calculate the equilibrium
and non-equilibrium properties of the model polymer chain; b = 2Ax, b = Ax as used by

Ahlrichs & Diinweg, [11] and b = 0.5Ax. The monomer friction and excluded volume were

0.6- 0 Dissipative
0 Fluctuating
0.5- A Stokes





0 2 4 6 8
Figure 3-6. Comparison of the perpendicular component of the two particle mobility
tensor, p12. Results obtained by dissipative (Eq. 3-10, circles) and fluctuating
(Eq. 3-11, squares) simulations are compared with an independent solution of
the periodic Stokes equations (Eq. 3-9, triangles).

also scaled in proportion to b. The computational time scales as b6, taking into account
both the change in number of grid points (oc b3) and the change in Zimm time (o b3).

The simulations were run for at least 20 Zimm times, which was sufficient to reduce the

statistical errors in the diffusion coefficient to a few percent. For the systems studied here,
the Zimm time ranges from 103At to 106At. Given the scaling of the computation time, it

is clearly advantageous to minimize the value of b/Ax.

Results for the radius of gyration Rg and the mean end-to-end distance Re are shown
in Fig. 3-7. The statistical properties of the polymer chains are independent of b and T,

and scale in the expected way for an excluded volume chain. The numerical values of the
radius of gyration fit the scaling relation Rg = 0.4bN', with v = 0.59. This is close in

absolute size to the renormalization group result for a self-avoiding random walk [42].


01 V

.-- R b~N.59-
RbN059 T=0.1, b=10

S10 T=0.1, b.- =0.5-

V T=I1.0, b=0.5Ax

16 64 256 1024
Figure 3-7. End-to-end distance and radius of gyration in a periodic unit cell. The scaling
of Re (upper points) and Rg (lower points) is shown for different monomer
separations, b, and temperatures, T. The dashed lines are the theoretical
scalings for an excluded volume chain, with exponent v = 0.59.

The lattice-Boltzmann method simulates a finite volume rather than an unbounded

domain, which has little effect on the distribution of polymer conformations but

significantly reduces the diffusivity because of hydrodynamic interactions between periodic

images. The box size was chosen to be a fixed multiple of the expected radius of gyration

to obtain a consistent set of diffusion coefficients, as shown in Fig. 3-8. The raw diffusion

coefficients, shown for two different box lengths, L ~ 3.5Rg and L ~ 7Rg, scale roughly as

N-A, but the diffusion coefficients increase with the ratio L/Rg. The finite-size effects can

be corrected by assuming that the polymer behaves hydrodynamically as a rigid sphere of

radius, Rhi, where

D = T/6lyRh. (3-12)



0 DL, L7R
o DL, L3.5R
A D,L-7R
* D, L-3.5R

s -0.59

0 5 q

0 o


o "

Figure 3-8.

Center of mass diffusion coefficient as a function of chain length in different
size cells. The diffusion coefficients are normalized by the monomer diffusion
coefficient, Dm = T/l, and are shown for periodic systems with lengths L ~
3.5Rg and L ~ 7Rg. In all cases the mean separation between the monomers
b = Ax and the temperature T = 0.1. Results are shown for the raw diffusion
coefficient (open symbols), as well as the corrected values (filled symbols).

Then we can use the mobility of a periodic array of spheres [43] to correct for the effects of
the image chains, [44]

D = DL [1 + 2.84(Rh/L) + O(Rh/L)3] .

(3 13)

Using the relation between D and Rh (Eq. 3-12), a self-consistent estimate of the diffusion
coefficient in an infinite system can be obtained. Figure 3-8 shows that the corrected
results obtained with L ~ 3 4Rg are similar to those with L ~ 6 7R,, and that
the corrected results fall precisely on the expected scaling of diffusivity with chain
length, D oc N-". The results agree well with the renormalization group prediction
for a self-avoiding random walk [42], D = 0.2T/IbN". The monomer friction in these




Sv 0 -059


0 T=0.1,b=0.5 Ax
V T=1.0, b=0.5 Ax
0 T=0.1, b=1.0 Ax
A T=1.0, b=1.0 Ax
$ T=1.0, b=2.0 Ax

10-2 I -----
16 64 256 1024
Figure 3-9. Center of mass diffusion coefficient, with finite-size corrections, for different
values of b and T. The length of the periodic cell L ~ 7Rg in each case. The
statistical uncertainties in the diffusion coefficient are smaller than the size of
the symbols, except for the longest chain, N = 1024, where they are shown

simulations was based on a hydrodynamic radius a b= /4, so the renormalization group

prediction for the ratio of polymer to monomer diffusion coefficients is D/Dm = O.N-'.
The largest simulations in Fig. 3-8 extend up to 256 monomer units and require

about 5 hours of computation per Zimm time. However the computational cost can be

drastically reduced by either reducing the separation between neighboring monomers along
the chain, b, or less dramatically, by increasing the temperature. Reducing b lowers the

accuracy of the hydrodynamic interactions between nearby monomers, while increasing

the temperature alters the Schmidt number and may prevent the Stokes flow field from
fully developing before the monomer positions change. Diffusion coefficients for different

values of b and T are shown in Fig. 3-9. In most simulations, we used b = Ax as in
Ahlrichs and Diinweg [11], but for larger chains we found that smaller values of b could be

used without affecting the accuracy of the calculated diffusion coefficient. For example,

the diffusion coefficients of chains longer than 128 monomers obtained with b = 0.5Ax

are indistinguishable from those with b = Ax. On the other hand, the smallest chains

(N < 64) require a larger value of b for a fully converged diffusion coefficient. Taken

together, these results demonstrate that an accurate calculation of the diffusion coefficient

is possible whenever the size of the polymer (R,) exceeds a few grid spacings; a reasonable

estimate is R9 > 5Ax. If this turns out to be true in general, hydrodynamically interacting

chains can be simulated with more or less constant computational cost, regardless of chain

length. Further research is needed on this point.

The scaled diffusion coefficients, D/Dm, are insensitive to temperature, except for

the smallest chains (N < 64) and the smallest mean monomer separations (b = 0.5Ax).

In these cases (N < 64, b = 0.5Ax, T = 1), the polymer diffusion coefficient exceeds

0.005AX2/At and is not small in comparison with the kinematic viscosity of the fluid

T]/p = 0.167Ax2/At. Experimentally Sc typically lies in the range 103 105, but the

simulation results show that a Schmidt number Sc > 30 is sufficient for an accurate

measure of the diffusion coefficient. This restriction is almost invariably satisfied in our

simulations, because of the large difference in diffusive time scales of the monomer and


3.2.4 Polymer Diffusion in a Channel

Unlike Brownian dynamics, where external boundaries present a serious complication,

the lattice-Boltzmann model can simulate the dynamics of confined polymers at no

additional cost. Here we consider a channel bounded in the y-direction by planar walls

separated by a distance H, with periodic boundary conditions in the x and z-directions;

much more complicated shapes can be implemented without difficulty. The monomer

units have an additional excluded volume interaction with the wall, but the hydrodynamic

calculation is unchanged. The method automatically takes account of the increased

friction of a monomer near a wall [7] and the additional hydrodynamic interactions




D, L 1H
AD L- 1H
* D, L 2H
* D, L=4H

Figure 3-10.

Center of mass diffusion coefficients parallel (Dll) and perpendicular (D') to
the channel walls as a function of channel width H; the radius of gyration of
the polymer (N = 128, b = 1) is approximately 7Ax. The confined diffusion
coefficients, normalized by the monomer diffusivity, are shown for a square
cell, L, with various ratios of L to H.

between pairs of monomers in the vicinity of the wall [45]. The simulations in the
confined geometry were run for comparable times to the bulk simulations, which typically
corresponds to 40-120 channel diffusion times, tH H2/4D. However, for the weakest
confinement, H/R, > 10, the run time was of the order of O1tH.
The diffusion coefficient for a confined polymer (N 1= 28, b = 1 and T 1= ) is
shown in Fig. 3-10 for a range of channel widths R, < H < 5R,, where R, is again the
radius of gyration in free solution. The diffusion coefficient in the perpendicular direction
is measured at the peak in the rate of growth of the mean-square displacement, after
sufficient time has elapsed for the development of the hydrodynamic interactions, but
before a significant displacement of the monomers has occurred. The no-slip boundaries
screen the hydrodynamic interactions in the direction parallel to the confining walls over







O g

distances greater than H. Therefore, the dependence of the diffusion coefficient on the

lateral dimensions L of the cell is smaller than in the fully periodic case. In general the

data shows that the diffusion coefficient perpendicular to the wall, D', is independent of

the area of the unit cell whenever L > 2H. The in-plane component of the diffusivity D1

is more sensitive to aspect ratio than the perpendicular diffusivity, requiring L > 4H for

convergent results. Moreover, when H ~ R, there is an additional constraint, L > 4Rg,

in order to prevent direct interactions between chains that are flattened by the narrow


The scaling of the parallel diffusion coefficient in a confined geometry is shown in

Fig. 3-11 for a range of confinement ratios, Rg/H, with a constant aspect ratio L/H = 4.

The diffusion coefficient DI /D of a weakly confined polymer, Rg/H
the reduction in mobility of a confined sphere [7] of radius Rh. In tighter confinement,

Rg/H > 1, there is a transition to a power-law decay, which is consistent with the

predictions of scaling theory [46], DI /D oc (Rg/H)-2/3. On the other hand Brownian

dynamics simulations in a square channel [35] -,-.- -1 an exponent closer to -1/2. This

result is supported by numerical mean-field calculations [47], which indicate that the

.,-i,! Ii I ic scaling is not reached in two-dimensional confinement until extremely long

chain lengths, in excess of 106 segments. However the confinement in a slit is less severe

than in a tube, and a mean-field calculation would be useful to confirm that scaling theory

applies to much shorter chains in one-dimensional confinement.

3.3 Conclusions

We have implemented and tested an inertial simulation of the dynamics of a single

polymer chain in solution. The method is based on a lattice-Boltzmann model of a

thermally fluctuating fluid and a point-particle representation of the interaction between

the polymer chain and the fluid. We tested the method using known results for the

dissipative and fluctuating dynamics of one and two-particle systems, and found that

the long-range Oseen level hydrodynamic interactions were reproduced quantitatively.


Faxen V \
A N=3 A H)-2/3
q 0 1N7 \
ON N15
VN= 31
A N= 63
N 127
N=255 2 \
V N= 511

0.01 0.1 1 10
R /H
Figure 3-11. Center of mass diffusion coefficient for a confined polymer relative to the
unconfined diffusivity, DH /D. The solid line indicates the Faxen correction [7]
for weak confinement, assuming the polymer lies in the center of the channel.
The dashed line is the theoretical (Rg/H)-2/3 scaling [46]. The aspect ratio of
the cell, L/H, was set to 4 in each case.

Surprisingly, the time step is comparable to typical Brownian dynamics simulations. This
fact, together with the more favorable scaling of the computational time with respect to
the size of the chain, has enabled us to perform dynamical simulations of longer chains
than has previously been possible. We obtained the expected scaling for the diffusion
constant in an unconfined system for chain lengths up to N = 1024 monomers. The
results -i .-.- -1 that a coarse graining of the hydrodynamic interactions is possible for
very long chains, keeping the radius of gyration at a fixed number of grid points of
the underlying fluid. In this case, the computational cost of calculating hydrodynamic
interactions is independent of chain length. The method can be readily extended to
complex confinement geometries, with or without an imposed flow field. Our results for
confined channels are qualitatively consistent with Brownian dynamics simulations in a

square duct [35, 36]. Although all these tests have been for a single chain, multi-chain

simulations are straightforward and require minimal additional computation if the system

volume remains fixed.


We investigate the lateral migration of a confined polymer under pressure driven

and uniform shear flows. We employ a hybrid algorithm which couples point particles to

a fluctuating lattice-Boltzmann fluid. We observe migration in both uniform shear and

pressure driven flows, supporting the idea that migration is driven by a combination of

shear and hydrodynamic interactions with the wall, rather than by the shear gradient.

Recent numerical and theoretical investigations have -, .-.- -1. 1 that polymers migrate

towards the centerline when hydrodynamic interactions are included, but our simulations

show that in sufficiently narrow channels there is a reversal of direction and the polymers

move towards the wall. We also report on differential velocities and reduced dispersion

coefficients in pressured driven flow.

4.1 Introduction

The transport of polymer solutions through microfluidic devices depends on polymer

conformation and the center-of-mass distribution across the channel. In the absence

of flow, steric repulsion between a globular polymer and the confining walls creates

depletion 1l.,-_iS on the order of the radius of gyration, R, [48-50]. In a pressure-driven

flow, the shear rate stretches and aligns the polymer along the flow direction, reducing its

configurational entropy. Thermodynamic arguments therefore predict that the polymer

will migrate to the centerline where the local shear rate is minimized [51, 52], creating

an apparent slip at the boundaries [53]. The thickness of the depletion 1 v, is observed

to increase in a pressure-driven flow [54, 55], which also implies a net migration of the

polymer towards the channel center. By contrast, it has usually been assumed that the

depletion 1-v-r thins in a uniform shear flow due to reduced steric hindrance between a

stretched polymer and the boundary [53].

Recent theory and simulations have underlined the importance of hydrodynamic

interactions in describing polymer migration. Kinetic theories [56, 57] attribute lateral

migration to the .,-viiI1.1 I1 y of the hydrodynamic interactions between a stretched

polymer and a planar no-slip boundary, which produces a net lift away from the wall.

Computer simulations [17, 36], including full hydrodynamic interactions, showed that the

polymer migrates towards the channel centerline in pressure driven flows. On the other

hand, simulations neglecting hydrodynamic interactions between polymer segments [39, 40]

find migration towards the walls.

There are qualitative differences between the predictions of kinetic and thermodynamic

theories of migration. In particular, in a uniform shear flow, no migration is expected on

thermodynamic grounds since the polymer extension is independent of position. However,

the kinetic theories predict migration towards the channel center in a uniform shear flow

as well as Poiseuille flow. To distinguish between these theories, we have used numerical

simulations, with hydrodynamic interactions, to calculate the direction and extent of the

lateral migration in pressure-driven and uniform shear flows. We discovered that migration

can occur in both directions, depending upon the degree of confinement; in wide channels

the polymers move towards the center, while in narrow channels they move towards

the walls. The polymer distribution is quantitatively different in shear and Poiseuille

flows, but the direction of migration is the same in both cases for the same degree of

confinement. This -~i-:.- -I that the primary driving force for migration is hydrodynamic

lift from the confining walls, driven by the shear rate rather than its gradient [56, 57].

4.2 Method

The polymer in this work is modeled by beads connected by stiff linear (Fraenkel)

springs of finite rest length, each spring representing a Kuhn length. The dynamics of the

fluid phase is obtained from the fluctuating lattice-Boltzmann model [4, 21], presented in

Ch'! pter 2. The polymer motion is coupled to the surrounding fluid using a frictional force,

based on the difference between the bead velocity and the local fluid velocity [11, 29].

The force exerted by the beads is redistributed to the surrounding fluid nodes so that the

total momentum of the coupled particle-fluid system is conserved. Extensive numerical

tests have shown that we obtain Oseen level hydrodynamic interactions between beads

and between beads and the solid boundaries [17]. In addition, there is a soft repulsion

between pairs of beads and between beads and confining walls. The excluded volume

sphere has a radius of about b/2, where b is the distance between beads. The fluctuating

lattice-Boltzmann approach has computational advantages over Brownian dynamics [58],

which has been the primary numerical method for simulating polymer solutions. First,

the computational cost scales linearly with the number of beads, allowing for simulations

of chains in excess of a thousand beads [17]. Second complex geometries, such as porous

media or micro devices, can be introduced by mapping the solid surfaces onto the lattice,

without introducing a new Green function for each specific geometry.

We study a single polymer chain, confined between two flat plates separated by a gap

H. Periodic boundaries are used in the other two directions and the periodic length in

these directions is sufficiently large (L = 2H) that the hydrodynamic interaction between

distant images is negligible [17]. The polymer chains are discretized into N beads with

most simulations reported in this paper using N = 16. Longer chains (N = 32 64) are

used to verify that the equilibrium radius of gyration (R.) and the free solution center

of mass diffusion coefficient (D) are the only important parameters and that the level of

discretization of the chain is sufficient. A pressure driven flow was created by applying a

uniform force to the fluid, while a uniform shear flow was generated by relative motion

of the bounding walls. The Peclet number is based on the shear rate and the equilibrium

radius of gyration, Pe -= R2/D, in both cases. In a uniform shear flow, the shear rate is

constant across the channel and given by = AV/H, where AV is the velocity difference

between the walls. In Poiseuille flow, the mean shear rate, 7 = 2Vmax/H, is used, where

Vmax is the centerline velocity. The degree of confinement is characterized by the ratio of

the gap width to the equilibrium radius of gyration of the polymer chain, H/R9. For the

smallest confinement ratio used in this work, H/Rg = 3, the polymer coil remains free to

rotate in the channel. The simulations are run for a sufficiently long time that a polymer,

under equilibrium conditions, would diffuse a distance of at least 50H. The statistical

errors in the figures are less than 5'. in the bulk region.

4.3 Results and Discussion

Simulation data for wide channels (H/R, > 5), shown in Fig. 4-la, confirm that

polymer chains migrate towards the channel centerline in a pressure driven flow. The

mechanism for this migration is not entirely agreed upon. The thermodynamic theory

argues that the shear gradient drives the polymer towards low shear regions to attain

maximum possible configurational entropy [51, 52]. We have tested this assumption

by repeating the simulations in a uniform shear flow. The simulations show that the

polymers still migrate towards the centerline (Fig. 4-2a), indicating that a shear gradient

is not required for migration. However, due to symmetry, no migration would occur in

an unbounded shear flow; therefore shear cannot by itself explain the migration either.

Hydrodynamic arguments -~'--.: -1 that the tension in the sheared polymer generates a

flow away from the wall, due to the .,-v:ii:ii. I ry of the hydrodynamic interactions between

a pair of beads and a planar boundary [55, 57]. A kinetic theory [57] of a dumbbell

near a planar wall predicts migration towards the channel center in both uniform shear

and Poiseuille flows when hydrodynamic interactions are included, in agreement with

simulations. A dumbbell is much less flexible than the multi-bead chains used in our

simulations. Moreover, the hydrodynamic center of a dumbbell is also its center of

mass, so an external torque leads to a pure rotation rather than tumbling. Despite

these limitations, the kinetic theory qualitatively explains the key features observed

in simulations of flow in wide channels. A feature of the polymer distribution in wide

channels is a symmetric double peak (Fig. 4-la) observed at high Peclet numbers [17, 36].

This is a consequence of the interplay of hydrodynamic forces with a spatially varying

diffusivity [57]. A sheared polymer undergoes significant extension and alignment along

the flow axis (Fig. 4-3a), hence the diffusivity along the confinement direction is reduced

with increasing shear rate. The competition between the spatially varying diffusivity

Pe=0 (a) 0.5 (b)
0.25 ....... Pe 50 / 0.4 -
.-.-. Pe-100 /..- --
0.2 -- -P e 200 g .......... ." 0.3-
C015 ./ ...20. 0
0.15 .2
0.1 -
-0.1L / ~ //

0.05- HR 8 0 HR 5
0 0I I
0 1 2 3 4 0 0.5 1 1.5 2 2.5

(c) (d)
0.5 0.8-

0.3 0.6-
0.3/ 0.4-
0.2- /.'4
0.1- H/R -4 0.2- H/R 3
0 0.5 1 1.5 2 0 0.5 1 1.5
y/R y/R
Y g Y g

Figure 4-1. Center of mass distribution for pressure driven flow as a function of Peclet
number and channel width. Results are for N = 16. The distributions are
normalized such that fJ R P(y/R,)d(y/R,) 1. Due to the symmetry of the
system, the distribution is plotted up to 0.5H. The bounding wall is at the

and the hydrodynamically induced migration leads to the double peak in the polymer

distribution function P(y/R,) at high Peclet numbers. Kinetic theory [57] would also

predict a double peak in the distribution function at high Peclet numbers if a spatially

varying diffusivity was used in the model.

Fig. 4-la shows that the peak position in the channel is independent of Peclet

number. Empirically, we observe that the position coincides with the polymer reaching

about 911'. of its maximum extension in the flow, as seen by comparing Figs. 4-la and

4-3a. The mobility along the confinement direction, which is a logarithmic function of

the aspect ratio, changes little with the marginal additional extension that occurs in the

vicinity of the walls. Consequently, the driving force away from the centerline, which

is proportional to the gradient of the diffusivity [57], diminishes very quickly beyond

.... Pe=35 A 0.8-
-- Pe=70 ._.9
02- *' 0.6 -

... ,"0.4 -
0.1- /.
,- H/R =8 0.2- /. H/R = 3

0 1 2 3 4 0 0.5 1 1.5
y/R y/R

Figure 4-2. Center of mass distribution for uniform shear flow as a function of Peclet
number and channel width. Results are shown for N = 16. The distributions
are normalized such that fj "g P(y/R,)d(y/R,) 1.

this point. In a uniform shear the distribution has at most a single maximum, since the

diffusivity is essentially constant.

As the channel width is reduced with respect to the size of the polymer (H/R, = 5),

the migration towards the center decreases (Fig. 4-1b). The extension data (Fig. 4-3b)

shows that the polymer shape is similar to that in the wider channels, indicating that the

reduced migration is a result of the smaller hydrodynamic force. At a channel width of

H/Rg = 4, migration towards the center is only observed for the highest Peclet number.

At lower Peclet numbers the migration is towards the walls (Fig. 4-1c). In still tighter

confinement (H/R, = 3), migration is towards the walls at all Peclet numbers. We observe

a similar outwards migration in uniform shear flow (Fig. 4-2b), although the extent of

the migration is lower than in Poiseuille flow. Migration towards the channel walls is

further supported by results for the second moment of the center-of-mass distribution [36]

(Fig. 6). Nevertheless, the evidence for outward migration in this work is by no means

unequivocal; indeed the authors conclude that the polymer alv-,- migrates towards to the

center if hydrodynamic interactions are included.

Polymer kinetic theory [57] accounts semi-quantitatively for migration in wide

channels, correcting earlier work [56] based on similar ideas. The key migration mechanism,

namely the lift from the wall induced by a sheared polymer, has been validated in previous


2 2- **o* 2 e-*******

r 1- s o g a o @ o o o gi 1- & o^oooooggooaooooo

HRg 8 o HR -5

0 1 2 3 4 0 0.5 1 1.5 2 2.5

3 (c) 3 (d)

< 2- *o- 01O o ooooooooo2

1 1- O0000000000000000000000

uf /HR =-4 H/R -3
C 1
0 HRg 0 0004 I I

0 0.5 1 1.5 2 0 0.5 1 1.5
y/R y/R
Y g Y g

Figure 4-3. Components of the radius of gyration along the flow, (x2), and confinement,
{(y2, axes. Results for chains of length N 16 are shown at two different
Peclet numbers, normalized by Rg ((x2) + (y2) + (z2))o05.

work [17, 36] and confirmed in the present study. However, the hydrodynamic interactions

in the channel were approximated by superposing the flow fields generated by each wall.

This is valid when the polymer is much closer to one wall than the other [7, 59], but

in a narrow channel the Green's function for a slit must be used. In narrow channels,

superposition over-predicts the extent of the hydrodynamic interactions, and therefore

the degree of migration towards the center. In reality, the hydrodynamic forces that

cause migration towards the center are screened by the closely spaced boundaries. We

have found that the polymer diffusion coefficient along the confinement axis shows

a cross-over to Rouse-like scaling in the range 8 > H/Rg > 4, which supports the

screening hypothesis. Moreover, computer simulations [39, 40] which neglect intra-bead

hydrodynamic interactions also show migration towards the walls, independent of the

channel size.

oPe 0,x *Pe 50,x (a)
nPe O,y oPe 50,y

0.2 CA

S0.5 1 1.5
Y g

Figure 4-4. Center of mass distribution for different levels of discretization of the polymer
in a tightly confined channel, H/R, = 3. The distributions are normalized
such that fJ1 P(y/R )d(y/R ) 1.

The shape of the polymer may also p1 i, an important role in the reversal of the

direction of migration. The extension curves (Figs. 4-3c and d) -,i--.- -1 that equilibrium

conformations become distorted within 1.5R, of the walls. Therefore highly confined

polymers are anisotropic even at equilibrium and there is no bulk region where the

polymer assumes its unconfined globular form. This indicates an analogy with anisotropic

rigid particles which also migrate towards the boundaries under flow. This migration is

due to the particle aligning in the direction of flow, thereby reducing the diffusivity in the

confinement direction [60, 61]. Alignment also reduces the thickness of the depletion l1 r

as a result of a reduced steric force.

The extension of the polymer and the distribution of polymer conformations depends

on chain length, especially at high Peclet numbers [62]. Resolution of such conformational

changes may require a finer discretization of the polymer chain. To assess the accuracy


1.3 -HR 8
A[s H/R=8
1.25 9g
A R =16

v 1.2- Fit


0 50 100 150 200 250 300
Figure 4-5. Average axial velocity of the center of mass of polymers normalized by the
average fluid flow. The results are presented with respect to the polymer
Peclet number for different channel sizes.

of the results presented, we have repeated some of the calculations using chains of 32
and 64 beads. A finer discretization allows for more conformational flexibility and at
the same time improves the accuracy of hydrodynamic interactions [17]. We scale the
system volume and pressure gradient to maintain constant confinement ratio and Peclet
number, which results in a finer discretization of the fluid grid in comparison to the radius
of gyration of the polymer. Longer chains at constant Peclet number also result in lower
Reynolds numbers, inversely proportional to their equilibrium radius of gyration. The
Reynolds numbers used in this work range from 0.6 to 2.3 for the shortest chains in the
largest channels and 0.15 to 0.5 for the longest chains. The results shown in Fig. 4-4
indicate that the level of discretization with N = 16 is sufficient for this problem. Longer
chains show identical distributions within the statistical errors.
Polymer migration also implies a differential velocity based on polymer chain length
[36, 53]. Two polymer chains of different length have different polymer Peclet numbers

under the same flow conditions. This ii--. -1i that the longer polymers will migrate

towards the channel center more than their shorter counterparts and as a result the

longer chains will have a higher average velocity than the shorter chains. Figure 4-5

shows the mean polymer velocities Up = (AX(t)) /t compared to the mean fluid velocity,

Uf = 2Umax/3. The polymer velocity, (Up) / (Uf), is a universal function of the polymer

Peclet number, Pe, at high Peclet numbers. However, at low flow rates the polymer

velocity is solely determined by the excluded volume.

The simulation trajectories were then used to calculate the first and second cumulants

of the the axial displacement of the center of mass of the chains, (AX(t)) and (AX(t)2)

(AX(t))2. The motion of the polymers were found to be diffusive on time scales larger
than H2/D and results for the axial dispersion coefficients are shown in Fig. 4-6 (solid

symbols). The Taylor-Aris calculation of the dispersion coefficient is an .1-i', i.d- I ic

solution of the convection-diffusion equation, assuming a uniform distribution of particles

across the channel,

+ Pe(4

where the coefficient 3 is 1/210 for tracer particle in a plane Poiseuille flow. It can be

seen that the numerical results are at least an order of magnitude smaller than Eq. 4- 1

even at low Peclet numbers. Note that the molecular contribution to D* is negligible

under all flow conditions studied here. If we take into a account the finite size of the

polymer, the Taylor-Aris result is modified because the polymer does not sample the high

shear regions near the wall. The center of mass is restricted to positions between Ld and

H Ld where Ld is the thickness of the depletion -1*v r. Assuming a uniform distribution

across a restricted channel width Ld < y < H Ld, it is straightforward to show that

3 = (1 2Ld/H)6/210. If we take the depletion 1 rv. r thickness, Ld, to be equal to the

radius of gyration, Rg, then we obtain the solid line shown in Fig. 4-6 in good agreement

with numerical simulations at the same confinement and at low Peclet numbers. We have





Figure 4-6.



Hydrodynamic dispersion coefficients of a polymer chain in Poiseuille flow.
Numerical simulations are shown as filled symbols, along with Taylor-Aris
theory (TA) for tracer particles (Ld = 0, dashed line) and for finite sized
particles with a fixed depletion lv.-r (Ld = Rg = H/8, solid line). Theoretical
calculations, accounting for the depletion l-iv.r thickness determined from
simulation, are shown as open symbols. The standard error in the simulation
data corresponds to the size of the symbols.

found that the leveling and saturation of the dispersion coefficient in these simulations are

due to inertial effects which become prominent only beyond Pe > 4000. On rerunning

similar simulations at lower inertia we have found out that the dispersion coefficients are

still greatly reduced but do not saturate at the high Peclet limit. We show the results

of these simulations in Fig. 4-7. We also show the effect of inertia on the center of mass

distribution in Fig. 4-8, but we do not find a great impact on the distribution apart from

a slight change in the dip in the middle of the distribution.

It is often r- .--. -1. 1 that the differential velocity of different chain lengths could be

used to separate chains of different lengths, but Taylor dispersion tends to spread the

distribution too much for it to be a useful separation technique. Here we briefly examine

op B

,'" Simulation, H/R 8

S TA, LdR HR 8
-- TA, Ld0O

10000- U


100 1000
Figure 4-7. Comparison of dispersion coefficients with different inertia. The highest
particle Reynolds number for kT = 0.1 is almost 3.4, whereas the channel
Reynolds numbers are up to 30 for the same simulations.

the consequences of the reduced hydrodynamic dispersion. If we take the saturated

dispersion coefficient from Fig. 4-6, D*/D ~ 100(H/Rg)2, then the expected width of

the axial distribution passing through a channel of length L is roughly 10/-LH/Pe1/2

A/L/H, since a typical polymer Peclet number is of order 100. On the other hand the

segregation due to the differential velocity is of the order L(AU/Uo) where AU is the

differential velocity of the chains. Thus chains with a differential velocity of the order of

Uo-H/L could be separated by this means. Unfortunately aspect ratios of L/H in excess

of 104 seem impractical, so at best chains of differential length of AN/N ~ 10-2 might be

separated by such means.

4.4 Conclusion

In this work we have shown that the extent and direction of the migration, in both

uniform shear and Poiseuille flow, depends on the degree of confinement as well as the

Peclet number. The dynamics in the highly confined channels differ significantly from

1.5 -

1 -


1 2 3 4
Y g

Figure 4-8. Comparison of center of mass distributions at 3 different temperatures at
Pe = 300. The three different temperatures correspond to different Reynolds
numbers. The highest particle Reynolds number at kT 0.1 is around 3.4.

that in wider channels, where hydrodynamic interactions are stronger. The hydrodynamic

screening from the walls reduces the migration towards the centerline at high confinements

and alignment in the flow can then produce migration towards the walls, similar to

what is observed in simulations without intra-bead hydrodynamic interactions [40]. We

have confirmed that the primary mechanism for migration towards the centerline is a

combination of shear and hydrodynamic interactions between the walls and the beads

[55-57]. We also found that the spatially varying diffusivity produces a double maximum

in Poiseuille flow [17, 36], but not in a uniform shear flow. Systematic experiments in

comparable parameter regimes would be an interesting confirmation of the results of

this study and other recent simulations [17, 36]. We have also reported on the dispersion

and differential velocities in such systems along with a discussion of the limitations of

separations based on our results.


We demonstrate that a polymer confined to a narrow channel migrates towards the

center when driven by an external force parallel to the channel walls. This migration

results from .,-viiiii, li ic hydrodynamic interactions between polymer segments and

the confining walls. A weak pressure-driven flow, applied in the same direction as the

external force, enhances the migration. However, when the pressure gradient and the

external force act in opposite directions the polymer can migrate towards the boundaries.

Nevertheless, for sufficiently strong forces the polymer alv- -, migrates towards the center.

A dumbbell kinetic theory explains these results qualitatively. A comparison of our results

with experimental measurements on DNA -, '- -- -I that hydrodynamic interactions in

polyelectrolytes are only partially screened. We -i-.-, -1 new experiments and theory to

test the extent of the screening in polyelectrolyte solutions.

5.1 Introduction

Flexible polymers in a pressure-driven flow field migrate towards the center of the

channel, because of hydrodynamic interactions. The local shear rate stretches the polymer

and the resulting tension in the chain generates an additional flow field around the

polymer. This flow field becomes .,-vmmetric near a no-slip boundary and results in a net

drift towards the center of the channel [36, 55, 57]. Recent simulations [57, 63] show that

hydrodynamic lift is the dominant migration mechanism in pressure-driven flow, rather

than spatial gradients in shear rate. Since the migration depends on chain length [36, 63],

it is in principle possible to fractionate polymers based on their length; longer polymers

migrate more strongly and therefore elute faster. However, Taylor dispersion precludes any

sharply peaked distribution, so there is considerable interest in using combinations of flow

and electric fields to improve separation efficiency [64-66]. We used numerical simulations

to investigate a flexible polymer driven by a combination of fluid flow and external force.

In this preliminary investigation we have not considered the effects of counterion screening.

We find that the polymer migrates towards the channel center under the action of a

body-force alone, while in combination with a pressure-driven flow, the polymer can move

either towards the channel wall or towards the channel center. A kinetic theory, based on

a dumbbell model of the polymer, gives a qualitative understanding of the results. The

simulations mimic recent experimental observations on the migration of DNA in combined

electric and pressure-driven flow fields [66, 67]. The similarities between these results

~-'.-.- -1 that hydrodynamic interactions in polyelectrolyte solutions are only partially

screened [68].

5.2 Methods

Numerical simulations were used to investigate the motion of a confined polymer

chain, driven by a combination of fluid flow and external forces. The system is bounded

in one direction by no-slip walls, separated by a gap H. Periodic boundaries were used

in the other two directions, with a repeat length (L = 2H) such that the hydrodynamic

interactions between periodic images are negligible [17]. The external field and pressure

gradient result in two different Peclet numbers: Pe = UR,/D and Pef -= Rg2/D. Here U

is the average polymer velocity with respect to the flow field, Rg is the equilibrium radius

of gyration, D is the free-solution center of mass coefficient, and 7 is the average shear

rate. We employ a hybrid algorithm, coupling point particles connected by stiff springs to

a Newtonian fluid [11]. The fluid is simulated by a fluctuating lattice-Boltzmann equation

[3], which accounts quantitatively for the dissipative and fluctuating hydrodynamic

interactions between polymer segments [17]. The polymer chains are discretized into

N 1 segments of length b, with most of the simulations using N = 16; additional

simulations with N = 32 were used to check the effects of chain discretization. Previous

work [63] has verified that the radius of gyration, R,, and the diffusion coefficient, D, scale

.I ',' ii ically (No'6) at this level of discretization.

5.3 Results

The numerical simulations show that a flexible polymer chain migrates towards

the channel center when subjected to a strong external force parallel to the channel


Figure 5-1. Center of mass distribution of a single polymer in a channel H/R, = 8 under
the application of an external body force. Half of the distribution is shown,
with the boundary at y/H = 0 and the center of the channel at y/H = 0.5.

boundaries. At Peclet numbers in excess of 100, this transverse migration results in a

nonuniform center of mass distribution (Fig. 5-1), with the polymer concentrated in the

middle of the channel. For small forces (Pe < 10), Brownian motion is dominant and

the distribution eventually becomes uniform, except for a small depletion 1v.-r near the

channel walls. However, if hydrodynamic interactions are neglected there is no migration,

and the distribution is uniform at all Peclet numbers. This observation emphasizes the

hydrodynamic origin of the migration.

In the limit of infinite Peclet number the distribution still has a finite width, which

shows that flexible polymers have an additional dispersive mechanism besides Brownian

motion. This hydrodynamic dispersion is due to the chaotic motion of the interacting

polymer segments and is similar to the dispersive mechanism in settling suspensions

[69, 70]. It p-vq an important role in limiting the migration. A polymer in the vicinity of

Pe =
-oPe = 1100
4 ----Pe = 110
O-OPe = 11
x-xPe = 110, No HI

a wall rotates due to the hydrodynamic interactions with the boundary, and drifts away

from the wall. However, the drifting polymer changes conformation through hydrodynamic

dispersion and eventually assumes a more spherical shape with limited further drift.

Without this dispersion the polymer would traverse back and forth across the channel,

similar to a rigid rod [71, 72], and the center of mass distribution would be uniform.

On the other hand a spherical particle would not migrate at all and the distribution

would again be uniform. Flexibility and dispersion are therefore crucial in determining

the distribution at large Peclet numbers. The hydrodynamic dispersion can be as large

as Brownian diffusion, as can be seen by comparing the width of the infinite Peclet

distribution with a finite Peclet number case, Pe 1100 (Fig. 5-1). Note that the

statistical errors for the infinite Peclet case are approximately 10'(. so these distributions

are statistically indistinguishable.

A theoretical approach can be used to gain further insight into the hydrodynamic

origin of the migration. We have used a kinetic theory similar to Ref. [57] and have

identified three mechanisms that lead to migration; two of these have been noted

previously [57, 66], one is new. We analyze the evolution of the distribution function of

a dumbbell near a planar boundary, incorporating a uniform external field in addition to

the hydrodynamic flow field [57]. The key assumption is that the distribution function can

be factored into a product of the center of mass distribution, P(y, t), and the orientation

distribution T'(qj, y, t), where the vector qj connects the monomers and y is the distance of

the center of mass from the boundary. The expression for the lateral flux (jy) is then

+- < Dy F > < Dyy >

PkT U- 9lnJ aln4
+ 2 < y > -P< D >, (5-1)

where the angle brackets denote an average over the distribution function T. In this

equation, ?9yjFF is the instantaneous lift velocity [57] of the polymer and Dij is the

Kirkwood diffusivity. The external force, FE, is directed along the central axis of the

r v


Figure 5-2.

Illustration of different migration mechanisms. In each figure the solid lines
represent the forces while the dashed lines represent the velocities generated
by these forces. The velocities can be calculated using an image system for
a point source near a planar boundary. A) The lift due to shear: the tension
in the polymer near a solid boundary generates a net velocity away from the
boundary. B) Rotation due to the external field: the velocity field due to the
external force on each bead results in a rotation about the center of mass of
the polymer. C) Drift of a rotated polymer: two beads oriented at an angle
to the external forces drift due to the hydrodynamic interactions between the

channel. We have further assumed that the dumbbell is connected by a linear spring,

F =- -kqj, and that the orientation distribution can be expanded in powers of the local

shear rate and external force. This simplifies the time-independent solution for the center

of mass distribution function, P(y),

0- [A + B P2 + CPF- kT aP
6y y2 y2 2( 6Oy


where 7 is the local shear rate, is the hydrodynamic resistance for a monomer and A,

B, and C are numerical coefficients. The calculation will be described more fully in a

subsequent paper which will include a more quantitative comparison between theory and

simulation. As in Ref. [57], this solution for a dumbbell near a single wall can be extended

to a channel flow by superposing the solutions for the two walls.

The first term in Eq. 5-2 describes the lift created for the center of mass of a

dumbbell by the .,i-v ii, I i ic hydrodynamic interactions between the polymer segments

and the walls [57]. The local shear rate 7 generates a tension in the polymer, which

couples with the linear (in 7) contribution to T. This coupling creates a nonuniform

distribution in pressure-driven flows as seen in recent experiments [55] and simulations

[17, 36, 63]. The second term describes the rotation and drift of the polymer towards the

channel center under the application of an external field. Although both terms arise from

hydrodynamic interactions with a boundary, the mechanisms are different in the motion

they create. In a shearing flow the polymer is in tension, and the opposing forces generate

a lift for the center of mass [57] (Fig. 5-3A), but in an external field the force on the

monomers has the same sign, which causes the polymer to rotate .iv from the boundary

(Fig. 5-3B). Subsequently the polymer drifts towards the center under the action of the

force. Finally, there is a cross term that arises from a coupling between the stretching

and rotation of the polymer by the flow and the external force (Fig. 5-3C). The direction

of the drift now depends on the sign of F, but is independent of the boundaries. This is

the mechanism proposed in Ref. [66], but it does not explain polymer migration in the

presence of a body force alone.

The numerical simulations -,-. -1 that the center of mass distribution across the

channel, P(y, Rg, H, D, U), can be written in a scaled form P(y/H, Pe). Figure 5-3a

shows center of mass distribution of a 16 bead polymer chain at various channel widths

and polymer Peclet numbers. Polymers with the same Peclet numbers have the same

distribution regardless of the size of the channel. Figure 5-3b shows that center-of-mass

distribution is independent of the number of segments used to discretize the chain.

However, the kinetic theory does not capture these scalings, which remains unexplained at


The external forces can be used in conjunction with a pressure driven flow such

that an axial velocity difference is created between migrating polymers based on size.

If the external force and pressure gradient drive the polymer in the same direction, the

theory predicts that the coupled term enhances the migration towards the channel center

(Eq. 5-2). The numerical results in Fig. 5-4 verify this prediction, as can be seen by

comparing the results (Pe = 110) in Fig. 5-1 with Fig. 5-4 (Pef = 12.5). Since the flow

Figure 5-3.

"0 0.1 0.2 0.3 0.4 0.5

a) Scaling of the distribution function with respect to the channel width
under constant polymer Peclet number, Pe. b) Scaling of the distribution
function with respect to polymer Peclet number using different levels of chain
discretization. The boundary is at y/H = 0 and the center of the channel is at
y/H 0.5.

is relatively weak (Pef 10), it makes only a 10 21'. contribution to the migration.

However, a mixture of polymers in such a system has significantly different elution rates;

larger polymers migrate more strongly and therefore sample higher velocity streamlines

near the center of the channel. Since the migration is primarily driven by the external

field, a weak hydrodynamic flow could be used for fractionation, thereby reducing the

Taylor dispersion.

In a countercurrent application of the two fields the polymer tends to orient itself

in different quadrants, depending on the relative magnitude of the two driving forces.

The polymer now drifts either towards the walls or towards the center depending on its

mean orientation. The results in Fig. 5-5 show migration towards the boundaries when

Figure 5-4.


Center of mass distributions for concurrent application of an external body
force and pressure-driven flow. The solid curve shows the level of migration
under the pressure-driven flow only. The flow Peclet number in all three cases
is Pef 12.5. The boundary is at y/H = 0 and the center of the channel is at
y/H =0.5.

the external force is small (Pe < 30), but increasing the force eventually reverses the

orientation of the polymer and the polymer again migrates towards the center (Pe > 100).

We note as an empirical observation that the peak of the polymer migration towards the

wall (Pe = 44) occurs when the polymer is driven with almost equal but opposite velocity

by the flow and force.

The numerical results and kinetic theory predictions are in qualitative agreement with

a series of experiments measuring the distribution of confined DNA under the combined

action of an electric field and a pressure-driven flow [66, 67]. These experiments showed

a strong migration towards the center in the concurrent application of force and flow

and also the reverse migration towards the boundaries in the countercurrent application.


SPePe = 22
0.5- / Pe = 44
0Pe = 66
=----Pe = 88
A--aPe = 110
0.1 0.2 0.3 0.4 0.5

Figure 5-5. Center of mass distributions for countercurrent application of an external body
force and pressure-driven flow. The solid curve shows the level of migration
under the pressure-driven flow only. The flow Peclet number in all cases is
Pef 12.5. The boundary is at y/H = 0 and the center of the channel is at
y/H 0.5.

However the experiments do not show a strong migration when the electric field is applied
by itself. This -i.--.- -1I that the dominant mechanism in the laboratory experiments
[66] is the coupling between the local shear-rate and the applied force (Eq. 5-2), rather
than rotation of the polymer by hydrodynamic interactions with the walls. The channel
dimensions in the experiments were about 10 times larger, in comparison with the polymer
size, than those in the simulations, which reduces the importance of wall effects. More
crucially, the hydrodynamic interactions in polyelectrolyte solutions driven by an electric
field are screened by the counterion motion [73, 74]. However, if there is no significant
hydrodynamic interaction between segments, then the polymer will follow the direction
of the external force without significant transverse motion, even when deformed by a

shear field. The rigid rod model cited in Ref. [67], only migrates in an external force

because of hydrodynamic interactions between segments. We have verified this behavior

with numerical simulations in which hydrodynamic interactions were excluded. Thus the

observation of migration in these experiments [66, 67], implies a degree of hydrodynamic

interaction on the scale of the polymer.

It has been predicted that the fluid velocity around a polyelectrolyte has a dipolar

component in an electric field, decaying as 1/r3 [68, 75]. This is in addition to the

angle-averaged velocity field, which decays exponentially [73]. For an anisotropic

polyelectrolyte, distorted by the shear flow, the dipolar contribution gives rise to a

transverse velocity of magnitude F/pi92R3 [68], where F is the total electric force on the

polymer, rI is the fluid viscosity, and K is the inverse Debye length. This velocity is small

compared with the electrophoretic velocity F/qIL; for A-DNA in 5mM salt solution, the

migration velocity is then predicted to be of the order of 1 of the electrophoretic velocity.

This is consistent with experimental measurements (E = 40 V/cm) [66], which -i-.-. -1

migration velocities of the order of 10-3 cm/sec, in comparison with electrophoretic

velocities of the order of 10-1 cm/sec. The theoretical estimates could be made more

quantitative using the kinetic theory outlined in this paper and described in more detail in

Ref. [57].

Additional experiments could be performed to probe the screening length associated

with the hydrodynamic interactions in polyelectrolyte solutions. We note that Fig. 3 in

Ref. [66] may indicate a weak force-induced migration, although the statistical errors

in the data make it inconclusive. A narrower channel would be helpful in promoting

hydrodynamic interactions between the polymer and the channel walls and thereby

strengthening any force-induced migration. The magnitude of the applied force could be

increased by using low-frequency AC fields; since the force-induced migration is quadratic

in the field, the direction of migration is independent of the sign. This should make it

possible to clearly decide if there is any force-induced migration or not. Screening could

also be investigated more directly by numerical simulations that include electrostatic

forces, either through explicit charges or through a Poisson-Boltzmann approximation

[76, 77].

5.4 Conclusion

In this work we have shown that confined polymers migrate towards the center of a

channel under the application of body forces, due to hydrodynamic interactions between

polymer segments and the channel walls. We have developed a kinetic theory that explains

these observations at least qualitatively. We note that our observations of polymer

migration mirror those of recent experiments with DNA [66, 67]. The similarities between

the migration of a neutral polymer driven by a body force and a polyelectrolyte driven by

an electric field SI,'-': -1 that hydrodynamic interactions in polyelectrolyte solutions are

only partially screened [68]. Indeed, if they were fully screened [73] the observed migration

would not occur. We have r-t'---: -1. 1 additional laboratory and numerical experiments to

further investigate this question.


We have implemented and tested a fluctuating lattice Boltzmann model for the

simulation of colloidal and polymeric suspensions. With this model, we have specifically

investigated the equilibrium and dynamic behavior of dilute polymer suspensions within

periodic and confined geometries.

Our simulations for the equilibrium properties of polymers agree well with both

previous computational studies and theoretical expressions from the renormalization group

calculations. We were able to recover the correct scalings for static properties, such as

radius of gyration and end-to-end vector, and also dynamic properties like the center of

mass diffusion coefficient. We have shown that the simulations can be greatly accelerated

by using a coarser description on the flow scale i.e. decreasing the bond length between

two beads compared to the lattice unit. However, as we have pointed out this comes at

the expense of the accuracy of hydrodynamic interactions. Nevertheless, it was shown that

for large polymers, for which the Schmidt number is large enough on the polymer scale,

the numerical method can still produce the correct hydrodynamic behavior. This allows

the use of long polymers where high resolution of different internal relaxation times of the

polymer is needed.

We have also investigated the flow of polymers through small channels where the

polymers are driven by either a unidirectional shear field, external forces, or a combination

of both. We have shown that a polymer in such conditions migrate thereby creating a

nonuniform center of mass distribution across the channel. Our findings in pressure driven

polymers agree with recent Brownian dynamics simulations [36] and a kinetic theory

for a dumbbell [57] and have supported the mechanism that the migration is a result of

.i-viiiii. Ii ic hydrodynamic interactions with the wall. In such flows we have shown that

migration can be towards the boundaries rather than the center in very narrow channels

due to the screening of hydrodynamic interactions. We were also able to clarify the correct

source of migration mechanism by simulations of a simple shear flow. We proved that

hydrodynamic interactions are the dominant mechanism for migration rather than the

entropic driving forces. Our results indicates that the idea that migration of polymers

greatly effect the macroscopic transport coefficients such as dispersion and elution rates.

Our simulations with external fields, in the slit geometry, have revealed some

important new physics. The polymers under external fields also migrate giving rise to

a nonuniform distribution. Our results and theoretical analysis have concluded that

although the fundamental reason for this migration is again hydrodynamic interactions

with boundaries, however the motion created by the external field differ fundamentally

from the shear flow. Hydrodynamic dispersion was also found to be a important part of

this migration mechanism in the limit of infinite Peclet numbers.

We have shown that combination of the external fields and unidirectional shear

flows give rise to both interesting fundamental physics and also greater control over both

the direction and magnitude of migration. This also results in a greater control over

the macroscopic transport properties. Size based separation in such systems are now

a possibility due to the highly reduced Taylor dispersion coefficients compared to pure

pressure driven flows.

This work enhances our understanding of fundamental physics involved in suspensions

of polymers in confined geometries. The use of a new and flexible algorithm opens the

doors for a wide variety of new problems previously prohibited by traditional algorithms.

These include, but are not limited to, flow of polymers in irregular geometries where

lattice Boltzmann offers inherent hydrodynamics as opposed to complicated calculations

of Green functions. Despite its drawbacks such as unnecessary fluid computations in very

dilute systems, and lattice artifacts for the mobility, lattice Boltzmann proves to be a very

efficient, accurate, competitive and most importantly flexible method for simulating

polymer systems. The improvement of the solid fluid coupling and incorporation

of charges are the most important two challenges to improve both the accuracy and

applicability of the method.


Numerical simulations show that the mean kinetic energy of a single monomer is

different from what would be expected from the level of the stress fluctuations in the fluid

(Fig. 3-3). Nevertheless, the fluctuation-dissipation relation between the effective friction

and diffusion coefficients holds within statistical errors (Fig. 3-1). It has been shown[21]

that fluctuations in a discrete system can be different from those in continuous systems

(c.f. Eqs. 2-15 and 2-16), and in order to understand why deviations from classical

equipartition do not affect the long-time dynamics, we here examine a simple model with

constant friction. The dynamics are described by the classical Langevin equation,


with a variance in the random force (R(t)R(t')) = 2T6(t t'). With this choice of random

force, Eq. A- 1 satisfies both the Stokes-Einstein relation D = T/l and equipartition

m (U2) = T. The real situation is more complex [20, 33], with a time-dependent friction

coefficient due to the temporal and spatial development of the hydrodynamic flow field.

Nevertheless, this simple model is sufficient to illustrate the deviations in mean kinetic

energy from classical equipartition.

The discrete equivalent of Eq. A- 1 can be written as

U(t + At)- (1 )U(t) + AU(t), (A-2)

where Q At/m and (AU(t)AU(t')) (2QT/m)6t,t,. The level of fluctuation in the

random velocity, AU, has been calculated by integrating the random force over the time

step At,[78]
AU(t) 1 t+At
A-(t) R(t')dt'. (A-3)
m Jt

After n time steps, Eq. A-2 has the solution

U(t + nAt) = (1 Q)U(t) + AU(t + kAt)(1 Q)(n-l-k), (A-4)

with the usual stability condition Q < 2. At long times, n > 1, the particle temperature is

related to the fluid temperature by

m (U2) = 2TO Z[(1 )21 k -T/2, (A-5)
tO 1 Q/2'

and violates the classical equipartition of energy. On the other hand, the diffusion

coefficient, which can be calculated from the discrete equivalent of the Green-Kubo

D = (U(0)U(0)) + At (U(nAt) U(0)) (A-6)
nl 1
satisfies the Stokes-Einstein relation. Thus we recover the correct long-time diffusive

behavior, despite the discrete time step, whereas the kinetic energy has an additional

factor (1 Q/2)-1. An accurate equipartition of energy requires that Q < 1, which

in our context means that the mass of the monomer must be large. However, numerical

calculations with Q < 0.5 show that the diffusive behavior of interacting particles is

unaffected by violations of equipartition.

An implicit update of Eq. A-1,

U(t + At) U(t)- QU(t + At) + AU(t), (A-7)

is often used to circumvent the stability constraint Q < 2 on the explicit solution

(Eq. A-2). It is straightforward to show that the Stokes-Einstein relation D = T/l is

again satisfied, but in this case the mean kinetic energy is given by

M KU2) =T (A-8)

The comparison of Eqs. A-5 and A-8 with the numerical simulations is shown in

Fig. 3-3. It can be seen that the simplified model presented here, which neglects the

time-dependence of the hydrodynamic interactions, still predicts the deviations from

equipartition quantitatively.


In order to test the fluctuations in the lattice Boltzmann scheme, we first derive

the correlation of momentum change on a particle. In Lattice Boltzmann the forces on

a particle are impulsive, instantaneous changes of momentum. Therefore, in order to

calculate the momentum change correlation, we start with the following equation,

pAt At
< APAP >= < F(t)F(t') >dt'dt (B-1)
2At J (1 ) < F(t)F(0) > dt (B-2)

To calculate the correlation function < F(t)F(0) > we follow the procedure mentioned

by [33]. We define Cf(t) =< F(t)F(0) >, the force autocorrelation function and

define Cp(w) as the power spectrum of F which is the Fourier transform of the force

autocorrelation function. We can now write

1 f
CF(t) = e-wtCp(w)dw, (B-3)

then equation B- 1 becomes

< APAP >- -. j CF(w)e- t( )dwdt. (B4)

We first start with the integration over t

=-> -it(1 t)dt= + (1 -e-^t (B-5)
o At Mw w2At
We now have to check if this quantity diverges as in the limit of small w (w -+ 0). To do

this we Taylor expand the exponential term which gives

1 1 (1 iAtw i2+t2w2
1 + 1. 2! (B-6)
S--. + O(w3). (B-7)
2 2

since there are no singularities we can nor rewrite this quantity in the following form
to simplify the calculation,

1 1- cos(wAt) + i sin(wAt)

wAt sin(wAt) 1 cos(wAt)
iw2A iw2A + w2AtB9)

The integral due to the total of the first two terms vanish and it remains only
to interpret the integral of the last term. Here we turn to the definition of the power
spectrum as given by [33]:

CF(w) kT[7(w) + (-w)] (B-10)

7(w) 6ij6wrlR[1 + R(-iw/v)112 (iR2/9v)]. (B- 11)

Here, 7(w) is what we call a friction kernel and R is the radius of the sphere. The
final result for the force spectrum function in frequency domain is found out to be,

CF(w) = ;i.iRkT [l1 + R ] (B-12)

Substituting this function we find,

< APAP >- :;:.-,/. TAt(R + R2) (B-13)


We first start with calculating the friction kernel for a plane wall that is sheared with

a oscillating velocity u0oet in one direction, in an infinite medium and we assume that

u = 0 at infinity from the wall.

We begin by solving the N ,1., r-Stokes equations with the above mentioned boundary

conditions. We find the velocity in x direction in the fluid as a function of y as

UX(y) Uoe-vy -v^ \ (C-t)

It now remains to calculate the shear stress aux at the wall(y = 0).

u(y =0) PW= (1 + )wt. (C-2)

The friction is the part that is proportional to u,

w)= (1+ i) (C-3)

Thus the friction kernel per unit area, the power spectrum of the forces, can be

written as

CF(w) kT [(w) + 7(-w)] (C-4)

CF(w) = kT 2pw. (C-5)

We then follow the same procedure for the sphere (see Appendix B) to calculate

< APAP > and the resulting equation for the plane wall is,

< APAP >- 4kT TAt-i (C-6)

Again it's important to note this result should be integrated once more over the

surface of interest.


The momentum correlation on a plane wall can also be calculated from the

information on the population densities and the use of bounce-back rules in lattice

Boltzmann method. We first have to identify the population densities hitting the wall

and bouncing back. In our model these are n3, n7, n8, n15, n17. The population density

n3 is in the normal direction (y) to the wall whereas n7 and ns have y and x (in opposite

directions) components and n15 and n17 have y and z (in opposite directions) components.

The momentum change on a stationary wall in the longitudinal directions can be

written as follows:

AP, 2(n n8)At (D-1)

APy 2(n15 17)At (D-2)

The population densities are given by the post collision distribution

S j- Ci (puu + [ineq') : (cc c)

since there is no flow and only fluctuations, puu can be neglected which leaves us with [neq'" (cici C'2)
ni* ac[p ci + c ] (D-4)

Here we will calculate the momentum change in x direction. The populations contributing

to this direction have a velocity of 4/2 thus have the weight factor a = The speed of

sound c = The post collision momentum flux [ineq' setting A = -1 is only the random

stress term thus, [Ineq' [If. These give

n7 = a[p + ] + + + ] (D-5)
4C 4 12 24
ns = a [p + ] + [- + ] (D-6)
1 4 t2 24

These give

< AP>2 > < > + < n2 > (D-7)

Since < j .j >= pAx3kT [2] and < j* > is a third of this, and using the equation 2-31, and

that p = 36 in our model, we find.

< AP2 > -kT + 12kT (D-8)

This is the momentum change correlation in x direction per unit area, and the same

expression is found for the z direction as well.


[1] C. W. J. Beenakker and P. Mazur, "Self-diffusion of spheres in a concentrated
suspension," Ph;-,:.i. A, vol. 120, pp. 388-410, 1983.

[2] A. J. C. Ladd and R. Verberg, "Lattice-Boltzmann simulations of particle-fluid
suspensions," J. Stat. Phys., vol. 56, pp. 286-311, 2001.

[3] A. J. C. Ladd, "Short-time motion of colloidal particles: Numerical simulation via a
fluctuating lattice-Boltzmann equation," Phys. Rev. Lett., vol. 70, pp. 1339, 1993.

[4] A. J. C. Ladd, "Numerical simulations of particulate suspensions via a discretized
Boltzmann equation Part I. Theoretical foundation," J. Fluid Mech., vol. 271, pp.
285, 1994.

[5] A. J. C. Ladd, "Hydrodynamic interactions in a suspension of spherical particles," J.
C7. in, Phys., vol. 88, pp. 5051, 1988.

[6] G. K. Batchelor, An Introduction to Fluid D,,,,i.. Cambridge University Press,
Cambridge, 1967.

[7] J. Happel and H. Brenner, Low-RB ,.,. ..1-1 Number H..l,.,1r.ii.,,*. Martinus Nijhoff,
Dordrecht, 1986.

[8] H. J. Masliyah and S. Bhattacharjee, Electrokinetic and Colloid Transport Phenom-
ena, Wiley, New York, 2006.

[9] A. J. C. Ladd, M. E. Colvin, and D. Frenkel, "Application of lattice-gas cellular
automata to the Brownian motion of solids in suspension," Phys. Rev. Lett., vol. 60,
pp. 975, 1988.

[10] P. Mazur and W. V. Saarloos, \ ,i y-sphere hydrodymamic interactions and
mobilities in a suspension," P,;;,'... A, vol. 115, pp. 21-57, 1982.

[11] P. Ahlrichs and B. Diinweg, "Simulation of a single polymer chain in solution by
combining lattice Boltzmann and molecular dynamics," J. C, ,,in Phys., vol. 111, no.
17, pp. 8225-8239, 1999.

[12] B. J. Alder and T. E. Wainwright, "Decay of the velocity autocorrelation function,"
Phys. Rev. A, vol. 1, no. 1, pp. 18-21, 1970.

[13] T. Phung, J. Brady, and G. Bossis, "Stokesian dynamcis simulation of brownian
suspensions," J. Fluid Mech., vol. 313, pp. 181-207, 1996.

[14] J. P. Boon and S. Yip, Molecular H. ,h..1 /,., ,i,:. Mc-Graw Hill, New York, 1980.

[15] A. J. C. Ladd, "Numerical simulations of particulate suspensions via a discretized
Boltzmann equation Part II. Numerical results," J. Fluid Mech., vol. 271, pp. 311,

[16] N. Q. Nguyen and A. J. C. Ladd, "Lubrication corrections for lattice-boltzmann
simulations of particle suspensions," Phys. Rev. E, vol. 66, no. 046708, 2002.

[17] 0. B. Usta, A. J. C. Ladd, and J. E. Butler, "Lattice-Boltzmann simulations of the
dynamics of polymer solutions in periodic and confined geometries," J. Ch. im Phys.,
vol. 122, pp. 094902-1, 2005.

[18] S. Succi, The Lattice Botlzmann Equation for Fluid D:in,,.:,- and B. :1 .,.l Oxford
University Press, 2001.

[19] K. Huang, Statistical Mechanics, Wiley, New York, 1987.

[20] L. D. Landau and E. M. Lifshitz, Fluid Mechnaics, Oxford University Press, New
York, 1966.

[21] A. J. C. Ladd and R. Verberg, "Lattice-Boltzmann simulations of particle-fluid
suspensions," J. Stat. Phys., vol. 104, pp. 1191-1251, 2001.

[22] A. J. C. Ladd, H. Gang, J. X. Zhu, and D. A. Weitz, "Time-dependent collective
diffusion of colloidal particles," Phys. Rev. Lett., vol. 74, pp. 318, 1995.

[23] A. J. C. Ladd, "Hydrodynamic screening in sedimenting suspensions of non-Brownian
spheres," Phys. Rev. Lett., vol. 76, pp. 1392, 1996.

[24] A. J. C. Ladd, "Effects of container walls on the velocity fluctuations of sedimenting
spheres," Phys. Rev. Lett., vol. 88, pp. 0 !- ;il, 2002.

[25] C. K. Aidun, Y. N. Lu, and E. Ding, "Direct ,in i,.i -i of particulate suspensions with
inertia using the discrete Boltzmann equation," J. Fluid Mech., vol. 373, pp. 287-311,

[26] F. M. Auzerais, J. Dunsmuir, B. B. Ferreol, N. Martys, J. Olson, T. S. Ramakrishnan,
D. H. Rothman, and L. M. Schwartz, "Transport in sandstone: A study based on
three dimensional microtomography," G. *'i i- Res. Lett., vol. 23, pp. 705-108, 1996.

[27] N. S. Martys and H. C'!. i,1 "Simulation of multicomponent fluids in complex
three-dimensional geometries by the lattice boltzmann method," Phys. Rev. E, vol.
53, no. 1, pp. 743-750, Jan 1996.

[28] L. W. P. Manz, B; Gladden, "Flow and dispersion in porous media:
Lattice-boltzmann and nmr studies," Aiche J., vol. 45, pp. 1845-1854, 1999.

[29] P. Ahlrichs, R. Everaers, and B. Diinweg, "Simulation of a single polymer chain in
solution by combining lattice Boltzmann and molecular dynamics," Phys. Rev. E, vol.
64, pp. 040501(R), 2001.

[30] U. Frisch, D. d'Humikres, B. Hasslacher, P. Lallemand, Y. Pomeau, and J.-P. Rivet,
"Lattice gas hydrodynamics in two and three dimensions," Complex S1-/4. n- vol. 1,
pp. 649, 1987.

[31] D. L. Ermak and J. A. McCammon, "Brownian dynamics with hydrodynamic
interactions," J. 0C1 iin Phys., vol. 69, pp. 1352, 1978.

[32] X. Fan, N. Phan-Thien, N. T. Yong, X. Wu, and D. Xu, \!icrochannel flow of a
macromolecular suspension," Phys. Fluids, vol. 15, no. 1, pp. 11-21, 2003.

[33] E. H. Hauge and A. Martin-Loef, "Fluctuating hydrodynamics and brownian
motion," Journal of Statistical Ph'-. vol. 7, no. 3, pp. 259-281, 1972.

[34] A. J. C. Ladd, "Dynamical simulations of sedimenting spheres," Phys. Fluids A, vol.
5, pp. 299, 1993.

[35] R. Jendi i 1: M. Graham, and J. dePablo, "Effect of confinement on DNA dynamics
in microfluidic devices," J. Ch( in Phys, vol. 119, pp. 1165-1173, 2003.

[36] R. M. Jendi. ..1.: D. C. Schwartz, J. J. de Pablo, and M. D. Graham, "Shear-induced
migration in flowing polymer solutions: Simulation of long-chain DNA in
microchannels," J. C(' ,in Phys., vol. 120, pp. 2513-2529, 2004.

[37] M. Fixman, "Construction of the Langevin forces in the simulation of hydrodynamic
interaction," Macromolecules, vol. 19, pp. 1204-1207, 1986.

[38] R. Jendi. i, .: M. Graham, and J. dePablo, "Hydrodynamic interactions in long
chain polymers: Application of the chebyshev polynomial approximation in stochastic
simulations," J. Ch'. il Phys, vol. 113, pp. 2894-2900, 2000.

[39] N. J. Woo, E. Shaqfeh, and B. Khomami, "Effect of confinement on dynamics and
rheology of dilute DNA solutions. I. Entropic spring force under confinement and
numerical algorithm," J. Rheol, vol. 48, pp. 281-298, 2004.

[40] N. J. Woo, E. Shaqfeh, and B. Khomami, "Effect of confinement on dynamics and
rheology of dilute DNA solutions. II. Effective rheology and single chain dynamics,"
J. Rheol, vol. 48, pp. 299-318, 2004.

[41] C.-C. Hsieh and R. G. Larson, "Modeling hydrodynamic interaction in Brownian
dynamics: Simulations of extensional and shear flows of dilute solutions of high
molecular weight polystyrene," J. Stat. Phys, vol. 36, pp. 717, 1984.

[42] Y. Oono and M. Kohmoto, "Renormalization group theory of transport properties of
polymer solutions. i. Dilute solutions," J. Ch'. in Phys., vol. 78, pp. 520-528, 1983.

[43] H. Hasimoto, "On the periodic fundamental solutions of the Stokes equations and
their application to viscous flow past a cubic array of spheres," J. Fluid Mech., vol. 5,
pp. 317, 1959.

[44] A. J. C. Ladd, "Hydrodynamic transport coefficients of random dispersions of hard
spheres," J. Ch'. i, Phys., vol. 93, pp. 3484, 1990.

[45] T. M. Squires and M. P. Brenner, "Like-charge attraction and hydrodynamic
interaction," Phys. Rev. Lett., vol. 85, no. 23, pp. 4976-4979, Dec 2000.

[46] F. Brochard and P. G. de Gennes, "Dynamics of confined polymers," J. Ch, in, Phys.,
vol. 67, pp. 52, 1977.

[47] H. J. L. and D. M., "Diffusion of macromolecules in narrow capillaries," j. Phys.
Ch. i, vol. 96, pp. 4046-4052, 1992.

[48] S. Asakura and F. Oosowa, "On interaction between two bodies immersed in a
solution of macromolecules," J. C. iin, Phys., vol. 22, pp. 1225-1256, 1954.

[49] J. H. Aubert and M. Tirrel, "Effective viscosity of dilute polymer solutions near
confining boundaries," J. Ch in, Phys., vol. 77, pp. 553-561, 1982.

[50] D. Aussere, H. Hervet, and F. Rondelez, "Concentration dependence of the interfacial
depletion 1,-r thickness for polymer solutions in contact with nonadsorbing walls,"
Macromolecules, vol. 19, pp. 85-88, 1986.

[51] W. F. Busse, "Two decades of high-polymer physics a survey and forecast," Phys.
T ./.i vol. 17, pp. 32, 1964.

[52] A. B. Metzner, Y. Cohen, and C. Rangel-Nafaile, "Inhomogeneous flows of
non-Newtonian fluids: Generation of spatial concentration gradients," J. Non-
Newtonian Fluid Mech., vol. 5, pp. 449-462, 1979.

[53] U. S. Agarwall, A. Dutta, and R. A. Mashelkar, \! i oi ,ii of macromolecules under
flow: The physical origin and engineering applications," C.i in, Eng. Sci., vol. 49, pp.
1693-1717, 1994.

[54] K. A. Dill and B. H. Zimm, "Rheological separator for very large DNA molecules,"
Nucleic Acids Research, vol. 7, pp. 735-749, 1979.

[55] L. Fang, H. Hu, and R. G. Larson, "DNA configurations and concentration in
shearing flow near a glass surface in a microchannel," J. Rheol., vol. 49, no. 1, pp.
127 -138, 2005.

[56] M. S. Jhon and K. F. Freed, "Polymer migration in Newtonian fluids," J.Polym.
Sci.Part B, vol. 23, pp. 955-971, 1985.

[57] H. Ma and M. Graham, "Theory of shear induced migration in dilute polymer
solutions near solid boundaries," Phys. Fluids., vol. 17, pp. 083103, 2005.

[58] M. Fixman, "Effects of fluctuating hydrodynamic interaction," J. Ch.l i Phys, vol.
78, pp. 1594, 1983.

[59] N. Liron and S. Mochon, "Stokes flow for a Stokeslet between two parallel flat
plates," J. Eng. Math., vol. 10, pp. 287 303, 1976.

[60] L. C. Nitsche and E. J. Hinch, "Shear induced lateral migration of Brownian rigid
rods in parabolic channel flow," J. Fluid Mech., vol. 332, pp. 1-21, 97.

[61] R. L. Schiek and E. S. G. Shaqfeh, "Cross-streamline migration of slender Brownian
fibres in plane Poiseuille flow," J. Fluid. Mech, vol. 332, pp. 23-39, 1997.

[62] R. G. Larson, "The rheology of dilute solutions of flexible polymers: Progress and
problems," J. Rheol., vol. 49, no. 1, pp. 1-70, 2005.

[63] 0. B. Usta, J. E. Butler, and A. J. C. Ladd, "Flow-induced migration of polymers in
dilute solution," Phys. Fluids., vol. 18, pp. 031703, 2006.

[64] J.-L. Viovy, "Electrophoresis of dna and other polyelectrolytes: Physical
mechanisms," Rev. Mod. Phys., vol. 72, no. 3, pp. 813-872, Jul 2000.

[65] Z. C'!i i, and A. C('! ,111, "DNA separation by EFFF in a microchannel ," J. Coll.
Int. Sci., vol. 2005.

[66] J. Zheng and E. S. Yeung, "Anomalous radial migration of single DNA molecules in
capillary electrophoresis," Anal. Ch.' vol. 74, pp. 4536-4547, 2002.

[67] J. Zheng and E. S. Yeung, \!. I, Ii ,,-in for the separation of large molecules based on
radial migration in capillary electrophoresis," Anal. Ch. i, vol. 75, pp. 3675, 2003.

[68] D. Long and A. Ajdari, "A note on the screening of hydrodynamic interactions, in
electrophoresis, and in porous media," Euro Phys. J. E, vol. 4, pp. 29-32, 2001.

[69] R. E. Caflisch and J. H. C. Luke, "Variance in the sedimentation speed of a
suspension," Phys. Fluids, vol. 28, pp. 759, 1985.

[70] J. Martin, N. Rakotomalala, and D. Salin, "Hydrodynamic Dispersion of Noncolloidal
Suspensions: Measurement from Einstein's Argument," Phys. Rev. Lett., vol. 74, pp.
1347, 1995.

[71] W. B. Russel, E. J. Hinch, L. G. Leal, and G. Tiefenbruck, "Rods falling near a
vertical wall," J. Fluid. Mech, vol. 83, pp. 273-287, 1977.

[72] T. N. Swaminathan, K. Mukundakrishnan, and H. Hu, "Sedimentation of an ellipsoid
inside an infinitely long tube at low and intermediate reynolds numbers," J. Fluid
Mech., vol. 551, pp. 357 385, 2006.

[73] G. S. M ,iii,,- "Limiting laws and counterion condensation in polylelectrolyte
solutions 7. electrophoretic mobility and conductance," J. Phys. Ch. i11, vol. 85, pp.
1506, 1981.

[74] J. Barrat and J. Joanny, "Theory of polyelectrolyte solutions," Adv. Ci. ,11 Phys.,
vol. 94, pp. 1-66, 1994.

Full Text








IwouldliketothankmyadvisorsJasonE.ButlerandAnthonyJ.C.Laddfortheirprofessionalism,tremendousknowledge,patience,supportandhelpduringthelastveyears.Ihavedisagreedwiththemonmanyoccasionsandstilldo.HoweverweallagreethatwecandisagreewhichIfoundoneofthemostvaluablethingsduringmydoctoraldegree.However,Idowishmydisagreementstohavebeenonmoresolidgrounds.Iamyettocompletelydisproveanythingtheyhavetoldmeexceptforsmallissues.Itwasagreatpleasureandalsoadistinctprivilegetohaveworkedwithtwosuchqualitypeoplesodierentyetsosimilar.IamgratefultoDr.JamesDuftynotonlyforhispresenceonmycommitteebutalsoforhislectures,kindness,sincerity,wisdomandgreatexampleofhowascientistshouldbe.IwouldalsoliketothankmyothercommitteemembersDr.RichardDickinson,andDr.SergeiObukhovandDr.Chang-WongParkfortheirtimeandadvice.IamalsogratefultomyundergraduateadvisorDr.TurkanHalilogluwhostillprovidesahervaluablesupportandfriendship.IwouldnowliketoturntomycolleaguesDr.PiotrSzymczak,Dr.JonghoonLeeandDr.ByoungjinChunforbeingmembersoftheoriginallunchgang.AsJonghoonputitinhisthesis"Wediscussedsubjectswithoutboundaries".Piotrhasprobablybeenoneofthemostinuentialpersoninmylifebesidesmyfamily,myfriendsandmyadvisors.IwishonedayIcanbehalfasgoodasheisinbeingabletoreducecomplexproblemsintosimpleterms,amongalltheotherthingshedoessowell.JonghoonandByoungjinhavebothbeengreatinspirationtomebytheircalmness,dedicationandhardwork.Iwillneverforgetneithertheirpersonalnortheirprofessionalhelpandsupport.IwouldalsoliketothankJonathanM.BrickerformanynightswehavespenttogetherinthelabsharingnotonlysciencebutourlivesandQuyenNguyenforhelpingmebothpersonallyandprofessionallywhenIrstjoinedthegroup.IwouldalsoliketoacknowledgeGauravMisrawho,withinthelimitedtimewekneweachother,impressedmegreatlywithhisfocusandloveforscience.IwouldalsoliketothankJoontaekPark,Hyun-OkPark,PhillipD.Cobb,MuraliRangarajan,andDazhiYu. 4




page ACKNOWLEDGMENTS ................................. 4 LISTOFFIGURES .................................... 8 ABSTRACT ........................................ 12 CHAPTER 1COMPLEXFLUIDSANDHYDRODYNAMICS ................. 14 1.1Introduction ................................... 14 1.2FluidMotionandHydrodynamicEquations ................. 15 1.3HydrodynamicInteractions ........................... 17 1.4NumericalMethodstoSolveHydrodynamicEquations ........... 19 2LATTICE-BOLTZMANNMETHODSFORCOLLOIDALANDPOLYMERICSUSPENSIONS .................................... 24 2.1BoltzmannEquation .............................. 24 2.2LatticeBoltzmannMethod ........................... 26 2.2.1AnOverviewofLattice-BoltzmannMethod .............. 26 2.2.2LatticeBoltzmannEquation ...................... 28 2.3Fluctuations ................................... 30 2.4SolidParticles .................................. 31 2.4.1Finite-SizeSphericalParticles ..................... 31 2.4.2PointParticles .............................. 33 2.4.3ParticleMotion ............................. 36 2.5VericationoftheFluctuations ........................ 37 3DYNAMICSOFPOLYMERSOLUTIONSINPERIODICANDCONFINEDGEOMETRIES .................................... 43 3.1Introduction ................................... 43 3.2Results ...................................... 45 3.2.1Fluctuation-Dissipation ......................... 45 3.2.2HydrodynamicInteractions ....................... 50 3.2.3DiusionofanIsolatedPolymer .................... 52 3.2.4PolymerDiusioninaChannel .................... 57 3.3Conclusions ................................... 59 4FLOW-INDUCEDMIGRATIONOFPOLYMERSINDILUTESOLUTION .. 62 4.1Introduction ................................... 62 4.2Method ..................................... 63 4.3ResultsandDiscussion ............................. 65 4.4Conclusion .................................... 73 6


... 75 5.1Introduction ................................... 75 5.2Methods ..................................... 76 5.3Results ...................................... 76 5.4Conclusion .................................... 85 6CONCLUSIONS ................................... 86 APPENDIX ADISCRETELANGEVINEQUATION ....................... 88 BTHEMOMENTUMCHANGECORRELATIONONASPHERE ........ 91 CTHEMOMENTUMCHANGECORRELATIONONAPLANEWALL,CONTINUUMTHEORY ....................................... 93 DTHEMOMENTUMCHANGECORRELATIONONPLANEWALL,DISCRETETHEORY ....................................... 94 REFERENCES ....................................... 96 BIOGRAPHICALSKETCH ................................ 102 7


Figure page 2-1Locationoftheboundarynodesforacircularobjectofradius2.5latticespacings.Thevelocitiesalonglinkscuttingtheboundarysurfaceareindicatedbyarrows.Thelocationsoftheboundarynodesareshownbysolidsquaresandthelatticenodesbysolidcircles. ................................. 32 2-2Momentumchangecorrelationsonaplanewall.P(t)isthetotalmomentumexchangebetweentheuidandtheparticleatanyinstantt.Resultsfromsimulationsarecomparedwithatheoreticalexpressionderivedfromcontinuumtheory.Thetheoreticalexpressionishigherthanthesimulationdataatthersttime-step,erroraccumulatesforacoupleoftime-stepsandthenthedierencebetweenthesetwocurvesturnintoaconstantoset.Therstpointforthiscorrelationiscalculatedcorrectly(1%)fromadiscretetheorystartingwiththeuctuationsinpopulationdensities(section D ). ......................... 39 2-3Momentumexchangecorrelationsonaninnitelymassivesphere.P(t)isthetotalmomentumexchangebetweentheuidandtheparticleatanyinstantt.Resultsfromsimulationsarecomparedwithatheoreticalexpressionderivedfromcontinuumtheory. ................................ 40 2-4VericationofOnsager'slinearresponsetheory.Thedecayofthevelocitycorrelationsiscomparedwiththevelocitydecayofaparticlewithaninitialgivenvelocity.Thetwocurvesleadtothesamediusioncoecientuponintegration. ..... 41 2-5Decayofthevelocityautocorrelationfunctionforasingleparticlewithmassfactorof1&10andaradiusof5.4,normalizedbytheequipartitionvalueU2=kT=mB.Thetwodierentcurvescorrespondtodierentintegrationschemes.Thesetwocurveseventuallyleadtothesamediusioncoecient. ........ 42 3-1Relationbetweentheeectivefrictioncoecientandtheinputfriction0(Eq. 2{21 ).Theeectivefrictioncoecienthasbeenobtainedfromthedragforceonasinglemonomer(circles)andthediusionofthemonomerinathermallyuctuatinguid(squares).Theslopeofthelineisunity,indicatingaconstantosetbetween1and10.TheunitcellsizeL=20x. ..................... 46 3-2Osetbetween1and10asafunctionofsystemvolume.Thevariationinthedimensionlessmobilityg(Eq. 3{1 )isshownforacubiccelloflength,L.Resultsareshownforinputhydrodynamicradiia0=0:32(circles)anda0=0:15(squares).Themonomerfriction0=6a0. ................. 48 3-3Meankineticenergyofasinglemonomerasafunctionofthedimensionlessfrictionalcoupling,t=m.Thecalculatedkineticenergyforimplicit(circles)andexplicit(squares)updatesiscomparedwiththemodeldescribedinAppendix A (lines). 49 8


................................... 51 3-5Comparisonoftheparallelcomponentofthetwoparticlemobilitytensor,k12.Resultsobtainedbydissipative(Eq. 3{10 ,circles)anductuating(Eq. 3{11 ,squares)simulationsarecomparedwithanindependentsolutionoftheperiodicStokesequations(Eq. 3{9 ,triangles). ........................ 52 3-6Comparisonoftheperpendicularcomponentofthetwoparticlemobilitytensor,?12.Resultsobtainedbydissipative(Eq. 3{10 ,circles)anductuating(Eq. 3{11 ,squares)simulationsarecomparedwithanindependentsolutionoftheperiodicStokesequations(Eq. 3{9 ,triangles). ........................ 53 3-7End-to-enddistanceandradiusofgyrationinaperiodicunitcell.ThescalingofRe(upperpoints)andRg(lowerpoints)isshownfordierentmonomerseparations,b,andtemperatures,T.Thedashedlinesarethetheoreticalscalingsforanexcludedvolumechain,withexponent=0:59. ....................... 54 3-8Centerofmassdiusioncoecientasafunctionofchainlengthindierentsizecells.Thediusioncoecientsarenormalizedbythemonomerdiusioncoecient,Dm=T=,andareshownforperiodicsystemswithlengthsL3:5RgandL7Rg.Inallcasesthemeanseparationbetweenthemonomersb=xandthetemperatureT=0:1.Resultsareshownfortherawdiusioncoecient(opensymbols),aswellasthecorrectedvalues(lledsymbols). ............. 55 3-9Centerofmassdiusioncoecient,withnite-sizecorrections,fordierentvaluesofbandT.ThelengthoftheperiodiccellL7Rgineachcase.Thestatisticaluncertaintiesinthediusioncoecientaresmallerthanthesizeofthesymbols,exceptforthelongestchain,N=1024,wheretheyareshownexplicitly. ..... 56 3-10Centerofmassdiusioncoecientsparallel(Dk)andperpendicular(D?)tothechannelwallsasafunctionofchannelwidthH;theradiusofgyrationofthepolymer(N=128,b=1)isapproximately7x.Theconneddiusioncoecients,normalizedbythemonomerdiusivity,areshownforasquarecell,L,withvariousratiosofLtoH. .......................... 58 3-11Centerofmassdiusioncoecientforaconnedpolymerrelativetotheunconneddiusivity,Dk=D.ThesolidlineindicatestheFaxencorrection[ 7 ]forweakconnement,assumingthepolymerliesinthecenterofthechannel.Thedashedlineisthetheoretical(Rg=H)2=3scaling[ 46 ].Theaspectratioofthecell,L=H,wassetto4ineachcase. ................................... 60 9


...... 66 4-2CenterofmassdistributionforuniformshearowasafunctionofPecletnumberandchannelwidth.ResultsareshownforN=16.ThedistributionsarenormalizedsuchthatRH=Rg0P(y=Rg)d(y=Rg)=1. ........................ 67 4-3Componentsoftheradiusofgyrationalongtheow,hx2i,andconnement,hy2i,axes.ResultsforchainsoflengthN=16areshownattwodierentPecletnumbers,normalizedbyRg=(hx2i+hy2i+hz2i)0:5. ............... 68 4-4Centerofmassdistributionfordierentlevelsofdiscretizationofthepolymerinatightlyconnedchannel,H=Rg=3.ThedistributionsarenormalizedsuchthatR10P(y=Rg)d(y=Rg)=1. ............................ 69 4-5Averageaxialvelocityofthecenterofmassofpolymersnormalizedbytheaverageuidow.TheresultsarepresentedwithrespecttothepolymerPecletnumberfordierentchannelsizes. .............................. 70 4-6HydrodynamicdispersioncoecientsofapolymerchaininPoiseuilleow.Numericalsimulationsareshownaslledsymbols,alongwithTaylor-Aristheory(TA)fortracerparticles(Ld=0,dashedline)andfornitesizedparticleswithaxeddepletionlayer(Ld=Rg=H=8,solidline).Theoreticalcalculations,accountingforthedepletionlayerthicknessdeterminedfromsimulation,areshownasopensymbols.Thestandarderrorinthesimulationdatacorrespondstothesizeofthesymbols. ...................................... 72 4-7Comparisonofdispersioncoecientswithdierentinertia.ThehighestparticleReynoldsnumberforkT=0:1isalmost3:4,whereasthechannelReynoldsnumbersareupto30forthesamesimulations. ....................... 73 4-8Comparisonofcenterofmassdistributionsat3dierenttemperaturesatPe=300.ThethreedierenttemperaturescorrespondtodierentReynoldsnumbers.ThehighestparticleReynoldsnumberatkT=0:1isaround3:4. ........ 74 5-1CenterofmassdistributionofasinglepolymerinachannelH=Rg=8undertheapplicationofanexternalbodyforce.Halfofthedistributionisshown,withtheboundaryaty=H=0andthecenterofthechannelaty=H=0:5. ...... 77 10


................... 79 5-3a)ScalingofthedistributionfunctionwithrespecttothechannelwidthunderconstantpolymerPecletnumber,Pe.b)ScalingofthedistributionfunctionwithrespecttopolymerPecletnumberusingdierentlevelsofchaindiscretization.Theboundaryisaty=H=0andthecenterofthechannelisaty=H=0:5. ... 81 5-4Centerofmassdistributionsforconcurrentapplicationofanexternalbodyforceandpressure-drivenow.Thesolidcurveshowsthelevelofmigrationunderthepressure-drivenowonly.TheowPecletnumberinallthreecasesisPef=12:5.Theboundaryisaty=H=0andthecenterofthechannelisaty=H=0:5. 82 5-5Centerofmassdistributionsforcountercurrentapplicationofanexternalbodyforceandpressure-drivenow.Thesolidcurveshowsthelevelofmigrationunderthepressure-drivenowonly.TheowPecletnumberinallcasesisPef=12:5.Theboundaryisaty=H=0andthecenterofthechannelisaty=H=0:5. ... 83 11






1 ]andconstitutesthefocusofthiswork.Thesehydrodynamicinteractionsarelongrangedandmany-body,thereforetheycreateagreatdealofcomplexityfortheoreticalcalculations.Inthisworkweinvestigatethestaticanddynamictransportpropertiesofhard-spherecolloidalsuspensionsanddilutepolymersolutionswithanumericalsimulationmethod.ThismethodcombinesNewtoniandynamicsofsolidparticleswithauctuatinglattice-Boltzmannuid[ 2 { 4 ].SuchanapproachleadstothesolutionofNavier-Stokesequations,withsuspendedparticles,onlargetimeandlengthscales.Thismethodiscloselyrelatedtoandisadescendentofthelattice-gascellularautomataforcolloidalsuspensions[ 5 ]. 14


1 ,Iwillbrieydiscussthehydrodynamicequationsgoverningthemotionofthesuspendinguid( 1.2 )alongwithadiscussionofthelowestorderdescriptionofhydrodynamicinteractionsbetweensuspendedsolidparticles( 1.3 ).Iwillalsopresentabriefhistoryandsummaryofnumericaltechniquesonthesubject.InChapter 2 ,Iwillpresenttheuctuatinglattice-Boltzmannmethod.Thisincludestwodierentapproachesforthesolidparticlescoupledtothelattice-Boltzmannuidforhard-spherecolloidalandpolymersuspensions.Iwillthenpresentstaticanddynamicproblems(Chapter 3 )involvingdilutepolymersolutions.ThisisfollowedbyanumericalinvestigationofthemigrationofpolymerchainsunderunidirectionalshearingowsandexternaleldsinChapters 4 and 5 6 ].Inphysicaltermsthisrequirestheratioofmeanfreepath()ofamoleculetothelineardimension(L)ofthevolumetobesmallsuchthat=L1.ThisratioisknownastheKnudsennumber. 15


@t=r(v);(1{1)whereisthelocaldensityoftheuidandvisthelocalmassaverageuidvelocity.Inthisstudywewillassumethattheuidisincompressible,inwhichcaseEq. 1{1 reducesto 7 ].Thereforethemomentumbalanceequationreads, 7 ].Thebulkviscosity,,canbeignoredforalowdensitymonatomicgas[ 7 ]andanincompressibleuid.SubstitutionofthisinformationandassumingincompressibilityinEq. 1{3 givestheNavier-Stokesequation: 16


,whereLisacharacteristicmacroscopiclength.ForsmallenoughReynoldsnumbers(Re1),theinertialtermscanbeneglected,andwiththecontinuityequationthisgivesrisetothecreepingmotionorStokesequationsofuidow[ 7 ]:r2v=rp 8 ],U(r)=U+(rR); 17


1{8 ,isv(r)=U+(rR)forjrRj=a; 9 ].Iftherelativevelocity(Uv(r))ofthesolidanduidphasesissmall,theeectsoftheuidcanbedescribedintermsofforcesandtorquesexertedonthesolidparticles[ 9 ].Theforces,F,andtorques,T,canbeexpressedintermsoftheparticlevelocities,U,and,throughcongurationdependentfrictiontensors,Fi=NXj=1[TTijUj+TRijj] (1{10)Ti=NXj=1[RTijUj+RRijj]: (1{12)i=NXj=1[RTijFj+RRijTj]: 10 ]whereadistributionofpointforcesonthesurfaceofeachsphereisassumedandthisdistributionisthenexpandedaboutthecenterofeachsphereinaseriesofsphericalharmonics[ 5 ].Theseinteractionsarelong-ranged,notpairwiseadditive,andoftenimpossibletocalculateanalyticallyexceptforspecialcasessuchaspairsofsuspended 18


11 ].Thisisagreatsimplicationandsuchconditionscanusuallybeassumedforapolymerchainthatconsistsofhydrodynamiccenters(beads)connectedbysprings.Theinteractionofsuchparticlescanbedescribedlocallyintermsoftherelativevelocitybetweentheuidandthesolidparticleandthecollectiveinuenceofparticlesbecomesadditiveinthecreepingowlimit.Eachparticleinthesuspensioninducesavelocityintheuidthatfollowsv(r)=Hf 8rhI+rr 12 ]anddierenceintheshear-ratedependenceofthehydrodynamicandBrowniancomponentsofstress[ 13 ].Computersimulationscanresolveproblemswhicharenotaccessiblebylaboratoryexperiments.Simulationsallowforinvestigationofidealsystemswhichhelpidentifysubtleindividualeects,oneofwhichishydrodynamicinteractions;thefocusofthisthesis.In 19


14 ].Althoughquitestraightforward,suchanalgorithmiscomputationallyveryineectivewhenoneisinterestedinthelargescaledynamicsofcolloidalparticleswhicharebiggerthanthesurroundinguidparticles.Thesimulationofsmalleruidparticlestakesupmostofthecomputationaltime,thereforeaverylimitedtimescalecanbesimulatedwithsuchanapproach.Thesizeseparationofsuspendedcolloidalparticlesandtheuidalsoresultsinatimescaleseparationforthemotionofthetwophases.Onecanexploitthesescaleseparationsanddescribetheuidasacontinuumwhichresultsinahugereductionofdegreesoffreedom.StokesiandynamicsandBrowniandynamicsaretwosuchmethodswherehydrodynamicinteractionsareassumedtobeinstantaneousduetotheabovementionedtime-scaleseparation.StokesiandynamicsisanaccuratesimulationtechniquewhichusesareformulationoftheStokesequations(Eq. 1{6 )asanintegralequationfortheuidvelocityeldgeneratedbyaprescribedforcedensityfind[ 2 ].Theforcedensityislocalizedtotheparticlesurfacesandforsphericalparticlesoneneedstoexpandtheforcedensityandresultingvelocity 20




1{18 canberewritten,withthesubstitutionforthehydrodynamicterm,whichyieldsaLangevintypeequation 22


2 .VericationtestsalongwithphysicalproblemswillbepresentedinChapters 3 4 ,and 5 23


4 ].Themostimportantfeaturesofthismethodarethenaturalevolutionofhydrodynamicinteractionsandalinearscalingofthecomputationalcostwithrespecttothenumberofparticles.However,itisimportanttonotethatthecomputationalcostinfactscalesdirectlywiththeuidvolume.ThisChapterincludesabriefreviewofthebasicconceptsofthelattice-Boltzmannmethod,startingwithadiscreteBoltzmannequationanddescriptionoftheeldsandthelattice.Theincorporationofuctuationsandthecouplingofsolidparticlestotheuidarealsodiscussed.ReadersinterestedinamoredetaileddiscussionshouldrefertopapersbyLaddandhisgroup[ 2 4 15 { 17 ]. 18 ].InthissectionIwillstartbydescribingsomekeypointsofthisequation.Theprobabilitydensityordistributionfunction,f(x;p;t),givestheprobabilityofndingamoleculearoundpositionxwithmomentumpattimet.Thisisthekeyobjectofthekinetictheory.TheBoltzmannequation,derivedin1872,describestheevolutionofthisdistributionfunction,f,intermsofmicrodynamicinteractions.Thisequationreadsasfollows. 24


18 ].Thecollisionoperatorinvolvesthetwobodydistributionfunctionf12expressingtheprobabilityofndingamolecule1aroundx1withvelocityv1andasecondmolecule2aroundx2withvelocityv2,bothattimet.Thedicultyinwritingf12isthatitcallsintoplaythethree-bodydistributionfunctionf123whichinturndependsonf1234andsoonaccordingtotheBBGKY(Bogoliubov,Born,Green,Kirkwood,Yvon)hierarchy[ 19 ].InordertoachieveaclosureinhisequationBoltzmannmadeafewstringentassumptionsonthenatureofthesystem( 2{1 );adilutegasofpoint-like,structurelessmoleculesinteractingviaashort-rangetwo-bodypotential.Undersuchconditionsintermolecularinteractionscanbedescribedsolely,intermsoflocalizedbinarycollisions,withmoleculesspendingmostoftheirlifespaninfreetrajectories(intheabsenceofexternalelds)perfectlyunawareofeachother.Withinthisapproximationwecansplitthecollisiontermintogainandlosscomponents: 25


18 ].Inaliquidthesituationiscompletelydierentduetomuchhigherdensity.Inthiscasetheparticlesareinconstantinteractionwithoneanother.ViolationsofBoltzmann'smolecularchaoscanoccurinrelativelydiluteuidsduetononlinearuid-particleinteractions.Oneimportantexampleisthelong-timetailsrstdetectedbyAlderandWainwright[ 12 ].Tosummarizealltheaboveinformation,wewritetheBoltzmannequationinitsnalform[ 18 ]. 2.2.1AnOverviewofLattice-BoltzmannMethodHistoricallyLatticeBoltzmannisthedescendentoftheLatticeGasAutomaton(LGA)wherethediscretepositionsandthevelocitiesofthesingleparticlesaretrackedonthelatticewhilesatisfyingtheconservationequationslocally.Latticegasmethodsprovedtobequitesimpleandusefulforhydrodynamics,thoughhadsomedrawbackssuchasthestatisticalnoisewhichrequireslongaveragingtimes.ThelatticeBoltzmannmethods,thoughtheycanbeindependentlyderivedfromtheBoltzmannequation,evolvedfromtheideaofxingthedrawbacksoftheLGAapproach.LatticeBoltzmannfollowstheevolutionofthedistributionfunctionfortheparticlesratherthantrackingindivdualuidparticles,whichresultsinpre-averaginginmomentum 26


20 ].Thisisachievedbyintroducingadiscreterandomstresstensorinagreementwithadiscreteuctuation-dissipationtheorem.ThecomputationalmethodIpresentinthisChapterisbasedonauctuatinglattice-Boltzmannmodeloftheuidphase,[ 3 4 21 ]whichincludesthermallydrivenuctuationsintheuidvelocityeldviarandomuctuationsaddedtothestresstensor.ThekeydierencewithBrowniandynamicsisthattheserandomuctuationsarelocalinspaceandtime,avoidinganyfactorizationofthecovariancematrix.Instead,initiallyuncorrelateductuationsdevelopinspaceandtimebyviscousdiusionofmomentum,givingrisetohydrodynamiccorrelationsatsucientlylargespatialandtemporalscales.Numericalsimulationshaveshownthatthedissipativeanductuatinghydrodynamicinteractionsbetweennite-sizeparticlescanbereproducedquantitatively[ 4 ].Computationally,thelattice-Boltzmannmethodisinherentlyandstraightforwardlyparallelizable,enablinglargesimulationstobespreadacrossaclusterofprocessorswitharelativelylownetworkbandwidth.Thelattice-Boltzmannmethodhasbeenusedforawidevarietyofcomplexows,includingsuspensionsofcolloidal[ 3 22 ]andnon-colloidal[ 23 { 25 ]particles,andporousmedia[ 26 { 28 ].Sincetheuidequationsaresolvedonagrid,externalboundariesandimposedoweldscanbeincludedatnoadditionalcomputationalcost.However,suspendedsolidparticlesoccupyasubstantialuidvolume,typicallyoftheorderof100gridpoints,andinteractwiththeuidthroughboundarynodesdistributedontheparticlesurface.Sincethecomputationtimeofthemethodisproportionaltotheuidvolume,thesuspensionmodelisunnecessarilyexpensiveforpolymersolutions,wherethemonomerscanbemodeledaspointforcescoupledtotheuidbyafrictioncoecient,.AhlrichsandDunwegintegratedthisfrictionalcouplingintotheuctuatinglattice-Boltzmann 27


11 29 ]. Thetimeevolutionofni(r;t)isdescribedbyadiscreteanalogueoftheBoltzmannequation,[ 30 ] 3;a1=1 18;ap 36:(2{7)Inthesesimulationsweusea3-parametercollisionoperator,allowingforseparaterelaxationofthe5shearmodes,1bulkmode,and9kineticmodes.Thepost-collision 28


2c4s;(2{8)wherethesoundspeedcs=x=p 2{8 shouldberetained,sinceitmaintainsGalileaninvarianceandpreventsanarticialcross-streamdrift[ 21 ].Empirically,wehavefoundthatthedomainofvalidityofStokesowisconsiderablylargerwiththenon-linearterminplace.Thenon-equilibriummomentumuxneq=Pinneqicicirelaxesduetocollisionsatthelatticenodes, 3(1+v)(neq:1)1;(2{9)whereneq=eq,eq=c2s+uuistheequilibriummomentumuxand 2;v=2c2s 2:(2{10)Thefactorof1=2correctsfornumericaldiusion,sothatviscousmomentumdiusesattheexpectedspeedforthegivenviscosity.Inthepresenceofanexternallyimposedforcedensityf,forexampleapressuregradientoragravitationaleld,thetimeevolutionofthelattice-Boltzmannmodelincludesanadditionalcontributionfi(r;t), 2c4st:(2{12) 29


21 ] 2ft:(2{13)Themacrodynamicalbehaviorarisingfromthelattice-Boltzmannequationcanbefoundfromamulti-scaleanalysis,[ 30 ]usinganexpansionparameter,denedastheratioofthelatticespacingtoacharacteristicmacroscopiclength;thehydrodynamiclimitcorrespondsto1.Itcanbeshownthatthelattice-BoltzmannequationreproducestheNavier-Stokesequationswithcorrectionsthatareoftheorderu2and2.Thus,atsucientlylowMachnumbers,themethodissecond-orderaccurateinspacewithrelativeerrorsproportionaltox2.Itisalsosecond-orderaccurateintimeiftheviscositiesaredenedaccordingtoEq. 2{10 3 4 ]incorporatesarandomcomponentintothemomentumuxduringthecollisionprocess(c.f.Eq. 2{9 ): 3(1+v)(neq:1)1+ff= whereRvisarandomvariablewithzeromeanandunitvariance,andRisasymmetricmatrixofrandomvariablesofzeromean.Theo-diagonalelementsofRhaveavarianceof1,whilethediagonalelementshaveavarianceof2.Inthisparticularimplementation,6randomnumbersarerequiredtogeneratethecomponentsofthesymmetricmatrixR,togetherwiththeconstraintthat 30


20 ]suggeststhat 4 21 ].Thecorrectcoecients,determinedbyrelatingthedecayoflongwavelengthstressuctuationstotheviscosityofthelattice-Boltzmannuidare[ 4 21 ] 2{15 ),butfor=v=1,anexactcorrespondencewiththeuctuation-dissipationrelationforcontinuoussystemsisobtained.Howeverthediscreteversion(Eq. 2{16 )canalsobeappliedtosystemswheretheviscouseigenvaluesarenotequalto1[ 4 ].Thismaybeadvantageousinproblemswhereatime-scaleseparationisrequiredbetweenhydrodynamicandkineticmodes. 2.4.1Finite-SizeSphericalParticlesInourimplementationoftheLattice-Boltzmannmodel,thesolidparticlesaredenedbyasurfacewhichcutssomeofthelinksbetweenlatticenodes(Fig. 2-1 ).Fluidparticlespropagatingalongthelinksinteractwiththesolidsurfaceatboundarynodesplacedhalfwayalongthelinks.Thisprovidesadiscreterepresentationoftheparticlesurfacewhichbecomesmorepreciseastheparticledimensiontogridsizeincreases.Adetailedexplanationofthetreatmentofsolid-uidboundariescanbefoundinreferences[ 2 4 16 ]. 31


Locationoftheboundarynodesforacircularobjectofradius2.5latticespacings.Thevelocitiesalonglinkscuttingtheboundarysurfaceareindicatedbyarrows.Thelocationsoftheboundarynodesareshownbysolidsquaresandthelatticenodesbysolidcircles. Amovingboundarycondition,withoutanyinterioruidisimplementedasfollows.Wetakethesetofuidnodesjustoutsidetheparticlesurface,andforeachnodeallthevelocitiescbsuchthatr+cbtliesinsidetheparticlesurface.AnexampleofasetofboundarynodesisshowninFig. 2-1 .Eachofthecorrespondingpopulationdensitiesisthenupdatedaccordingtoasimplerulewhichtakesintoaccountthemotionoftheparticlesurface.[ 4 ] 2cbtisthelocationoftheboundarynode.Weusethemeandensity0inEq. 2{17 insteadofthelocaldensity(r;t)sinceitsimpliestheupdateprocedure.Thedierencesbetweenthetwomethodsaresmall,ofthesameorder(u3) 32


2{17 2t)=x3 16 ]. 2{20 byneglectingtheinertialtermandintegratingtherandomcomponentofthehydrodynamicforceoverthedurationofthetimestep[ 31 ].InthisthesisevidencewillbepresentedtosupportourcontentionthatitismoreecienttointegrateEq. 2{20 directly.Thelargetime-scale 33


21 ].Thecouplingbetweenthebeadsandtheuidisderivedfromafrictionalforcebasedonthedierenceinvelocitybetweenthebead,U,andthesurroundinguid,u(r;t),[ 11 ] 11 ].Sincetheuidsatisesitsownuctuation-dissipationrelation,Frhasalocalcovariancematrix 3.2.2 .Thetimescaleforhydrodynamicinteractionstodevelop,R2g=,mustbemuchshorterthanthediusiontimeforthepolymer,R2g=D.Inotherwords,viscousmomentummustdiuseconsiderablyfasterthan 34


32 ]andtheeectofvariationsinSchmidtnumberisusuallysmall.TheEulermethodwasusedtointegratetheequationsofmotionofthemonomerswithatimestepthatwastypically1/100ofthemonomerdiusiontime.Whennecessary,thethermodynamicforceswereintegratedwithasmallertimestepthanthehydrodynamicforcestomaintainstability.Ahard-spheremoleculardynamicscodewasusedateachsteptodetectcollisionsbetweenpairsofmonomers.Sincethemonomersmovecontinuouslyoverthegrid,whilethevelocityeldisevaluatedatadiscretesetofpoints,weemployatrilinearinterpolationtoevaluateu(X;t),usingtheuidvelocitiesatthenodesofthecubesurroundingthemonomer[ 11 ].Aftercalculatingtheweights,wi,foreachofthenodes,themassdensity,,andthemomentumdensityuareinterpolatedtodeterminethevelocityeldatthebeadlocationX, 35


(2{25) (2{26) whereisahighfrequencyresistance,andsubscripts\imp"and\exp"standforimplicitandexplicitupdateschemesrespectively.Thusfortheimplicitupdatetheequationbecomes mt)=U(t)+F0 mt)=U(t)+t mU(t)t WeeasilyrecognizetheterminbracketsasFexp U=(+ m)1Fexp 2 16 ]. 36


4 ].Neverthelessinthissameworkitwasshownthatare-scalingofthevelocityautocorrelationfunctionshowsagoodagrementwiththevelocitydecayofaparticlesetinmotionwithanimpulsivekick,accordingtoOnsager's"LinearResponseTheory"[ 15 ].Thecurrentstudyextendstheresultsofthesepapers[ 4 15 ]andthatof[ 2 ]withtheimprovementsinthealgorithmsuchasbetterhandlingoftheboundarynodes,additionofastationarymodetothe18velocitymodelandanadditionalparameter(v)forthedecayofthenormalstressmodesinaneorttostudytheself-diusionincolloidalsuspensions.Thissectionconsistsofthetestswehaveperformedtoverifyoursimulationmethod.Thelattice-Boltzmannmethodusedinthisstudyisanisothermalmodelandfortheuctuatingschemeandthetemperatureischaracterizedbythestressuctuations.Wehavecalculatedtheseuctuationsoveranode,overaplane,andoverthewholeperiodicbox.Theaveragestressuctuationsarelocalinspaceandtimewithavariancegivenby 2{31 upontheuseofdierentstressrelaxationparameters(andv)andalsointhepresenceofsolidparticlesofconsiderablevolumefraction.Inthecurrentmodel,theonlyinteractionbetweenastationaryplanarwallandauidismomentumexchangebetweenthetwoarisingfromtheuctuationsintheuid.Aphysicallysimilarproblemisanoscillatingplanarwallinastationaryuid, 37


C ).AdiscretetheorycanalsobeconstructedbasedsimplyontheinformationfromlatticeBoltzmanntheorytocalculatetheinstantaneouscorrelations(Appendix D ).Figure 2-2 illustratestheagreementofthesimulationresultsandthetheoreticalexpression.Thetheoreticalexpressionisslightlyinerrorattherstfewtimestepsandthiserrorturnsintoaconstantosetaftertherstfewsteps,implyingasimilarrstandsecondderivativeforthetwocurveswhichdeterminethetransportcoecients.Asimilarandpracticalproblemforthesimulationofcolloidalsuspensionsistheinteractionofaninnitelymassive,immobile,spherewiththeuctuatinguid.Herethedierentfrictionkernelisdierentduetothedierencesingeometry.Usingthefrictionkernel[ 33 ]themomentumexchangecorrelationsbetweentheuidandthespherecanbeanalyticallycalculated(Appendix B ).TheagreementbetweenthetheoreticalexpressionandsimulationvaluesaregiveninFigure 2-3 .Argumentssimilartothecaseoftheplanarwallforthetransportcoecientsalsoapplyhere.Thesetwoexampleswithdierentgeometriesshowthatthemethodiscapableofproducingthecorrecttransportcoecients.Thevelocityauto-correlationfunctionofaparticlerevealstransportpropertiessuchasthediusioncoecientandcharacteristictimescalesfortherelaxationoftheparticles.AccordingtothelinearresponsetheoryofOnsager,thevelocityrelaxationofasingleparticlewhichisgivenanimpulsivekickisequivalenttotherelaxationofthevelocity 38


Momentumchangecorrelationsonaplanewall.P(t)isthetotalmomentumexchangebetweentheuidandtheparticleatanyinstantt.Resultsfromsimulationsarecomparedwithatheoreticalexpressionderivedfromcontinuumtheory.Thetheoreticalexpressionishigherthanthesimulationdataatthersttime-step,erroraccumulatesforacoupleoftime-stepsandthenthedierencebetweenthesetwocurvesturnintoaconstantoset.Therstpointforthiscorrelationiscalculatedcorrectly(1%)fromadiscretetheorystartingwiththeuctuationsinpopulationdensities(section D ). uctuationsofaBrownianparticlegoingthroughmanyrandomkickssincethemeansofrelaxationisthedissipationoftheenergyinthesurroundinguidinbothcases.Figure 2-4 demonstratesthisforoursimulationtechniqueandisyetanothermeanofvericationcode.AsillustratedinFigures 2-4 and 2-5 ,asingleBrownianparticleneverreachesitsequipartitionvelocityuctuationsduetotheartifactsofthediscreteLangevinlikeequationutilizedinoursimulations.Uponintegrationofthisdiscreteequation,usingaconstantfrictionterm,tocalculatethevelocityuctuationsandtheautocorrelation 39


Momentumexchangecorrelationsonaninnitelymassivesphere.P(t)isthetotalmomentumexchangebetweentheuidandtheparticleatanyinstantt.Resultsfromsimulationsarecomparedwithatheoreticalexpressionderivedfromcontinuumtheory. itisfoundthatthethevelocityuctuationsdierfromtheequipartitionvaluebyafactorof2 2t.Thesigncorrespondstotheexplicitandimplicitintegrationschemes.ThustheratiooftheBrownianrelaxationtimesteptothesimulationstepplayanimportantroleforthevelocityuctuations.Figure 2-5 showscalculationsofthevelocityautocorrelationfunctionfordierentvaluesofandillustratestheeectoftheartifact.Yetoursimplifyingassumptionofaconstantfrictiontermdoesnotapplytothesimulationsandthetheoryisonlyqualitativelycorrect.Nevertheless,thediusioncoecientcalculatedbyintegratingthevelocityautocorrelationfunctionaccordingtotheGreen-KuborelationshipisfreeofthisartifactandgivesthecorrectStokes-Einsteinvalue(equation 2{32 )asveriedbysimulations.ThesecalculationsaresummarizedinAppendix 40


VericationofOnsager'slinearresponsetheory.Thedecayofthevelocitycorrelationsiscomparedwiththevelocitydecayofaparticlewithaninitialgivenvelocity.Thetwocurvesleadtothesamediusioncoecientuponintegration. A .Thesepreliminarytestsshowthattheuctuatinglattice-BoltzmannschemeiscapableofproducingthecorrecttransportcoecientsforsphericalBrownianparticlesof :(2{32)Thisconcludesourtestforcolloidalnitesizedparticles.FurthertestsonthehydrodynamicinteractionsandscalingbehaviorofpointparticleswillbediscussedinChapter 3 41


Decayofthevelocityautocorrelationfunctionforasingleparticlewithmassfactorof1&10andaradiusof5.4,normalizedbytheequipartitionvalueU2=kT=mB.Thetwodierentcurvescorrespondtodierentintegrationschemes.Thesetwocurveseventuallyleadtothesamediusioncoecient. 42


2 ,tosimulatethedynamicsofpolymersolutionsinconnedgeometrieshasbeenimplementedandtested.Themethodcombinesauctuatinglattice-Boltzmannmodelofthesolvent[ 34 ]withapoint-particlemodelofthepolymerchains.Africtiontermcouplesthemonomerstotheuid[ 11 ],providingboththehydrodynamicinteractionsbetweenthemonomersandthecorrelatedrandomforces.Thecoupledequationsforparticlesanduidaresolvedonaninertialtimescale,whichprovestobesurprisinglysimpleandecient,avoidingthecostlylinearalgebraassociatedwithBrowniandynamics.Complexconnedgeometriescanberepresentedbyastraightforwardmappingoftheboundarysurfacesontoaregularthree-dimensionalgrid.ThehydrodynamicinteractionsbetweenmonomersareshowntocomparewellwithsolutionsoftheStokesequationsdowntodistancesoftheorderofthegridspacing.Numericalresultsarepresentedfortheradiusofgyration,end-to-enddistance,anddiusioncoecientofanisolatedpolymerchain,rangingfrom16to1024monomersinlength.Thesimulationsareinexcellentagreementwithrenormalizationgroupcalculationsforanexcludedvolumechain.Weshowthathydrodynamicinteractionsinlargepolymerscanbesystematicallycoarse-grainedtosubstantiallyreducethecomputationalcostofthesimulation.Finally,weexaminetheeectsofconnementandowonthepolymerdistributionanddiusionconstantinanarrowchannel.OurresultssupportthequalitativeconclusionsofrecentBrowniandynamicssimulationsofconnedpolymers[ 35 36 ]. 43


36 37 ]reducesthistoN3:8,butwithpossiblelossofpositivedeniteness[ 36 ].Still,themaximumchainlengththatcanbestudiedbyBrowniandynamicsisoftheorderof100monomers[ 35 36 38 { 40 ],whichisnotalwayssucienttodetermineaccuraterheologicalpropertiesofexiblepolymers[ 41 ].Inasemi-dilutesolutionofNcchains,thecomputationalcostofBrowniandynamicsincreasesfurther,inproportiontoN2corN3cdependingonthemethodoffactorization.Inadditiontotheunfortunatescalingofthecomputationalcost,Browniandynamicsisnoteasilyextendedtoincludeexternalboundariesandowelds.No-slipboundariesrequireacomplicatedGreenfunctionevenforchannelows[ 36 ],whilemoreinvolvedgeometriescanonlybehandledbynite-elementorboundary-elementmethods.Browniandynamicshasnostraightforwardwayofincludinganyexternaloweldbeyondasimplelinearshearow.Morecomplexowsrequireacoupledcalculationoftheevolutionoftheuidvelocityeldwiththedynamicsoftheparticlechains.EventhegreatadvantageofBrowniandynamics,whichisthecapabilitytosolvedirectlyforthedynamicsonthediusivetimescale,isnegatedbytheverysmalltimestepthatisnecessarytointegratethestochasticequations,typically105tZformoderatelengthchains(N100).Wewillseethatcomparabletimestepscanbeusedininertialsimulations. 44


2 ,forindividualmonomersaswellasapplicationstothediusionofpolymerchainsinperiodicandconnedgeometries.WeshowthattheOseenlimitofthedissipativeanductuatinghydrodynamicinteractionsbetweentwomonomersisobtainedforseparationslargerthanx,thegridspacingofthelattice-Boltzmannmodel.Wealsoshowthattheradiusofgyration,end-to-enddistance,anddiusioncoecientofaexibleexcludedvolumechainagreequantitativelywithrenormalizationgroupcalculations[ 42 ]forchainsupto1024monomers.Acoarse-grainingofthehydrodynamicinteractionshasbeeninvestigated,whichsuggeststhepossibilityofaconsiderablespeedupinthecomputationforlargechains.Despitetheextrainertialtimescale,oursimulationsusetimestepsthatarecomparabletoBrowniandynamics,andthereforeachieveamorefavorablescalingofthecomputationaltimewithchainlength.Webrieyinvestigatetheeectsofconnementandowonthepolymerdistributioninatwo-dimensionalchannel. 2{16 )and0isthemonomerfriction.However,asAhlrichsandDunweg[ 11 ]originallyindicated,thediusioncoecientofanisolatedpointparticleisnot 45


Relationbetweentheeectivefrictioncoecientandtheinputfriction0(Eq. 2{21 ).Theeectivefrictioncoecienthasbeenobtainedfromthedragforceonasinglemonomer(circles)andthediusionofthemonomerinathermallyuctuatinguid(squares).Theslopeofthelineisunity,indicatingaconstantosetbetween1and10.TheunitcellsizeL=20x. preciselyproportionaltotheinverseofthefrictioncoecient,butisosetbyaconstantvaluethatisindependentoftheinputfriction0, x:(3{1)WethereforeattemptedtoverifytheStokes-Einsteinrelationforasinglemonomerbycalculatingthevelocityofamonomerunderanexternalbodyforceaswellasthediusivityofthemonomerduetothermaluctuationsintheuid.Fig. 3-1 summarizesthemostimportantnding,namelythatthefriction-coecientdeterminedfromthemean-squaredisplacementofthemonomeratnitetemperatures,1=Dm=T,matchesthefrictioncoecientdeterminedfromthedragforce,Fd=U.Figure 3-1 alsoshows 46


2{21 ).Inessence,thereisarenormalizationofthefrictioncoecientbytheoweldinducedbythemovingparticle.ApointforceFappliedattheoriginofacontinuumuidcreatesavelocityeld, x;(3{3)wheregisanumericalconstant.Asteady-stateforcebalancebetweentheimposedforceFandtheparticlevelocity, 3{1 .Furthertestsverifythattheparametergisindependentofsolventviscosity,conrmingtherelationgiveninEq. 3{1 .However,gdoesdependweaklyonsystemsizeasshowninFig. 3-2 ,butnotontheinputfriction,0.Theoset,g,iseasilycalculatedforanycellgeometry,usingameasurementoftheuidvelocityattheoriginofapointforcetogetherwithEq. 3{3 .Whentheeectivesizeofthemonomerissmallcomparedtothegridspacing(weakfriction),x=1andtheeectofthelatticerenormalizationissmall.Inthiscasethediusivityiscontrolledbytheinputparticleradius,Dm=T!1=6a0.Theoppositelimit,wheretheparticleislargecomparedtothegridspacing,violatestheunderlyingassumptionofthepoint-particlemodel.However,evenherethediusivityisnitebutcontrolledbythegridspacing,Dm=T!g=x.Inthesenumericalsimulationsthelatticecontributiontothediusivityisoftheorderof20%.Theuctuation-dissipationtheoremforadiscretemodeldoesnotnecessarilyimplyequipartitionofenergy,asisthecaseforacontinuoussystemdescribedbytheclassicalLangevinequation.Thedatain 47


Osetbetween1and10asafunctionofsystemvolume.Thevariationinthedimensionlessmobilityg(Eq. 3{1 )isshownforacubiccelloflength,L.Resultsareshownforinputhydrodynamicradiia0=0:32(circles)anda0=0:15(squares).Themonomerfriction0=6a0. Fig. 3-3 showsthatthemeankineticenergyofasinglemonomerinathermallyuctuatinguidvariesbetween1%and10%fromequipartition.TheresultscanbeunderstoodbyanalyzingtheLangevinequationforasingleparticlewithconstantfriction, 48


Meankineticenergyofasinglemonomerasafunctionofthedimensionlessfrictionalcoupling,t=m.Thecalculatedkineticenergyforimplicit(circles)andexplicit(squares)updatesiscomparedwiththemodeldescribedinAppendix A (lines). or 3{6 correspondstoarst-orderexplicitintegrationofEq. 3{5 ,whileEq. 3{7 correspondstoanimplicitintegration.InAppendix A ,weshowthatthischoiceofrandomforce,deriveddirectlyfromthecontinuouscase,stillsatisestheStokes-EinsteinrelationD=T=,regardlessoftheintegrationmethod,butclassicalequipartitionofenergyismodied, 3{8 ,whileimplicitintegrationleadstotheplussign.Figure 3-3 showsthatthissimplemodelaccounts 49


2 providesanaccuratesolutionfortheow-eldaroundapointforce,ascanbeseeninFig. 3-4 whereasectionoftheuidvelocityeldisshownforaperiodiccelloflengthL=100x.Thisexampleistypicalofthemethod,whichquantitativelydescribestheoweldondistanceslargerthanasinglegridspacing,x.ThereisasignicantdierenceatlargedistancesbetweentheOseeneldforaninnitesystem(dashedline)andStokesowinaperiodiccell(solidline),whichwascalculatedfromtheFourierrepresentation,[ 43 ] 50


Velocityduetoapointforce.Thevelocitycomponentparalleltothedirectionoftheforce,vx(r;0;0),isplottedasafunctionofdistancefromtheforce.ResultsareshownfortheOseenoweld(dashed),theStokesoweldinaperiodiccelloflength100x(solidline)andlattice-Boltzmannsimulations(circles)inthesamesizecell. ThelinearityofStokesowallowsustosimplifytheanalysisbykeepingtheparticlesxedinplaceduringthecourseofthesimulation,eventhoughtheparticlevelocitiesarenon-zero.Thesingle-particlemobilitywasfoundtobeunaectedbythepresenceofaneighboringparticle,aswouldbeexpectedforpointforces.Thepairmobility,whichcanbedecomposedintocomponentsparallel(k)andperpendicular(?)totheseparationvector,isshowninFigs. 3-5 and 3-6 .TheresultsfordissipativeanductuatingcalculationsofthemobilityagreewithoneanotherandwithanexactStokesowcalculationinthesameperiodicgeometry,Eq. 3{9 ,downtoseparationsoftheorderofx.Thediscrepancieslargelyreecterrorsininterpolatingtheoweldtoo-lattice 51


Comparisonoftheparallelcomponentofthetwoparticlemobilitytensor,k12.Resultsobtainedbydissipative(Eq. 3{10 ,circles)anductuating(Eq. 3{11 ,squares)simulationsarecomparedwithanindependentsolutionoftheperiodicStokesequations(Eq. 3{9 ,triangles). locations;theresultsshowninFig. 3-4 wereevaluatedatthelatticenodesandaremuchclosertotheStokesowresultatsmalldistances.Improvedinterpolationshouldleadtomoreaccuratevaluesofthemobilityofcloselyspacedmonomers. 11 ]andb=0:5x.Themonomerfrictionandexcludedvolumewere 52


Comparisonoftheperpendicularcomponentofthetwoparticlemobilitytensor,?12.Resultsobtainedbydissipative(Eq. 3{10 ,circles)anductuating(Eq. 3{11 ,squares)simulationsarecomparedwithanindependentsolutionoftheperiodicStokesequations(Eq. 3{9 ,triangles). alsoscaledinproportiontob.Thecomputationaltimescalesasb6,takingintoaccountboththechangeinnumberofgridpoints(/b3)andthechangeinZimmtime(/b3).Thesimulationswererunforatleast20Zimmtimes,whichwassucienttoreducethestatisticalerrorsinthediusioncoecienttoafewpercent.Forthesystemsstudiedhere,theZimmtimerangesfrom103tto106t.Giventhescalingofthecomputationtime,itisclearlyadvantageoustominimizethevalueofb=x.ResultsfortheradiusofgyrationRgandthemeanend-to-enddistanceReareshowninFig. 3-7 .ThestatisticalpropertiesofthepolymerchainsareindependentofbandT,andscaleintheexpectedwayforanexcludedvolumechain.ThenumericalvaluesoftheradiusofgyrationtthescalingrelationRg=0:4bN,with=0:59.Thisiscloseinabsolutesizetotherenormalizationgroupresultforaself-avoidingrandomwalk[ 42 ]. 53


End-to-enddistanceandradiusofgyrationinaperiodicunitcell.ThescalingofRe(upperpoints)andRg(lowerpoints)isshownfordierentmonomerseparations,b,andtemperatures,T.Thedashedlinesarethetheoreticalscalingsforanexcludedvolumechain,withexponent=0:59. Thelattice-Boltzmannmethodsimulatesanitevolumeratherthananunboundeddomain,whichhaslittleeectonthedistributionofpolymerconformationsbutsignicantlyreducesthediusivitybecauseofhydrodynamicinteractionsbetweenperiodicimages.Theboxsizewaschosentobeaxedmultipleoftheexpectedradiusofgyrationtoobtainaconsistentsetofdiusioncoecients,asshowninFig. 3-8 .Therawdiusioncoecients,shownfortwodierentboxlengths,L3:5RgandL7Rg,scaleroughlyasN,butthediusioncoecientsincreasewiththeratioL=Rg.Thenite-sizeeectscanbecorrectedbyassumingthatthepolymerbehaveshydrodynamicallyasarigidsphereofradius,Rh,where 54


Centerofmassdiusioncoecientasafunctionofchainlengthindierentsizecells.Thediusioncoecientsarenormalizedbythemonomerdiusioncoecient,Dm=T=,andareshownforperiodicsystemswithlengthsL3:5RgandL7Rg.Inallcasesthemeanseparationbetweenthemonomersb=xandthetemperatureT=0:1.Resultsareshownfortherawdiusioncoecient(opensymbols),aswellasthecorrectedvalues(lledsymbols). Thenwecanusethemobilityofaperiodicarrayofspheres[ 43 ]tocorrectfortheeectsoftheimagechains,[ 44 ] 3{12 ),aself-consistentestimateofthediusioncoecientinaninnitesystemcanbeobtained.Figure 3-8 showsthatthecorrectedresultsobtainedwithL34RgaresimilartothosewithL67Rg,andthatthecorrectedresultsfallpreciselyontheexpectedscalingofdiusivitywithchainlength,D/N.Theresultsagreewellwiththerenormalizationgrouppredictionforaself-avoidingrandomwalk[ 42 ],D=0:2T=bN.Themonomerfrictioninthese 55


Centerofmassdiusioncoecient,withnite-sizecorrections,fordierentvaluesofbandT.ThelengthoftheperiodiccellL7Rgineachcase.Thestatisticaluncertaintiesinthediusioncoecientaresmallerthanthesizeofthesymbols,exceptforthelongestchain,N=1024,wheretheyareshownexplicitly. simulationswasbasedonahydrodynamicradiusa=b=4,sotherenormalizationgrouppredictionfortheratioofpolymertomonomerdiusioncoecientsisD=Dm=0:9N.ThelargestsimulationsinFig. 3-8 extendupto256monomerunitsandrequireabout5hoursofcomputationperZimmtime.Howeverthecomputationalcostcanbedrasticallyreducedbyeitherreducingtheseparationbetweenneighboringmonomersalongthechain,b,orlessdramatically,byincreasingthetemperature.Reducingblowerstheaccuracyofthehydrodynamicinteractionsbetweennearbymonomers,whileincreasingthetemperaturealterstheSchmidtnumberandmaypreventtheStokesoweldfromfullydevelopingbeforethemonomerpositionschange.DiusioncoecientsfordierentvaluesofbandTareshowninFig. 3-9 .Inmostsimulations,weusedb=xasinAhlrichsandDunweg[ 11 ],butforlargerchainswefoundthatsmallervaluesofbcouldbe 56


7 ]andtheadditionalhydrodynamicinteractions 57


Centerofmassdiusioncoecientsparallel(Dk)andperpendicular(D?)tothechannelwallsasafunctionofchannelwidthH;theradiusofgyrationofthepolymer(N=128,b=1)isapproximately7x.Theconneddiusioncoecients,normalizedbythemonomerdiusivity,areshownforasquarecell,L,withvariousratiosofLtoH. betweenpairsofmonomersinthevicinityofthewall[ 45 ].Thesimulationsintheconnedgeometrywererunforcomparabletimestothebulksimulations,whichtypicallycorrespondsto40-120channeldiusiontimes,tH=H2=4D.However,fortheweakestconnement,H=Rg>10,theruntimewasoftheorderof10tH.Thediusioncoecientforaconnedpolymer(N=128,b=1andT=1)isshowninFig. 3-10 forarangeofchannelwidthsRg

3-11 forarangeofconnementratios,Rg=H,withaconstantaspectratioL=H=4.ThediusioncoecientDk=Dofaweaklyconnedpolymer,Rg=H1,isconsistentwiththereductioninmobilityofaconnedsphere[ 7 ]ofradiusRh.Intighterconnement,Rg=H>1,thereisatransitiontoapower-lawdecay,whichisconsistentwiththepredictionsofscalingtheory[ 46 ],Dk=D/(Rg=H)2=3.OntheotherhandBrowniandynamicssimulationsinasquarechannel[ 35 ]suggestanexponentcloserto1=2.Thisresultissupportedbynumericalmean-eldcalculations[ 47 ],whichindicatethattheasymptoticscalingisnotreachedintwo-dimensionalconnementuntilextremelylongchainlengths,inexcessof106segments.Howevertheconnementinaslitislessseverethaninatube,andamean-eldcalculationwouldbeusefultoconrmthatscalingtheoryappliestomuchshorterchainsinone-dimensionalconnement. 59


Centerofmassdiusioncoecientforaconnedpolymerrelativetotheunconneddiusivity,Dk=D.ThesolidlineindicatestheFaxencorrection[ 7 ]forweakconnement,assumingthepolymerliesinthecenterofthechannel.Thedashedlineisthetheoretical(Rg=H)2=3scaling[ 46 ].Theaspectratioofthecell,L=H,wassetto4ineachcase. Surprisingly,thetimestepiscomparabletotypicalBrowniandynamicssimulations.Thisfact,togetherwiththemorefavorablescalingofthecomputationaltimewithrespecttothesizeofthechain,hasenabledustoperformdynamicalsimulationsoflongerchainsthanhaspreviouslybeenpossible.WeobtainedtheexpectedscalingforthediusionconstantinanunconnedsystemforchainlengthsuptoN=1024monomers.Theresultssuggestthatacoarsegrainingofthehydrodynamicinteractionsispossibleforverylongchains,keepingtheradiusofgyrationataxednumberofgridpointsoftheunderlyinguid.Inthiscase,thecomputationalcostofcalculatinghydrodynamicinteractionsisindependentofchainlength.Themethodcanbereadilyextendedtocomplexconnementgeometries,withorwithoutanimposedoweld.OurresultsforconnedchannelsarequalitativelyconsistentwithBrowniandynamicssimulationsina 60


35 36 ].Althoughallthesetestshavebeenforasinglechain,multi-chainsimulationsarestraightforwardandrequireminimaladditionalcomputationifthesystemvolumeremainsxed. 61


48 { 50 ].Inapressure-drivenow,theshearratestretchesandalignsthepolymeralongtheowdirection,reducingitscongurationalentropy.Thermodynamicargumentsthereforepredictthatthepolymerwillmigratetothecenterlinewherethelocalshearrateisminimized[ 51 52 ],creatinganapparentslipattheboundaries[ 53 ].Thethicknessofthedepletionlayerisobservedtoincreaseinapressure-drivenow[ 54 55 ],whichalsoimpliesanetmigrationofthepolymertowardsthechannelcenter.Bycontrast,ithasusuallybeenassumedthatthedepletionlayerthinsinauniformshearowduetoreducedsterichindrancebetweenastretchedpolymerandtheboundary[ 53 ].Recenttheoryandsimulationshaveunderlinedtheimportanceofhydrodynamicinteractionsindescribingpolymermigration.Kinetictheories[ 56 57 ]attributelateralmigrationtotheasymmetryofthehydrodynamicinteractionsbetweenastretched 62


17 36 ],includingfullhydrodynamicinteractions,showedthatthepolymermigratestowardsthechannelcenterlineinpressuredrivenows.Ontheotherhand,simulationsneglectinghydrodynamicinteractionsbetweenpolymersegments[ 39 40 ]ndmigrationtowardsthewalls.Therearequalitativedierencesbetweenthepredictionsofkineticandthermodynamictheoriesofmigration.Inparticular,inauniformshearow,nomigrationisexpectedonthermodynamicgroundssincethepolymerextensionisindependentofposition.However,thekinetictheoriespredictmigrationtowardsthechannelcenterinauniformshearowaswellasPoiseuilleow.Todistinguishbetweenthesetheories,wehaveusednumericalsimulations,withhydrodynamicinteractions,tocalculatethedirectionandextentofthelateralmigrationinpressure-drivenanduniformshearows.Wediscoveredthatmigrationcanoccurinbothdirections,dependinguponthedegreeofconnement;inwidechannelsthepolymersmovetowardsthecenter,whileinnarrowchannelstheymovetowardsthewalls.ThepolymerdistributionisquantitativelydierentinshearandPoiseuilleows,butthedirectionofmigrationisthesameinbothcasesforthesamedegreeofconnement.Thissuggeststhattheprimarydrivingforceformigrationishydrodynamicliftfromtheconningwalls,drivenbytheshearrateratherthanitsgradient[ 56 57 ]. 4 21 ],presentedinChapter 2 .Thepolymermotioniscoupledtothesurroundinguidusingafrictionalforce,basedonthedierencebetweenthebeadvelocityandthelocaluidvelocity[ 11 29 ].Theforceexertedbythebeadsisredistributedtothesurroundinguidnodessothatthetotalmomentumofthecoupledparticle-uidsystemisconserved.ExtensivenumericaltestshaveshownthatweobtainOseenlevelhydrodynamicinteractionsbetweenbeads 63


17 ].Inaddition,thereisasoftrepulsionbetweenpairsofbeadsandbetweenbeadsandconningwalls.Theexcludedvolumespherehasaradiusofaboutb=2,wherebisthedistancebetweenbeads.Theuctuatinglattice-BoltzmannapproachhascomputationaladvantagesoverBrowniandynamics[ 58 ],whichhasbeentheprimarynumericalmethodforsimulatingpolymersolutions.First,thecomputationalcostscaleslinearlywiththenumberofbeads,allowingforsimulationsofchainsinexcessofathousandbeads[ 17 ].Secondcomplexgeometries,suchasporousmediaormicrodevices,canbeintroducedbymappingthesolidsurfacesontothelattice,withoutintroducinganewGreenfunctionforeachspecicgeometry.Westudyasinglepolymerchain,connedbetweentwoatplatesseparatedbyagapH.Periodicboundariesareusedintheothertwodirectionsandtheperiodiclengthinthesedirectionsissucientlylarge(L=2H)thatthehydrodynamicinteractionbetweendistantimagesisnegligible[ 17 ].ThepolymerchainsarediscretizedintoNbeadswithmostsimulationsreportedinthispaperusingN=16.Longerchains(N=3264)areusedtoverifythattheequilibriumradiusofgyration(Rg)andthefreesolutioncenterofmassdiusioncoecient(D)aretheonlyimportantparametersandthatthelevelofdiscretizationofthechainissucient.Apressuredrivenowwascreatedbyapplyingauniformforcetotheuid,whileauniformshearowwasgeneratedbyrelativemotionoftheboundingwalls.ThePecletnumberisbasedontheshearrateandtheequilibriumradiusofgyration,Pe=_R2g=D,inbothcases.Inauniformshearow,theshearrateisconstantacrossthechannelandgivenby_=V=H,whereVisthevelocitydierencebetweenthewalls.InPoiseuilleow,themeanshearrate,_=2Vmax=H,isused,whereVmaxisthecenterlinevelocity.Thedegreeofconnementischaracterizedbytheratioofthegapwidthtotheequilibriumradiusofgyrationofthepolymerchain,H=Rg.Forthesmallestconnementratiousedinthiswork,H=Rg=3,thepolymercoilremainsfreetorotateinthechannel.Thesimulationsarerunforasucientlylongtimethatapolymer, 64


4-1 a,conrmthatpolymerchainsmigratetowardsthechannelcenterlineinapressuredrivenow.Themechanismforthismigrationisnotentirelyagreedupon.Thethermodynamictheoryarguesthatthesheargradientdrivesthepolymertowardslowshearregionstoattainmaximumpossiblecongurationalentropy[ 51 52 ].Wehavetestedthisassumptionbyrepeatingthesimulationsinauniformshearow.Thesimulationsshowthatthepolymersstillmigratetowardsthecenterline(Fig. 4-2 a),indicatingthatasheargradientisnotrequiredformigration.However,duetosymmetry,nomigrationwouldoccurinanunboundedshearow;thereforeshearcannotbyitselfexplainthemigrationeither.Hydrodynamicargumentssuggestthatthetensionintheshearedpolymergeneratesaowawayfromthewall,duetotheasymmetryofthehydrodynamicinteractionsbetweenapairofbeadsandaplanarboundary[ 55 57 ].Akinetictheory[ 57 ]ofadumbbellnearaplanarwallpredictsmigrationtowardsthechannelcenterinbothuniformshearandPoiseuilleowswhenhydrodynamicinteractionsareincluded,inagreementwithsimulations.Adumbbellismuchlessexiblethanthemulti-beadchainsusedinoursimulations.Moreover,thehydrodynamiccenterofadumbbellisalsoitscenterofmass,soanexternaltorqueleadstoapurerotationratherthantumbling.Despitetheselimitations,thekinetictheoryqualitativelyexplainsthekeyfeaturesobservedinsimulationsofowinwidechannels.Afeatureofthepolymerdistributioninwidechannelsisasymmetricdoublepeak(Fig. 4-1 a)observedathighPecletnumbers[ 17 36 ].Thisisaconsequenceoftheinterplayofhydrodynamicforceswithaspatiallyvaryingdiusivity[ 57 ].Ashearedpolymerundergoessignicantextensionandalignmentalongtheowaxis(Fig. 4-3 a),hencethediusivityalongtheconnementdirectionisreducedwithincreasingshearrate.Thecompetitionbetweenthespatiallyvaryingdiusivity 65


CenterofmassdistributionforpressuredrivenowasafunctionofPecletnumberandchannelwidth.ResultsareforN=16.ThedistributionsarenormalizedsuchthatRH=Rg0P(y=Rg)d(y=Rg)=1.Duetothesymmetryofthesystem,thedistributionisplottedupto0:5H.Theboundingwallisattheorigin. andthehydrodynamicallyinducedmigrationleadstothedoublepeakinthepolymerdistributionfunctionP(y=Rg)athighPecletnumbers.Kinetictheory[ 57 ]wouldalsopredictadoublepeakinthedistributionfunctionathighPecletnumbersifaspatiallyvaryingdiusivitywasusedinthemodel.Fig. 4-1 ashowsthatthepeakpositioninthechannelisindependentofPecletnumber.Empirically,weobservethatthepositioncoincideswiththepolymerreachingabout90%ofitsmaximumextensionintheow,asseenbycomparingFigs. 4-1 aand 4-3 a.Themobilityalongtheconnementdirection,whichisalogarithmicfunctionoftheaspectratio,changeslittlewiththemarginaladditionalextensionthatoccursinthevicinityofthewalls.Consequently,thedrivingforceawayfromthecenterline,whichisproportionaltothegradientofthediusivity[ 57 ],diminishesveryquicklybeyond 66


CenterofmassdistributionforuniformshearowasafunctionofPecletnumberandchannelwidth.ResultsareshownforN=16.ThedistributionsarenormalizedsuchthatRH=Rg0P(y=Rg)d(y=Rg)=1. thispoint.Inauniformshearthedistributionhasatmostasinglemaximum,sincethediusivityisessentiallyconstant.Asthechannelwidthisreducedwithrespecttothesizeofthepolymer(H=Rg=5),themigrationtowardsthecenterdecreases(Fig. 4-1 b).Theextensiondata(Fig. 4-3 b)showsthatthepolymershapeissimilartothatinthewiderchannels,indicatingthatthereducedmigrationisaresultofthesmallerhydrodynamicforce.AtachannelwidthofH=Rg=4,migrationtowardsthecenterisonlyobservedforthehighestPecletnumber.AtlowerPecletnumbersthemigrationistowardsthewalls(Fig. 4-1 c).Instilltighterconnement(H=Rg=3),migrationistowardsthewallsatallPecletnumbers.Weobserveasimilaroutwardsmigrationinuniformshearow(Fig. 4-2 b),althoughtheextentofthemigrationislowerthaninPoiseuilleow.Migrationtowardsthechannelwallsisfurthersupportedbyresultsforthesecondmomentofthecenter-of-massdistribution[ 36 ](Fig.6).Nevertheless,theevidenceforoutwardmigrationinthisworkisbynomeansunequivocal;indeedtheauthorsconcludethatthepolymeralwaysmigratestowardstothecenterifhydrodynamicinteractionsareincluded.Polymerkinetictheory[ 57 ]accountssemi-quantitativelyformigrationinwidechannels,correctingearlierwork[ 56 ]basedonsimilarideas.Thekeymigrationmechanism,namelytheliftfromthewallinducedbyashearedpolymer,hasbeenvalidatedinprevious 67


Componentsoftheradiusofgyrationalongtheow,hx2i,andconnement,hy2i,axes.ResultsforchainsoflengthN=16areshownattwodierentPecletnumbers,normalizedbyRg=(hx2i+hy2i+hz2i)0:5. work[ 17 36 ]andconrmedinthepresentstudy.However,thehydrodynamicinteractionsinthechannelwereapproximatedbysuperposingtheoweldsgeneratedbyeachwall.Thisisvalidwhenthepolymerismuchclosertoonewallthantheother[ 7 59 ],butinanarrowchanneltheGreen'sfunctionforaslitmustbeused.Innarrowchannels,superpositionover-predictstheextentofthehydrodynamicinteractions,andthereforethedegreeofmigrationtowardsthecenter.Inreality,thehydrodynamicforcesthatcausemigrationtowardsthecenterarescreenedbythecloselyspacedboundaries.Wehavefoundthatthepolymerdiusioncoecientalongtheconnementaxisshowsacross-overtoRouse-likescalingintherange8>H=Rg>4,whichsupportsthescreeninghypothesis.Moreover,computersimulations[ 39 40 ]whichneglectintra-beadhydrodynamicinteractionsalsoshowmigrationtowardsthewalls,independentofthechannelsize. 68


Centerofmassdistributionfordierentlevelsofdiscretizationofthepolymerinatightlyconnedchannel,H=Rg=3.ThedistributionsarenormalizedsuchthatR10P(y=Rg)d(y=Rg)=1. Theshapeofthepolymermayalsoplayanimportantroleinthereversalofthedirectionofmigration.Theextensioncurves(Figs. 4-3 candd)suggestthatequilibriumconformationsbecomedistortedwithin1:5Rgofthewalls.Thereforehighlyconnedpolymersareanisotropicevenatequilibriumandthereisnobulkregionwherethepolymerassumesitsunconnedglobularform.Thisindicatesananalogywithanisotropicrigidparticleswhichalsomigratetowardstheboundariesunderow.Thismigrationisduetotheparticlealigninginthedirectionofow,therebyreducingthediusivityintheconnementdirection[ 60 61 ].Alignmentalsoreducesthethicknessofthedepletionlayerasaresultofareducedstericforce.Theextensionofthepolymerandthedistributionofpolymerconformationsdependsonchainlength,especiallyathighPecletnumbers[ 62 ].Resolutionofsuchconformationalchangesmayrequireanerdiscretizationofthepolymerchain.Toassesstheaccuracy 69


Averageaxialvelocityofthecenterofmassofpolymersnormalizedbytheaverageuidow.TheresultsarepresentedwithrespecttothepolymerPecletnumberfordierentchannelsizes. oftheresultspresented,wehaverepeatedsomeofthecalculationsusingchainsof32and64beads.Anerdiscretizationallowsformoreconformationalexibilityandatthesametimeimprovestheaccuracyofhydrodynamicinteractions[ 17 ].WescalethesystemvolumeandpressuregradienttomaintainconstantconnementratioandPecletnumber,whichresultsinanerdiscretizationoftheuidgridincomparisontotheradiusofgyrationofthepolymer.LongerchainsatconstantPecletnumberalsoresultinlowerReynoldsnumbers,inverselyproportionaltotheirequilibriumradiusofgyration.TheReynoldsnumbersusedinthisworkrangefrom0:6to2:3fortheshortestchainsinthelargestchannelsand0:15to0:5forthelongestchains.TheresultsshowninFig. 4-4 indicatethatthelevelofdiscretizationwithN=16issucientforthisproblem.Longerchainsshowidenticaldistributionswithinthestatisticalerrors.Polymermigrationalsoimpliesadierentialvelocitybasedonpolymerchainlength[ 36 53 ].TwopolymerchainsofdierentlengthhavedierentpolymerPecletnumbers 70


4-5 showsthemeanpolymervelocitiesUp=hX(t)i=tcomparedtothemeanuidvelocity,Uf=2Umax=3.Thepolymervelocity,hUpi=hUfi,isauniversalfunctionofthepolymerPecletnumber,Pe,athighPecletnumbers.However,atlowowratesthepolymervelocityissolelydeterminedbytheexcludedvolume.Thesimulationtrajectorieswerethenusedtocalculatetherstandsecondcumulantsofthetheaxialdisplacementofthecenterofmassofthechains,hX(t)iandhX(t)2ihX(t)i2.ThemotionofthepolymerswerefoundtobediusiveontimescaleslargerthanH2=DandresultsfortheaxialdispersioncoecientsareshowninFig. 4-6 (solidsymbols).TheTaylor-Ariscalculationofthedispersioncoecientisanasymptoticsolutionoftheconvection-diusionequation,assumingauniformdistributionofparticlesacrossthechannel, 4{1 evenatlowPecletnumbers.NotethatthemolecularcontributiontoDisnegligibleunderallowconditionsstudiedhere.Ifwetakeintoaaccountthenitesizeofthepolymer,theTaylor-Arisresultismodiedbecausethepolymerdoesnotsamplethehighshearregionsnearthewall.ThecenterofmassisrestrictedtopositionsbetweenLdandHLdwhereLdisthethicknessofthedepletionlayer.AssumingauniformdistributionacrossarestrictedchannelwidthLdyHLd,itisstraightforwardtoshowthat=(12Ld=H)6=210.Ifwetakethedepletionlayerthickness,Ld,tobeequaltotheradiusofgyration,Rg,thenweobtainthesolidlineshowninFig. 4-6 ingoodagreementwithnumericalsimulationsatthesameconnementandatlowPecletnumbers.Wehave 71


HydrodynamicdispersioncoecientsofapolymerchaininPoiseuilleow.Numericalsimulationsareshownaslledsymbols,alongwithTaylor-Aristheory(TA)fortracerparticles(Ld=0,dashedline)andfornitesizedparticleswithaxeddepletionlayer(Ld=Rg=H=8,solidline).Theoreticalcalculations,accountingforthedepletionlayerthicknessdeterminedfromsimulation,areshownasopensymbols.Thestandarderrorinthesimulationdatacorrespondstothesizeofthesymbols. foundthatthelevelingandsaturationofthedispersioncoecientinthesesimulationsareduetoinertialeectswhichbecomeprominentonlybeyondPe>4000.OnrerunningsimilarsimulationsatlowerinertiawehavefoundoutthatthedispersioncoecientsarestillgreatlyreducedbutdonotsaturateatthehighPecletlimit.WeshowtheresultsofthesesimulationsinFig. 4-7 .WealsoshowtheeectofinertiaonthecenterofmassdistributioninFig. 4-8 ,butwedonotndagreatimpactonthedistributionapartfromaslightchangeinthedipinthemiddleofthedistribution.Itisoftensuggestedthatthedierentialvelocityofdierentchainlengthscouldbeusedtoseparatechainsofdierentlengths,butTaylordispersiontendstospreadthedistributiontoomuchforittobeausefulseparationtechnique.Herewebrieyexamine 72


Comparisonofdispersioncoecientswithdierentinertia.ThehighestparticleReynoldsnumberforkT=0:1isalmost3:4,whereasthechannelReynoldsnumbersareupto30forthesamesimulations. theconsequencesofthereducedhydrodynamicdispersion.IfwetakethesaturateddispersioncoecientfromFig. 4-6 ,D=D100(H=Rg)2,thentheexpectedwidthoftheaxialdistributionpassingthroughachanneloflengthLisroughly10p 73


Comparisonofcenterofmassdistributionsat3dierenttemperaturesatPe=300.ThethreedierenttemperaturescorrespondtodierentReynoldsnumbers.ThehighestparticleReynoldsnumberatkT=0:1isaround3:4. thatinwiderchannels,wherehydrodynamicinteractionsarestronger.Thehydrodynamicscreeningfromthewallsreducesthemigrationtowardsthecenterlineathighconnementsandalignmentintheowcanthenproducemigrationtowardsthewalls,similartowhatisobservedinsimulationswithoutintra-beadhydrodynamicinteractions[ 40 ].Wehaveconrmedthattheprimarymechanismformigrationtowardsthecenterlineisacombinationofshearandhydrodynamicinteractionsbetweenthewallsandthebeads[ 55 { 57 ].WealsofoundthatthespatiallyvaryingdiusivityproducesadoublemaximuminPoiseuilleow[ 17 36 ],butnotinauniformshearow.Systematicexperimentsincomparableparameterregimeswouldbeaninterestingconrmationoftheresultsofthisstudyandotherrecentsimulations[ 17 36 ].Wehavealsoreportedonthedispersionanddierentialvelocitiesinsuchsystemsalongwithadiscussionofthelimitationsofseparationsbasedonourresults. 74


36 55 57 ].Recentsimulations[ 57 63 ]showthathydrodynamicliftisthedominantmigrationmechanisminpressure-drivenow,ratherthanspatialgradientsinshearrate.Sincethemigrationdependsonchainlength[ 36 63 ],itisinprinciplepossibletofractionatepolymersbasedontheirlength;longerpolymersmigratemorestronglyandthereforeelutefaster.However,Taylordispersionprecludesanysharplypeakeddistribution,sothereisconsiderableinterestinusingcombinationsofowandelectriceldstoimproveseparationeciency[ 64 { 66 ].Weusednumericalsimulationstoinvestigateaexiblepolymerdrivenbyacombinationofuidowandexternalforce.Inthispreliminaryinvestigationwehavenotconsideredtheeectsofcounterionscreening.Wendthatthepolymermigratestowardsthechannelcenterundertheactionofa 75


66 67 ].Thesimilaritiesbetweentheseresultssuggestthathydrodynamicinteractionsinpolyelectrolytesolutionsareonlypartiallyscreened[ 68 ]. 17 ].TheexternaleldandpressuregradientresultintwodierentPecletnumbers:Pe= 11 ].Theuidissimulatedbyauctuatinglattice-Boltzmannequation[ 3 ],whichaccountsquantitativelyforthedissipativeanductuatinghydrodynamicinteractionsbetweenpolymersegments[ 17 ].ThepolymerchainsarediscretizedintoN1segmentsoflengthb,withmostofthesimulationsusingN=16;additionalsimulationswithN=32wereusedtochecktheeectsofchaindiscretization.Previouswork[ 63 ]hasveriedthattheradiusofgyration,Rg,andthediusioncoecient,D,scaleasymptotically(N0:6)atthislevelofdiscretization. 76


CenterofmassdistributionofasinglepolymerinachannelH=Rg=8undertheapplicationofanexternalbodyforce.Halfofthedistributionisshown,withtheboundaryaty=H=0andthecenterofthechannelaty=H=0:5. boundaries.AtPecletnumbersinexcessof100,thistransversemigrationresultsinanonuniformcenterofmassdistribution(Fig. 5-1 ),withthepolymerconcentratedinthemiddleofthechannel.Forsmallforces(Pe<10),Brownianmotionisdominantandthedistributioneventuallybecomesuniform,exceptforasmalldepletionlayernearthechannelwalls.However,ifhydrodynamicinteractionsareneglectedthereisnomigration,andthedistributionisuniformatallPecletnumbers.Thisobservationemphasizesthehydrodynamicoriginofthemigration.InthelimitofinnitePecletnumberthedistributionstillhasanitewidth,whichshowsthatexiblepolymershaveanadditionaldispersivemechanismbesidesBrownianmotion.Thishydrodynamicdispersionisduetothechaoticmotionoftheinteractingpolymersegmentsandissimilartothedispersivemechanisminsettlingsuspensions[ 69 70 ].Itplaysanimportantroleinlimitingthemigration.Apolymerinthevicinityof 77


71 72 ],andthecenterofmassdistributionwouldbeuniform.Ontheotherhandasphericalparticlewouldnotmigrateatallandthedistributionwouldagainbeuniform.FlexibilityanddispersionarethereforecrucialindeterminingthedistributionatlargePecletnumbers.ThehydrodynamicdispersioncanbeaslargeasBrowniandiusion,ascanbeseenbycomparingthewidthoftheinnitePecletdistributionwithanitePecletnumbercase,Pe=1100(Fig. 5-1 ).NotethatthestatisticalerrorsfortheinnitePecletcaseareapproximately10%,sothesedistributionsarestatisticallyindistinguishable.Atheoreticalapproachcanbeusedtogainfurtherinsightintothehydrodynamicoriginofthemigration.WehaveusedakinetictheorysimilartoRef.[ 57 ]andhaveidentiedthreemechanismsthatleadtomigration;twoofthesehavebeennotedpreviously[ 57 66 ],oneisnew.Weanalyzetheevolutionofthedistributionfunctionofadumbbellnearaplanarboundary,incorporatingauniformexternaleldinadditiontothehydrodynamicoweld[ 57 ].Thekeyassumptionisthatthedistributionfunctioncanbefactoredintoaproductofthecenterofmassdistribution,P(y;t),andtheorientationdistribution(qj;y;t),wherethevectorqjconnectsthemonomersandyisthedistanceofthecenterofmassfromtheboundary.Theexpressionforthelateralux(jy)isthenjy=P kT@P @y+PkT yjFsjistheinstantaneousliftvelocity[ 57 ]ofthepolymerandDijistheKirkwooddiusivity.Theexternalforce,Fx,isdirectedalongthecentralaxisofthe 78


Illustrationofdierentmigrationmechanisms.Ineachgurethesolidlinesrepresenttheforceswhilethedashedlinesrepresentthevelocitiesgeneratedbytheseforces.Thevelocitiescanbecalculatedusinganimagesystemforapointsourcenearaplanarboundary.A)Theliftduetoshear:thetensioninthepolymernearasolidboundarygeneratesanetvelocityawayfromtheboundary.B)Rotationduetotheexternaleld:thevelocityeldduetotheexternalforceoneachbeadresultsinarotationaboutthecenterofmassofthepolymer.C)Driftofarotatedpolymer:twobeadsorientedatanangletotheexternalforcesdriftduetothehydrodynamicinteractionsbetweenthebeads. channel.Wehavefurtherassumedthatthedumbbellisconnectedbyalinearspring,Fsj=kqj,andthattheorientationdistributioncanbeexpandedinpowersofthelocalshearrateandexternalforce.Thissimpliesthetime-independentsolutionforthecenterofmassdistributionfunction,P(y),0=@ @y[AP y22+BP y2F2+CPFkT @y]; 57 ],thissolutionforadumbbellnearasinglewallcanbeextendedtoachannelowbysuperposingthesolutionsforthetwowalls.TherstterminEq. 5{2 describestheliftcreatedforthecenterofmassofadumbbellbytheasymmetrichydrodynamicinteractionsbetweenthepolymersegmentsandthewalls[ 57 ].Thelocalshearrategeneratesatensioninthepolymer,whichcoupleswiththelinear(in)contributionto.Thiscouplingcreatesanonuniform 79


55 ]andsimulations[ 17 36 63 ].Thesecondtermdescribestherotationanddriftofthepolymertowardsthechannelcenterundertheapplicationofanexternaleld.Althoughbothtermsarisefromhydrodynamicinteractionswithaboundary,themechanismsaredierentinthemotiontheycreate.Inashearingowthepolymerisintension,andtheopposingforcesgeneratealiftforthecenterofmass[ 57 ](Fig. 5-3 A),butinanexternaleldtheforceonthemonomershasthesamesign,whichcausesthepolymertorotateawayfromtheboundary(Fig. 5-3 B).Subsequentlythepolymerdriftstowardsthecenterundertheactionoftheforce.Finally,thereisacrosstermthatarisesfromacouplingbetweenthestretchingandrotationofthepolymerbytheowandtheexternalforce(Fig. 5-3 C).ThedirectionofthedriftnowdependsonthesignofF,butisindependentoftheboundaries.ThisisthemechanismproposedinRef.[ 66 ],butitdoesnotexplainpolymermigrationinthepresenceofabodyforcealone.Thenumericalsimulationssuggestthatthecenterofmassdistributionacrossthechannel,P(y;Rg;H;D;U),canbewritteninascaledformP(y=H;Pe).Figure 5-3 ashowscenterofmassdistributionofa16beadpolymerchainatvariouschannelwidthsandpolymerPecletnumbers.PolymerswiththesamePecletnumbershavethesamedistributionregardlessofthesizeofthechannel.Figure 5-3 bshowsthatcenter-of-massdistributionisindependentofthenumberofsegmentsusedtodiscretizethechain.However,thekinetictheorydoesnotcapturethesescalings,whichremainsunexplainedatpresent.Theexternalforcescanbeusedinconjunctionwithapressuredrivenowsuchthatanaxialvelocitydierenceiscreatedbetweenmigratingpolymersbasedonsize.Iftheexternalforceandpressuregradientdrivethepolymerinthesamedirection,thetheorypredictsthatthecoupledtermenhancesthemigrationtowardsthechannelcenter(Eq. 5{2 ).ThenumericalresultsinFig. 5-4 verifythisprediction,ascanbeseenbycomparingtheresults(Pe=110)inFig. 5-1 withFig. 5-4 (Pef=12:5).Sincetheow 80


a)ScalingofthedistributionfunctionwithrespecttothechannelwidthunderconstantpolymerPecletnumber,Pe.b)ScalingofthedistributionfunctionwithrespecttopolymerPecletnumberusingdierentlevelsofchaindiscretization.Theboundaryisaty=H=0andthecenterofthechannelisaty=H=0:5. isrelativelyweak(Pef10),itmakesonlya1020%contributiontothemigration.However,amixtureofpolymersinsuchasystemhassignicantlydierentelutionrates;largerpolymersmigratemorestronglyandthereforesamplehighervelocitystreamlinesnearthecenterofthechannel.Sincethemigrationisprimarilydrivenbytheexternaleld,aweakhydrodynamicowcouldbeusedforfractionation,therebyreducingtheTaylordispersion.Inacountercurrentapplicationofthetwoeldsthepolymertendstoorientitselfindierentquadrants,dependingontherelativemagnitudeofthetwodrivingforces.Thepolymernowdriftseithertowardsthewallsortowardsthecenterdependingonitsmeanorientation.TheresultsinFig. 5-5 showmigrationtowardstheboundarieswhen 81


Centerofmassdistributionsforconcurrentapplicationofanexternalbodyforceandpressure-drivenow.Thesolidcurveshowsthelevelofmigrationunderthepressure-drivenowonly.TheowPecletnumberinallthreecasesisPef=12:5.Theboundaryisaty=H=0andthecenterofthechannelisaty=H=0:5. theexternalforceissmall(Pe<30),butincreasingtheforceeventuallyreversestheorientationofthepolymerandthepolymeragainmigratestowardsthecenter(Pe>100).Wenoteasanempiricalobservationthatthepeakofthepolymermigrationtowardsthewall(Pe=44)occurswhenthepolymerisdrivenwithalmostequalbutoppositevelocitybytheowandforce.ThenumericalresultsandkinetictheorypredictionsareinqualitativeagreementwithaseriesofexperimentsmeasuringthedistributionofconnedDNAunderthecombinedactionofanelectriceldandapressure-drivenow[ 66 67 ].Theseexperimentsshowedastrongmigrationtowardsthecenterintheconcurrentapplicationofforceandowandalsothereversemigrationtowardstheboundariesinthecountercurrentapplication. 82


Centerofmassdistributionsforcountercurrentapplicationofanexternalbodyforceandpressure-drivenow.Thesolidcurveshowsthelevelofmigrationunderthepressure-drivenowonly.TheowPecletnumberinallcasesisPef=12:5.Theboundaryisaty=H=0andthecenterofthechannelisaty=H=0:5. Howevertheexperimentsdonotshowastrongmigrationwhentheelectriceldisappliedbyitself.Thissuggeststhatthedominantmechanisminthelaboratoryexperiments[ 66 ]isthecouplingbetweenthelocalshear-rateandtheappliedforce(Eq. 5{2 ),ratherthanrotationofthepolymerbyhydrodynamicinteractionswiththewalls.Thechanneldimensionsintheexperimentswereabout10timeslarger,incomparisonwiththepolymersize,thanthoseinthesimulations,whichreducestheimportanceofwalleects.Morecrucially,thehydrodynamicinteractionsinpolyelectrolytesolutionsdrivenbyanelectriceldarescreenedbythecounterionmotion[ 73 74 ].However,ifthereisnosignicanthydrodynamicinteractionbetweensegments,thenthepolymerwillfollowthedirectionoftheexternalforcewithoutsignicanttransversemotion,evenwhendeformedbya 83


67 ],onlymigratesinanexternalforcebecauseofhydrodynamicinteractionsbetweensegments.Wehaveveriedthisbehaviorwithnumericalsimulationsinwhichhydrodynamicinteractionswereexcluded.Thustheobservationofmigrationintheseexperiments[ 66 67 ],impliesadegreeofhydrodynamicinteractiononthescaleofthepolymer.Ithasbeenpredictedthattheuidvelocityaroundapolyelectrolytehasadipolarcomponentinanelectriceld,decayingas1=r3[ 68 75 ].Thisisinadditiontotheangle-averagedvelocityeld,whichdecaysexponentially[ 73 ].Forananisotropicpolyelectrolyte,distortedbytheshearow,thedipolarcontributiongivesrisetoatransversevelocityofmagnitudeF=2R3g[ 68 ],whereFisthetotalelectricforceonthepolymer,istheuidviscosity,andistheinverseDebyelength.ThisvelocityissmallcomparedwiththeelectrophoreticvelocityF=L;for-DNAin5mMsaltsolution,themigrationvelocityisthenpredictedtobeoftheorderof1%oftheelectrophoreticvelocity.Thisisconsistentwithexperimentalmeasurements(E=40V/cm)[ 66 ],whichsuggestmigrationvelocitiesoftheorderof103cm/sec,incomparisonwithelectrophoreticvelocitiesoftheorderof101cm/sec.ThetheoreticalestimatescouldbemademorequantitativeusingthekinetictheoryoutlinedinthispaperanddescribedinmoredetailinRef.[ 57 ].Additionalexperimentscouldbeperformedtoprobethescreeninglengthassociatedwiththehydrodynamicinteractionsinpolyelectrolytesolutions.WenotethatFig.3inRef.[ 66 ]mayindicateaweakforce-inducedmigration,althoughthestatisticalerrorsinthedatamakeitinconclusive.Anarrowerchannelwouldbehelpfulinpromotinghydrodynamicinteractionsbetweenthepolymerandthechannelwallsandtherebystrengtheninganyforce-inducedmigration.Themagnitudeoftheappliedforcecouldbeincreasedbyusinglow-frequencyACelds;sincetheforce-inducedmigrationisquadraticintheeld,thedirectionofmigrationisindependentofthesign.Thisshouldmakeitpossibletoclearlydecideifthereisanyforce-inducedmigrationornot.Screeningcould 84


76 77 ]. 66 67 ].Thesimilaritiesbetweenthemigrationofaneutralpolymerdrivenbyabodyforceandapolyelectrolytedrivenbyanelectriceldsuggestthathydrodynamicinteractionsinpolyelectrolytesolutionsareonlypartiallyscreened[ 68 ].Indeed,iftheywerefullyscreened[ 73 ]theobservedmigrationwouldnotoccur.Wehavesuggestedadditionallaboratoryandnumericalexperimentstofurtherinvestigatethisquestion. 85


36 ]andakinetictheoryforadumbbell[ 57 ]andhavesupportedthemechanismthatthemigrationisaresultofasymmetrichydrodynamicinteractionswiththewall.Insuchowswehaveshownthatmigrationcanbetowardstheboundariesratherthanthecenterinverynarrowchannelsduetothescreeningofhydrodynamicinteractions.Wewerealsoabletoclarifythecorrectsourceofmigrationmechanismbysimulationsofasimpleshearow.Weprovedthat 86




3-3 ).Nevertheless,theuctuation-dissipationrelationbetweentheeectivefrictionanddiusioncoecientsholdswithinstatisticalerrors(Fig. 3-1 ).Ithasbeenshown[ 21 ]thatuctuationsinadiscretesystemcanbedierentfromthoseincontinuoussystems(c.f.Eqs. 2{15 and 2{16 ),andinordertounderstandwhydeviationsfromclassicalequipartitiondonotaectthelong-timedynamics,wehereexamineasimplemodelwithconstantfriction.ThedynamicsaredescribedbytheclassicalLangevinequation, dt=U+R(t);(A{1)withavarianceintherandomforcehR(t)R(t0)i=2T(tt0).Withthischoiceofrandomforce,Eq. A{1 satisesboththeStokes-EinsteinrelationD=T=andequipartitionmhU2i=T.Therealsituationismorecomplex[ 20 33 ],withatime-dependentfrictioncoecientduetothetemporalandspatialdevelopmentofthehydrodynamicoweld.Nevertheless,thissimplemodelissucienttoillustratethedeviationsinmeankineticenergyfromclassicalequipartition.ThediscreteequivalentofEq. A{1 canbewrittenas 78 ] U(t)=1 88


A{2 hasthesolution ;(A{6)satisestheStokes-Einsteinrelation.Thuswerecoverthecorrectlong-timediusivebehavior,despitethediscretetimestep,whereasthekineticenergyhasanadditionalfactor(1=2)1.Anaccurateequipartitionofenergyrequiresthat1,whichinourcontextmeansthatthemassofthemonomermustbelarge.However,numericalcalculationswith<0:5showthatthediusivebehaviorofinteractingparticlesisunaectedbyviolationsofequipartition.AnimplicitupdateofEq. A{1 A{2 ).ItisstraightforwardtoshowthattheStokes-EinsteinrelationD=T=isagainsatised,butinthiscasethemeankineticenergyisgivenby A{5 and A{8 withthenumericalsimulationsisshowninFig. 3-3 .Itcanbeseenthatthesimpliedmodelpresentedhere,whichneglectsthe 89




33 ].WedeneCF(t)=,theforceautocorrelationfunctionanddeneCF(w)asthepowerspectrumofFwhichistheFouriertransformoftheforceautocorrelationfunction.Wecannowwrite 2Z1eiwtCF(w)dw;(B{3)thenequation B{1 becomes Z1Zt0CF(w)eiwt(1t 91


iw2tsin(wt) 33 ]:CF(w)=kT[(w)+(w)] (B{10)(w)=ij6R[1+R(iw=)1=2(iwR2=9)]: 92


(C{4) Wethenfollowthesameprocedureforthesphere(seeAppendix B )tocalculateandtheresultingequationfortheplanewallis, t1 2:(C{6)Againit'simportanttonotethisresultshouldbeintegratedoncemoreoverthesurfaceofinterest. 93


Px=2(n7n8)t Py=2(n15n17)t Thepopulationdensitiesaregivenbythepostcollisiondistribution 2c4s;(D{3)sincethereisnoowandonlyuctuations,uucanbeneglectedwhichleavesuswith 2c4s;(D{4)Herewewillcalculatethemomentumchangeinxdirection.Thepopulationscontributingtothisdirectionhaveavelocityofp 36.Thespeedofsoundc2s=Thepostcollisionmomentumuxneq0setting=1isonlytherandomstresstermthus,neq0=f.Thesegive (D{5) (D{6) 94


36+(D{7)Since=x3kT[ 2 ]andisathirdofthis,andusingtheequation 2{31 ,andthat=36inourmodel,wend. 3kT+12kT(D{8)Thisisthemomentumchangecorrelationinxdirectionperunitarea,andthesameexpressionisfoundforthezdirectionaswell. 95


[1] C.W.J.BeenakkerandP.Mazur,\Self-diusionofspheresinaconcentratedsuspension,"PhysicaA,vol.120,pp.388{410,1983. [2] A.J.C.LaddandR.Verberg,\Lattice-Boltzmannsimulationsofparticle-uidsuspensions,"J.Stat.Phys.,vol.56,pp.286{311,2001. [3] A.J.C.Ladd,\Short-timemotionofcolloidalparticles:Numericalsimulationviaauctuatinglattice-Boltzmannequation,"Phys.Rev.Lett.,vol.70,pp.1339,1993. [4] A.J.C.Ladd,\NumericalsimulationsofparticulatesuspensionsviaadiscretizedBoltzmannequationPartI.Theoreticalfoundation,"J.FluidMech.,vol.271,pp.285,1994. [5] A.J.C.Ladd,\Hydrodynamicinteractionsinasuspensionofsphericalparticles,"J.Chem.Phys.,vol.88,pp.5051,1988. [6] G.K.Batchelor,AnIntroductiontoFluidDynamics,CambridgeUniversityPress,Cambridge,1967. [7] J.HappelandH.Brenner,Low-ReynoldsNumberHydrodynamics,MartinusNijho,Dordrecht,1986. [8] H.J.MasliyahandS.Bhattacharjee,ElectrokineticandColloidTransportPhenom-ena,Wiley,NewYork,2006. [9] A.J.C.Ladd,M.E.Colvin,andD.Frenkel,\Applicationoflattice-gascellularautomatatotheBrownianmotionofsolidsinsuspension,"Phys.Rev.Lett.,vol.60,pp.975,1988. [10] P.MazurandW.V.Saarloos,\Many-spherehydrodymamicinteractionsandmobilitiesinasuspension,"PhysicaA,vol.115,pp.21{57,1982. [11] P.AhlrichsandB.Dunweg,\SimulationofasinglepolymerchaininsolutionbycombininglatticeBoltzmannandmoleculardynamics,"J.Chem.Phys.,vol.111,no.17,pp.8225{8239,1999. [12] B.J.AlderandT.E.Wainwright,\Decayofthevelocityautocorrelationfunction,"Phys.Rev.A,vol.1,no.1,pp.18{21,1970. [13] T.Phung,J.Brady,andG.Bossis,\Stokesiandynamcissimulationofbrowniansuspensions,"J.FluidMech.,vol.313,pp.181{207,1996. [14] J.P.BoonandS.Yip,MolecularHydrodynamics,Mc-GrawHill,NewYork,1980. [15] A.J.C.Ladd,\NumericalsimulationsofparticulatesuspensionsviaadiscretizedBoltzmannequationPartII.Numericalresults,"J.FluidMech.,vol.271,pp.311,1994. 96


N.Q.NguyenandA.J.C.Ladd,\Lubricationcorrectionsforlattice-boltzmannsimulationsofparticlesuspensions,"Phys.Rev.E,vol.66,no.046708,2002. [17] O.B.Usta,A.J.C.Ladd,andJ.E.Butler,\Lattice-Boltzmannsimulationsofthedynamicsofpolymersolutionsinperiodicandconnedgeometries,"J.Chem.Phys.,vol.122,pp.094902{1,2005. [18] S.Succi,TheLatticeBotlzmannEquationforFluidDynamicsandBeyond,OxfordUniversityPress,2001. [19] K.Huang,StatisticalMechanics,Wiley,NewYork,1987. [20] L.D.LandauandE.M.Lifshitz,FluidMechnaics,OxfordUniversityPress,NewYork,1966. [21] A.J.C.LaddandR.Verberg,\Lattice-Boltzmannsimulationsofparticle-uidsuspensions,"J.Stat.Phys.,vol.104,pp.1191{1251,2001. [22] A.J.C.Ladd,H.Gang,J.X.Zhu,andD.A.Weitz,\Time-dependentcollectivediusionofcolloidalparticles,"Phys.Rev.Lett.,vol.74,pp.318,1995. [23] A.J.C.Ladd,\Hydrodynamicscreeninginsedimentingsuspensionsofnon-Brownianspheres,"Phys.Rev.Lett.,vol.76,pp.1392,1996. [24] A.J.C.Ladd,\Eectsofcontainerwallsonthevelocityuctuationsofsedimentingspheres,"Phys.Rev.Lett.,vol.88,pp.048301,2002. [25] C.K.Aidun,Y.N.Lu,andE.Ding,\DirectanalysisofparticulatesuspensionswithinertiausingthediscreteBoltzmannequation,"J.FluidMech.,vol.373,pp.287{311,1998. [26] F.M.Auzerais,J.Dunsmuir,B.B.Ferreol,N.Martys,J.Olson,T.S.Ramakrishnan,D.H.Rothman,andL.M.Schwartz,\Transportinsandstone:Astudybasedonthreedimensionalmicrotomography,"Geophys.Res.Lett.,vol.23,pp.705{108,1996. [27] N.S.MartysandH.Chen,\Simulationofmulticomponentuidsincomplexthree-dimensionalgeometriesbythelatticeboltzmannmethod,"Phys.Rev.E,vol.53,no.1,pp.743{750,Jan1996. [28] L.W.P.Manz,B;Gladden,\Flowanddispersioninporousmedia:Lattice-boltzmannandnmrstudies,"AicheJ.,vol.45,pp.1845{1854,1999. [29] P.Ahlrichs,R.Everaers,andB.Dunweg,\SimulationofasinglepolymerchaininsolutionbycombininglatticeBoltzmannandmoleculardynamics,"Phys.Rev.E,vol.64,pp.040501(R),2001. [30] U.Frisch,D.d'Humieres,B.Hasslacher,P.Lallemand,Y.Pomeau,andJ.-P.Rivet,\Latticegashydrodynamicsintwoandthreedimensions,"ComplexSystems,vol.1,pp.649,1987. 97


D.L.ErmakandJ.A.McCammon,\Browniandynamicswithhydrodynamicinteractions,"J.Chem.Phys.,vol.69,pp.1352,1978. [32] X.Fan,N.Phan-Thien,N.T.Yong,X.Wu,andD.Xu,\Microchannelowofamacromolecularsuspension,"Phys.Fluids,vol.15,no.1,pp.11{21,2003. [33] E.H.HaugeandA.Martin-Loef,\Fluctuatinghydrodynamicsandbrownianmotion,"JournalofStatisticalPhysics,vol.7,no.3,pp.259{281,1972. [34] A.J.C.Ladd,\Dynamicalsimulationsofsedimentingspheres,"Phys.FluidsA,vol.5,pp.299,1993. [35] R.Jendrejack,M.Graham,andJ.dePablo,\EectofconnementonDNAdynamicsinmicrouidicdevices,"J.Chem.Phys,vol.119,pp.1165{1173,2003. [36] R.M.Jendrejack,D.C.Schwartz,J.J.dePablo,andM.D.Graham,\Shear-inducedmigrationinowingpolymersolutions:Simulationoflong-chainDNAinmicrochannels,"J.Chem.Phys.,vol.120,pp.2513{2529,2004. [37] M.Fixman,\ConstructionoftheLangevinforcesinthesimulationofhydrodynamicinteraction,"Macromolecules,vol.19,pp.1204{1207,1986. [38] R.Jendrejack,M.Graham,andJ.dePablo,\Hydrodynamicinteractionsinlongchainpolymers:Applicationofthechebyshevpolynomialapproximationinstochasticsimulations,"J.Chem.Phys,vol.113,pp.2894{2900,2000. [39] N.J.Woo,E.Shaqfeh,andB.Khomami,\EectofconnementondynamicsandrheologyofdiluteDNAsolutions.I.Entropicspringforceunderconnementandnumericalalgorithm,"J.Rheol,vol.48,pp.281{298,2004. [40] N.J.Woo,E.Shaqfeh,andB.Khomami,\EectofconnementondynamicsandrheologyofdiluteDNAsolutions.II.Eectiverheologyandsinglechaindynamics,"J.Rheol,vol.48,pp.299{318,2004. [41] C.-C.HsiehandR.G.Larson,\ModelinghydrodynamicinteractioninBrowniandynamics:Simulationsofextensionalandshearowsofdilutesolutionsofhighmolecularweightpolystyrene,"J.Stat.Phys,vol.36,pp.717,1984. [42] Y.OonoandM.Kohmoto,\Renormalizationgrouptheoryoftransportpropertiesofpolymersolutions.i.Dilutesolutions,"J.Chem.Phys.,vol.78,pp.520{528,1983. [43] H.Hasimoto,\OntheperiodicfundamentalsolutionsoftheStokesequationsandtheirapplicationtoviscousowpastacubicarrayofspheres,"J.FluidMech.,vol.5,pp.317,1959. [44] A.J.C.Ladd,\Hydrodynamictransportcoecientsofrandomdispersionsofhardspheres,"J.Chem.Phys.,vol.93,pp.3484,1990. 98


T.M.SquiresandM.P.Brenner,\Like-chargeattractionandhydrodynamicinteraction,"Phys.Rev.Lett.,vol.85,no.23,pp.4976{4979,Dec2000. [46] F.BrochardandP.G.deGennes,\Dynamicsofconnedpolymers,"J.Chem.Phys.,vol.67,pp.52,1977. [47] H.J.L.andD.M.,\Diusionofmacromoleculesinnarrowcapillaries,"j.Phys.Chem.,vol.96,pp.4046{4052,1992. [48] S.AsakuraandF.Oosowa,\Oninteractionbetweentwobodiesimmersedinasolutionofmacromolecules,"J.Chem.Phys.,vol.22,pp.1225{1256,1954. [49] J.H.AubertandM.Tirrel,\Eectiveviscosityofdilutepolymersolutionsnearconningboundaries,"J.Chem.Phys.,vol.77,pp.553{561,1982. [50] D.Aussere,H.Hervet,andF.Rondelez,\Concentrationdependenceoftheinterfacialdepletionlayerthicknessforpolymersolutionsincontactwithnonadsorbingwalls,"Macromolecules,vol.19,pp.85{88,1986. [51] W.F.Busse,\Twodecadesofhigh-polymerphysics-asurveyandforecast,"Phys.Today,vol.17,pp.32,1964. [52] A.B.Metzner,Y.Cohen,andC.Rangel-Nafaile,\Inhomogeneousowsofnon-Newtonianuids:Generationofspatialconcentrationgradients,"J.Non-NewtonianFluidMech.,vol.5,pp.449{462,1979. [53] U.S.Agarwall,A.Dutta,andR.A.Mashelkar,\Migrationofmacromoleculesunderow:Thephysicaloriginandengineeringapplications,"Chem.Eng.Sci.,vol.49,pp.1693{1717,1994. [54] K.A.DillandB.H.Zimm,\RheologicalseparatorforverylargeDNAmolecules,"NucleicAcidsResearch,vol.7,pp.735{749,1979. [55] L.Fang,H.Hu,andR.G.Larson,\DNAcongurationsandconcentrationinshearingownearaglasssurfaceinamicrochannel,"J.Rheol.,vol.49,no.1,pp.127{138,2005. [56] M.S.JhonandK.F.Freed,\PolymermigrationinNewtonianuids,"J.Polym.Sci.PartB,vol.23,pp.955{971,1985. [57] H.MaandM.Graham,\Theoryofshearinducedmigrationindilutepolymersolutionsnearsolidboundaries,"Phys.Fluids.,vol.17,pp.083103,2005. [58] M.Fixman,\Eectsofuctuatinghydrodynamicinteraction,"J.Chem.Phys,vol.78,pp.1594,1983. [59] N.LironandS.Mochon,\StokesowforaStokesletbetweentwoparallelatplates,"J.Eng.Math.,vol.10,pp.287{303,1976. 99

PAGE 100

L.C.NitscheandE.J.Hinch,\ShearinducedlateralmigrationofBrownianrigidrodsinparabolicchannelow,"J.FluidMech.,vol.332,pp.1{21,97. [61] R.L.SchiekandE.S.G.Shaqfeh,\Cross-streamlinemigrationofslenderBrownianbresinplanePoiseuilleow,"J.Fluid.Mech,vol.332,pp.23{39,1997. [62] R.G.Larson,\Therheologyofdilutesolutionsofexiblepolymers:Progressandproblems,"J.Rheol.,vol.49,no.1,pp.1{70,2005. [63] O.B.Usta,J.E.Butler,andA.J.C.Ladd,\Flow-inducedmigrationofpolymersindilutesolution,"Phys.Fluids.,vol.18,pp.031703,2006. [64] J.-L.Viovy,\Electrophoresisofdnaandotherpolyelectrolytes:Physicalmechanisms,"Rev.Mod.Phys.,vol.72,no.3,pp.813{872,Jul2000. [65] Z.ChenandA.Chauhan,\DNAseparationbyEFFFinamicrochannel,"J.Coll.Int.Sci.,vol.,2005. [66] J.ZhengandE.S.Yeung,\AnomalousradialmigrationofsingleDNAmoleculesincapillaryelectrophoresis,"Anal.Chem.,vol.74,pp.4536{4547,2002. [67] J.ZhengandE.S.Yeung,\Mechanismfortheseparationoflargemoleculesbasedonradialmigrationincapillaryelectrophoresis,"Anal.Chem.,vol.75,pp.3675,2003. [68] D.LongandA.Ajdari,\Anoteonthescreeningofhydrodynamicinteractions,inelectrophoresis,andinporousmedia,"EuroPhys.J.E,vol.4,pp.29{32,2001. [69] R.E.CaischandJ.H.C.Luke,\Varianceinthesedimentationspeedofasuspension,"Phys.Fluids,vol.28,pp.759,1985. [70] J.Martin,N.Rakotomalala,andD.Salin,\HydrodynamicDispersionofNoncolloidalSuspensions:MeasurementfromEinstein'sArgument,"Phys.Rev.Lett.,vol.74,pp.1347,1995. [71] W.B.Russel,E.J.Hinch,L.G.Leal,andG.Tiefenbruck,\Rodsfallingnearaverticalwall,"J.Fluid.Mech,vol.83,pp.273{287,1977. [72] T.N.Swaminathan,K.Mukundakrishnan,andH.Hu,\Sedimentationofanellipsoidinsideaninnitelylongtubeatlowandintermediatereynoldsnumbers,"J.FluidMech.,vol.551,pp.357{385,2006. [73] G.S.Manning,\Limitinglawsandcounterioncondensationinpolylelectrolytesolutions7.electrophoreticmobilityandconductance,"J.Phys.Chem,vol.85,pp.1506,1981. [74] J.BarratandJ.Joanny,\Theoryofpolyelectrolytesolutions,"Adv.Chem.Phys.,vol.94,pp.1{66,1994. 100

PAGE 101

D.Long,J.-L.Viovy,andA.Ajdari,\Simultaneousactionofelectriceldsandnonelectricforcesonapolyelectrolyte:Motionanddeformation,"Phys.Rev.Lett.,vol.76,no.20,pp.3858{3861,May1996. [76] J.HorbachandD.Frenkel,\Lattice-boltzmannmethodforthesimulationoftransportphenomenainchargedcolloids,"Phys.Rev.E,vol.64,no.6,pp.061507,Nov2001. [77] K.Kim,Y.Nakayama,andR.Yamamoto,\Directnumericalsimulationsofelectrophoresisofchargedcolloids,"Phys.Rev.Lett.,vol.96,no.20,pp.208302,2006. [78] H.C.Ottinger,StochasticProcessesinPolymericFluids,Springer-Verlag,Berlin,1996. 101

PAGE 102

O.BerkUstawasborninAdana,TurkeyandwasraisedinMersin,TurkeywhereheattendedtheIcelAnadoluHighSchool.HethenmovedtoIstanbulwherehegothisBachelorsdegreeinchemicalengineeringfromtheBogaziciUniversity.Sincethenhe'sbeenworkingonhisdoctoraldegreeinchemicalengineeringattheUniversityofFlorida.BerkwillbenishinginDecember2006.He'sthankfulforthebeautifulnatureandatmosphereprovidedtohimforhisstudiesinFlorida.Herehealsotookupphotographyduetothedazzlingbeautyofeverythingaroundhim. 102