UFDC Home  myUFDC Home  Help 
Title Page  
Dedication  
Acknowledgement  
Table of Contents  
List of Figures  
Abstract  
Complex fluids and hydrodynami...  
Latticeboltzmann methods for colloidal...  
Dynamics of polymer solutions in...  
Flowinduced migration of polymers...  
Transverse migration of polymers...  
Conclusions  
Appendices  
References  
Biographical sketch 



Table of Contents  
Title Page
Page 1 Page 2 Dedication Page 3 Acknowledgement Page 4 Page 5 Table of Contents Page 6 Page 7 List of Figures Page 8 Page 9 Page 10 Page 11 Abstract Page 12 Page 13 Complex fluids and hydrodynamics Page 14 Page 15 Page 16 Page 17 Page 18 Page 19 Page 20 Page 21 Page 22 Page 23 Latticeboltzmann methods for colloidal and polymeric suspensions Page 24 Page 25 Page 26 Page 27 Page 28 Page 29 Page 30 Page 31 Page 32 Page 33 Page 34 Page 35 Page 36 Page 37 Page 38 Page 39 Page 40 Page 41 Page 42 Dynamics of polymer solutions in periodic and confined geometries Page 43 Page 44 Page 45 Page 46 Page 47 Page 48 Page 49 Page 50 Page 51 Page 52 Page 53 Page 54 Page 55 Page 56 Page 57 Page 58 Page 59 Page 60 Page 61 Flowinduced migration of polymers in dilute solution Page 62 Page 63 Page 64 Page 65 Page 66 Page 67 Page 68 Page 69 Page 70 Page 71 Page 72 Page 73 Page 74 Transverse migration of polymers in external fields Page 75 Page 76 Page 77 Page 78 Page 79 Page 80 Page 81 Page 82 Page 83 Page 84 Page 85 Conclusions Page 86 Page 87 Appendices Page 88 Page 89 Page 90 Page 91 Page 92 Page 93 Page 94 Page 95 References Page 96 Page 97 Page 98 Page 99 Page 100 Page 101 Biographical sketch Page 102 

Full Text  
DYNAMICS OF COLLOIDAL AND POLYMERIC SUSPENSIONS VIA A FLUCTUATING LATTICEBOLTZMANN MODEL By 0. BERK USTA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2007 Copyright 2007 by 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 ACKNOWLEDGMENTS 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, HyunOk 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, ChiChuing Li, Vicky Huang, Darren McDuff, Colin Sturm, Young Sun, YoungSeok Kim, ChiaYi 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. TABLE OF CONTENTS page ACKNOW LEDGMENTS ................................. 4 LIST OF FIGURES .. .. .. ... .. .. .. .. ... .. .. .. ... .. .. .. 8 A B ST R A C T . . . . . . . . . . 12 CHAPTER 1 COMPLEX FLUIDS AND HYDRODYNAMICS ................. 14 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 2 LATTICEBOLTZMANN METHODS FOR COLLOIDAL AND POLYMERIC SU SPEN SIO N S . . . . . . . . . 24 2.1 Boltzm ann Equation .............................. 24 2.2 Lattice Boltzmann Method ........................... 26 2.2.1 An Overview of LatticeBoltzmann Method ............. 26 2.2.2 Lattice Boltzmann Equation ...................... 28 2.3 Fluctuations . . . . . . . . 30 2.4 Solid Particles . . . . . . . . 31 2.4.1 FiniteSize Spherical Particles . . . . . 31 2.4.2 Point Particles . . . . . . . 33 2.4.3 Particle M otion . . . . . . . 36 2.5 Verification of the Fluctuations . . . . . . 37 3 DYNAMICS OF POLYMER SOLUTIONS IN PERIODIC AND CONFINED GEOMETRIES . . . . . . . . . 43 3.1 Introduction . . . . . . . . 43 3.2 R results . . . . . . . . 45 3.2.1 FluctuationDissipation . . . . . . 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 FLOWINDUCED MIGRATION OF POLYMERS IN DILUTE SOLUTION 62 4.1 Introduction . . . . . . . . 62 4.2 M ethod . . . . . . . . . 63 4.3 Results and Discussion . . . . . . . 65 4.4 Conclusion . . . . . . . . . 73 5 TRANSVERSE MIGRATION OF POLYMERS IN EXTERNAL FIELDS ... 75 5.1 Introduction . . . . . . . . 75 5.2 M ethods . . . . . . . . . 76 5.3 R results . . . . . . . . . 76 5.4 C conclusion . . . . . . . . . 85 6 CONCLUSIONS .. .. .. ... .. .. .. ... .. .. .. .. ... .. .. 86 APPENDIX A DISCRETE LANGEVIN EQUATION ....................... 88 B THE MOMENTUM CHANGE CORRELATION ON A SPHERE ........ 91 C THE MOMENTUM CHANGE CORRELATION ON A PLANE WALL, CONTINUUM T H E O RY . . . . . . . . . 93 D THE MOMENTUM CHANGE CORRELATION ON PLANE WALL, DISCRETE T H E O RY . . . . . . . . . 94 REFERENCES ....................................... 96 BIOGRAPHICAL SKETCH ................................ 102 LIST OF FIGURES Figure page 21 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 22 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 timestep, error accumulates for a couple of timesteps 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 23 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 24 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 25 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 31 Relation between the effective friction coefficient and the input friction o (Eq. 221). 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 32 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 33 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 34 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 latticeBoltzmann simulations (circles) in the sam e size cell . . . . . . . . 51 35 Comparison of the parallel component of the two particle mobility tensor, pl2. Results obtained by dissipative (Eq. 310, circles) and fluctuating (Eq. 311, squares) simulations are compared with an independent solution of the periodic Stokes equations (Eq. 39, triangles). . . . . . . 52 36 Comparison of the perpendicular component of the two particle mobility tensor, pf2. Results obtained by dissipative (Eq. 310, circles) and fluctuating (Eq. 311, squares) simulations are compared with an independent solution of the periodic Stokes equations (Eq. 39, triangles). . . . . . . 53 37 Endtoend 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 38 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 39 Center of mass diffusion coefficient, with finitesize 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 310 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 311 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 41 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 42 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 43 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 44 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 45 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 46 Hydrodynamic dispersion coefficients of a polymer chain in Poiseuille flow. Numerical simulations are shown as filled symbols, along with TaylorAris 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 l1 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 47 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 48 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 51 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 52 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 .1iv, 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 53 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 54 Center of mass distributions for concurrent application of an external body force and pressuredriven flow. The solid curve shows the level of migration under the pressuredriven 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 55 Center of mass distributions for countercurrent application of an external body force and pressuredriven flow. The solid curve shows the level of migration under the pressuredriven 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 DYNAMICS OF COLLOIDAL AND POLYMERIC SUSPENSIONS VIA A FLUCTUATING LATTICEBOLTZMANN MODEL By 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 latticeBoltzmann fluid. This results in the solution of the NavierStokes 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 Stokeslike 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 latticeBoltzmann 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. CHAPTER 1 COMPLEX FLUIDS AND HYDRODYNAMICS 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 interparticle interactions such as electrostatic interactions and hardsphere collisions which also exist in simple fluids. In addition to these, submicron 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 ,niv1 body, therefore they create a great deal of complexity for theoretical calculations. In this work we investigate the static and dynamic transport properties of hardsphere colloidal suspensions and dilute polymer solutions with a numerical simulation method. This method combines Newtonian dynamics of solid particles with a fluctuating latticeBoltzmann fluid [24]. 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 latticegas 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 latticeBoltzmann method. This includes two different approaches for the solid particles coupled to the latticeBoltzmann fluid for hardsphere 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 NavierStokes 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 equation: Op ( ) = V.(pv), (11) 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. 11 reduces to V v 0. (12) 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 outflux of momentum through its surface. [7]. Therefore the momentum balance equation reads, a + V (pvv) V H + pF, (13) 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)'. (14) 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. 13 gives the NavierStokes equation: p (O + v Vv Vp + /V2v + pF. (15) 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 [7]: pV2v = Vp (16) V v 0. (17) The Stokes equations constitutes a great simplification over NavierStokes 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), (18) 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. Noslip (stick) boundary conditions apply at the solidfluid 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. 18, is v(r) U + x (r R) for r R = a, (19) 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 (, N Fi  U + Q] (110) [13 Uj i (1310) j1 N Ti I j + Ci Qj]. (111) j1 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, N Ui [ TFj + F Tj] (112) j1 N i  [i Fj + Tj]. (113) j1 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 longranged, 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 r1 and from the n iiv1 body nature of the interactions. When the particles in the suspension are small enough such that they are much smaller than the average interparticle 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 (114) where H is the Oseen tensor. This tensor is given by H I + ]. (115) 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 shearrate 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, (116) 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 timecorrelation 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 timescale separation. Stokesian dynamics is an accurate simulation technique which uses a reformulation of the Stokes equations (Eq. 16) 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 fastmultipole 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 nonspherical 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 longtime 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, (117) where H denotes hydrodynamic forces, B Brownian forces, and NH nonhydrodynamic 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: dr. where ( is the friction coefficient and ur(ri) is the unperturbed fluid velocity at the position ri of the particle. Equation 118 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 configurationdependent 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 NavierStokes equations, latticegas 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 scales. In this study we have chosen to study a coupling of a latticeBoltzmann 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. CHAPTER 2 LATTICEBOLTZMANN METHODS FOR COLLOIDAL AND POLYMERIC SUSPENSIONS The latticeBoltzmann method is a basis for simulating solidfluid 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 latticeBoltzmann 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, (21) 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(twobody) collisions taking molecules in and out of the streaming trajectory [18]. 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 threebody 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 pointlike, structureless molecules interacting via a shortri,. twobody 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 (22) 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 (23) 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 veryshort 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 fluidparticle interactions. One important example is the longtime 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. (24) The left hand side is similar to Newtonian singleparticle 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,. rStokes equations. 2.2 Lattice Boltzmann Method 2.2.1 An Overview of LatticeBoltzmann 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 preaveraging 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 fluctuationdissipation theorem. The computational method I present in this Ch'! pter is based on a fluctuating latticeBoltzmann 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 finitesize particles can be reproduced quantitatively [4]. Computationally, the latticeBoltzmann method is inherently and straightforwardly parallelizable, enabling large simulations to be spread across a cluster of processors with a relatively low network bandwidth. The latticeBoltzmann method has been used for a wide variety of complex flows, including suspensions of colloidal [3, 22] and noncolloidal [2325] particles, and porous media [2628]. 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 latticeBoltzmann model and investigated the static and dynamical scaling laws of individual polymers and semidilute solutions [11, 29]. 2.2.2 Lattice Boltzmann Equation The fundamental quantity in the latticeBoltzmann model is the discretized oneparticle 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. (25) 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)], (26) 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 secondorder accurate in space and time. The numerical diffusion that usually accompanies a loworder gridbased 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 (27) 3' 18' 36 (27) In these simulations we use a 3parameter collision operator, allowing for separate relaxation of the 5 shear modes, 1 bulk mode, and 9 kinetic modes. The postcollision distribution nt = nt + A, is written as nt (= ac/ C 2 (28) 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 nonlinear puu term in Eq. 28 should be retained, since it maintains Galilean invariance and prevents an artificial crossstream drift [21]. Empirically, we have found that the domain of validity of Stokes flow is considerably larger with the nonlinear term in place. The nonequilibrium momentum flux fIneq neqC C. relaxes due to collisions at the lattice nodes, Fneq,* (1 + A) f + (1 + A,)(fI" : 1)1, (29) 3 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 ( + (210) 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 latticeBoltzmann 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. (212) 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. (213) The macrodynamical behavior arising from the latticeBoltzmann equation can be found from a multiscale 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 NavierStokes equations with corrections that are of the order u2 and c2. Thus, at sufficiently low Mach numbers, the method is secondorder accurate in space with relative errors proportional to Ax2. It is also secondorder accurate in time if the viscosities are defined according to Eq. 210. 2.3 Fluctuations The latticeBoltzmann model can be extended to include thermal fluctuations, which leads to Brownian motion of colloidal particles and polymers. The fluctuating latticeBoltzmann model [3, 4] incorporates a random component into the momentum flux during the collision process (c.f. Eq. 29): neq, (1 + A)eq+q + }(1 + )(fleq 1)1 +If 3 ni ( + (R 1, (214) 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 offdiagonal 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 secondorder accuracy in time. An analogy with fluctuating hydrodynamics [20] sI 1 that ( 2T 1/2 ( 2T2, 1/2 AXAt, AXAt, (215) 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 latticeBoltzmann fluid are [4, 21] 2TA2 1/2 2TA 1/2 ( 16) There is an additional multiplicative factor of A2 (c.f. Eq. 215), but for A = A, = 1, an exact correspondence with the fluctuationdissipation relation for continuous systems is obtained. However the discrete version (Eq. 216) can also be applied to systems where the viscous eigenvalues are not equal to 1 [4]. This may be advantageous in problems where a timescale separation is required between hydrodynamic and kinetic modes. 2.4 Solid Particles 2.4.1 FiniteSize Spherical Particles In our implementation of the LatticeBoltzmann model, the solid particles are defined by a surface which cuts some of the links between lattice nodes (Fig. 21). 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 solidfluid boundaries can be found in references [2, 4, 16]. Boundary Conditions. Figure 21. 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. 21. 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) (217) Cs where n (r, t) is the postcollision distribution at (r,t) in the direction cb and c cb. The local velocity of the particle surface, Ub =U + QX (rb R) (218) is determined from the particle centerofmass velocity U, angular velocity fl, and centerof mass R; rb r + oCbAt is the location of the boundary node. We use the mean density po in Eq.217 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 latticeBoltzmann 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.217, f(r, t + At) 2n (r, t) 2abpOUb (219) 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 latticeBoltzmann simulations of hydrodynamic interactions are given in Ref. [16]. 2.4.2 Point Particles A polymer chain is represented by a beadspring 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 hardsphere excludedvolume interaction between the monomers was included, with collision diameter b/2. Other excluded volume and intrachain potentials could be used. The equations of motion for the monomers are written in inertial form, d2X, ViIJXN) + F, (2 20) dt2 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. 220 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. 220 directly. The large timescale separation between the dynamics of the polymer and the individual monomers allows time for the hydrodynamic interactions to reach a quasisteady 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 latticeBoltzmann model is just one possible framework for such calculations, which could also be carried out within a fluctuating finitedifference or finiteelement simulation of the fluid. The hydrodynamic interactions between the beads are determined by coupling them individually to the fluctuating latticeBoltzmann 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. (221) The random force F' is introduced to balance the additional dissipation caused by using a frictional coupling instead of a noslip boundary condition on the bead surfaces [11]. Since the fluid satisfies its own fluctuationdissipation relation, F' has a local covariance matrix (Fr(t)Fr(t')) 2ToJ6(t t')1. (222) 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 pointparticle 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 hardsphere 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 X, u(X, t) = (223) 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 (224) 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) (225) (226) 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 + iAt) U(t) + Fo At mn m U(t + At)(1 + At) U(t) + U(t) + [ FAt U(t)At] n rn rn r (227) (228) 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 (229) 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]. 230) te These equations coupled with the forces generated at the solidfluid surface due to the fluctuating stresses, results in a discrete Langevinlike equation for the update of the particle velocity. 2.5 Verification of the Fluctuations The fluctuating LatticeBoltzmann 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 rescaling 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 selfdiffusion in colloidal suspensions. This section consists of the tests we have performed to verify our simulation method. The latticeBoltzmann 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 (23 ) The fluctuations we calculated all corresponded to the right input temperature, according to equation 231 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 22 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 coefficients. 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 23. 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 autocorrelation 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 22. 0 200 400 600 800 1000 t 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 timestep, error accumulates for a couple of timesteps 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 24 demonstrates this for our simulation technique and is yet another mean of verification code. As illustrated in Figures 24 and 25, 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 7000 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 t Figure 23. 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 25 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 GreenKubo relationship is free of this artifact and gives the correct StokesEinstein value (equation 2 32) as verified by simulations. These calculations are summarized in Appendix Figure 24. 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. A. These preliminary tests show that the fluctuating latticeBoltzmann scheme is capable of producing the correct transport coefficients for spherical Brownian particles of D = BT (232) 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 explicit S implicit  0.25 t/(a vV) Figure 25. 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 1 A 0.5 ''05 V 0.25 t (a v) CHAPTER 3 DYNAMICS OF POLYMER SOLUTIONS IN PERIODIC AND CONFINED GEOMETRIES 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 latticeBoltzmann model of the solvent [34] with a pointparticle 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 threedimensional 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, endtoend 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 coarsegrained 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 excludedvolume 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, 3840], which is not ..v ,i, sufficient to determine accurate theological properties of flexible polymers [41]. In a semidilute 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. Noslip boundaries require a complicated Green function even for channel flows [36], while more involved geometries can only be handled by finiteelement or boundaryelement 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 lOstz 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 latticeBoltzmann 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 latticeBoltzmann model. We also show that the radius of gyration, endtoend distance, and diffusion coefficient of a flexible excluded volume chain agree quantitatively with renormalization group calculations [42] for chains up to 1024 monomers. A coarsegraining 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 twodimensional 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 finitesize 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 FluctuationDissipation The fluctuationdissipation 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. 216) and o is the monomer friction. However, as Ahlrichs and Diinweg[ll] originally indicated, the diffusion coefficient of an isolated point particle is not 0.5 0.4 S0.3 0.2 0.1 0 Figure 31. O Dissipative * Fluctuating 0 0.2 0.4 0.6 frAx/0 Relation between the effective friction coefficient and the input friction o (Eq. 221). 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 StokesEinstein 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. 31 summarizes the most important finding, namely that the frictioncoefficient determined from the meansquare displacement of the monomer at finite temperatures, 1 = Dm/T, matches the friction coefficient determined from the drag force, Fd = U. Figure 31 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. 221). 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, (32) that is singular at the origin. However, the discrete nature of the latticeBoltzmann velocity field leads instead to a finite result, u(0) Fg (33) where g is a numerical constant. A steadystate force balance between the imposed force F and the particle velocity, F Fd o [U u(0)] = U, (34) leads directly to Eq. 31. Further tests verify that the parameter g is independent of solvent viscosity, confirming the relation given in Eq. 31. However, g does depend weakly on system size as shown in Fig. 32, 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. 33. 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 pointparticle 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 fluctuationdissipation 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 0.04 0.03 Figure 32. 00 0 ao=0.32 a =0.15 0.021 1 1 1 7 0 20 40 60 80 L/Ax Offset between 1 and iO1 as a function of system volume. The variation in the dimensionless mobility g (Eq. 31) 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. 33 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 , (35) 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 time, U(t + At) U(t) QU(t) + AU(t), 0 (36)  Implicit  Explicit 4 r 41r I I I I I 0 0.05 0.1 0.15 0.2 At/m Figure 33. 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). or U(t + At) U(t) QU(t + At) + AU(t), (37) where Q = (At/m and ((AU)2) 2T//m. Equation 36 corresponds to a firstorder explicit integration of Eq. 35, while Eq. 37 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 StokesEinstein relation D = T/l, regardless of the integration method, but classical equipartition of energy is modified, T mKU2 li /2' (3 8) explicit integration leads to the minus sign in the denominator of Eq. 38, while implicit integration leads to the plus sign. Figure 33 shows that this simple model accounts 1.1 S1.05 A V 1 0.95 0.9 for the kinetic energy measured in the simulations, although here the fluctuations in particle velocity are largely driven by the fluid. The longtime 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 flowfield around a point force, as can be seen in Fig. 34 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. (39) kO0 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 (310) and also from the autocorrelation function of the particle velocity in a fluctuating fluid, Ai (U<(t)U (0)) dt. (311) T Jo 0.01 U 1u ^ 01 U . Oseen Stokes 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 ..ee " r/Ax Figure 34. 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 latticeBoltzmann 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 nonzero. The singleparticle 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. 35 and 36. 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. 39, down to separations of the order of Ax. The discrepancies largely reflect errors in interpolating the flow field to offlattice 0.6 A 0 Dissipative D Fluctuating 0.5 A Stokes 0.4 0.3 O A 0.2 6 0.1 0 I I 0 2 4 6 8 r/Ax Figure 35. Comparison of the parallel component of the two particle mobility tensor, /12. Results obtained by dissipative (Eq. 310, circles) and fluctuating (Eq. 311, squares) simulations are compared with an independent solution of the periodic Stokes equations (Eq. 39, triangles). locations; the results shown in Fig. 34 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 latticeBoltzmann model, Ax. Consequently, three different values of b were used to calculate the equilibrium and nonequilibrium 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.4 0 4"0.3 0.2 0.1 0 0 2 4 6 8 r/Ax Figure 36. Comparison of the perpendicular component of the two particle mobility tensor, p12. Results obtained by dissipative (Eq. 310, circles) and fluctuating (Eq. 311, squares) simulations are compared with an independent solution of the periodic Stokes equations (Eq. 39, 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 endtoend distance Re are shown in Fig. 37. 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 selfavoiding random walk [42]. 102 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 N Figure 37. Endtoend 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 latticeBoltzmann 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. 38. The raw diffusion coefficients, shown for two different box lengths, L ~ 3.5Rg and L ~ 7Rg, scale roughly as NA, but the diffusion coefficients increase with the ratio L/Rg. The finitesize effects can be corrected by assuming that the polymer behaves hydrodynamically as a rigid sphere of radius, Rhi, where D = T/6lyRh. (312) O 0 D 0 DL, L7R o DL, L3.5R A D,L7R * D, L3.5R g s 0.59 N 0 5 q 0 o 0] o " 0 0 Figure 38. 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. 312), a selfconsistent estimate of the diffusion coefficient in an infinite system can be obtained. Figure 38 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 selfavoiding random walk [42], D = 0.2T/IbN". The monomer friction in these 1 101 102 102 256 Sv 0 059 S~N 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 102 102 I  16 64 256 1024 N Figure 39. Center of mass diffusion coefficient, with finitesize 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 explicitly. 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. 38 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. 39. 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 polymer. 3.2.4 Polymer Diffusion in a Channel Unlike Brownian dynamics, where external boundaries present a serious complication, the latticeBoltzmann model can simulate the dynamics of confined polymers at no additional cost. Here we consider a channel bounded in the ydirection by planar walls separated by a distance H, with periodic boundary conditions in the x and zdirections; 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 A A p .1 D, L 1H SDII L 2H O DI L 4H AD L 1H * D, L 2H * D, L=4H Figure 310. H/Ax 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 40120 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. 310 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 meansquare displacement, after sufficient time has elapsed for the development of the hydrodynamic interactions, but before a significant displacement of the monomers has occurred. The noslip boundaries screen the hydrodynamic interactions in the direction parallel to the confining walls over 2.5 2 1.5 01 0.5 0 O g O0 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 inplane 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 channel. The scaling of the parallel diffusion coefficient in a confined geometry is shown in Fig. 311 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 powerlaw 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 meanfield calculations [47], which indicate that the .,i,! Ii I ic scaling is not reached in twodimensional 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 meanfield calculation would be useful to confirm that scaling theory applies to much shorter chains in onedimensional 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 latticeBoltzmann model of a thermally fluctuating fluid and a pointparticle 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 twoparticle systems, and found that the longrange Oseen level hydrodynamic interactions were reproduced quantitatively. Ar^A 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.1 0.01 0.1 1 10 R /H g Figure 311. 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, multichain simulations are straightforward and require minimal additional computation if the system volume remains fixed. CHAPTER 4 FLOWINDUCED MIGRATION OF POLYMERS IN DILUTE SOLUTION 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 latticeBoltzmann 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 centerofmass 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, [4850]. In a pressuredriven 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 pressuredriven 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 1vr 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 noslip 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 pressuredriven 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 latticeBoltzmann 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 particlefluid 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 latticeBoltzmann 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. 4la, 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. 42a), 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 multibead 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. 4la) 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. 43a), 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  ... Pe100 /..  0.2  P e 200 g .......... ." 0.3 C015 ./ ...20. 0 0.15 .2 0.22 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 0.5 1 1.5 2 0 0.5 1 1.5 y/R y/R Y g Y g Figure 41. 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 origin. 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. 4la 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. 4la and 43a. 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 42. 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. 41b). The extension data (Fig. 43b) 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. 41c). 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. 42b), 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 centerofmass 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 semiquantitatively 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 (b) 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 43. 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 overpredicts 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 crossover to Rouselike scaling in the range 8 > H/Rg > 4, which supports the screening hypothesis. Moreover, computer simulations [39, 40] which neglect intrabead 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 A A 0.2 CA S0.5 1 1.5 y/R Y g Figure 44. 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. 43c 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.35 1.3 HR 8 A[s H/R=8 0HIR=12 1.25 9g A R =16 v 1.2 Fit 1.15 1.1 0 50 100 150 200 250 300 Pe Figure 45. 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. 44 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 45 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. 46 (solid symbols). The TaylorAris calculation of the dispersion coefficient is an .1i', i.d I ic solution of the convectiondiffusion equation, assuming a uniform distribution of particles across the channel, D* + 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 TaylorAris 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. 46 in good agreement with numerical simulations at the same confinement and at low Peclet numbers. We have 10000 q * q 1000 100 Figure 46. 1000 10000 Hydrodynamic dispersion coefficients of a polymer chain in Poiseuille flow. Numerical simulations are shown as filled symbols, along with TaylorAris 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 liv.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. 47. We also show the effect of inertia on the center of mass distribution in Fig. 48, 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 U U op B * ,'" Simulation, H/R 8 S TA, LdR HR 8  TA, Ld0O 10000 U kT=0.1 kT=0.05 1000 100 1000 Pe P Figure 47. 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. 46, 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 UoH/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 ~ 102 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  0.5 1 2 3 4 y/R Y g Figure 48. 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 intrabead 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 [5557]. 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. CHAPTER 5 TRANSVERSE MIGRATION OF POLYMERS IN EXTERNAL FIELDS 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 pressuredriven 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 pressuredriven 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 noslip 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 pressuredriven 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 [6466]. 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 bodyforce alone, while in combination with a pressuredriven 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 pressuredriven 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 noslip 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 freesolution 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 latticeBoltzmann 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 y/H Figure 51. 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. 51), 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 pvq an important role in limiting the migration. A polymer in the vicinity of Pe = oPe = 1100 4 Pe = 110 OOPe = 11 xxPe = 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. 51). 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 P 2P OP PkT U 9lnJ aln4 + 2 < y > P< D >, (51) 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 r Figure 52. 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 beads. 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 timeindependent solution for the center of mass distribution function, P(y), 0 [A + B P2 + CPF kT aP 6y y2 y2 2( 6Oy (52) 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. 52 describes the lift created for the center of mass of a dumbbell by the .,iv 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 pressuredriven 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. 53A), 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. 53B). 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. 53C). 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 53a 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 53b shows that centerofmass 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 present. 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. 52). The numerical results in Fig. 54 verify this prediction, as can be seen by comparing the results (Pe = 110) in Fig. 51 with Fig. 54 (Pef = 12.5). Since the flow Figure 53. "0 0.1 0.2 0.3 0.4 0.5 y/H 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. 55 show migration towards the boundaries when Figure 54. y/H Center of mass distributions for concurrent application of an external body force and pressuredriven flow. The solid curve shows the level of migration under the pressuredriven 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 pressuredriven 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. 1.5 .Pe=0 SPePe = 22 0.5 / Pe = 44 0Pe = 66 =Pe = 88 AaPe = 110 0.1 0.2 0.3 0.4 0.5 y/H Figure 55. Center of mass distributions for countercurrent application of an external body force and pressuredriven flow. The solid curve shows the level of migration under the pressuredriven 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 shearrate and the applied force (Eq. 52), 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 angleaveraged 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 ADNA 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 103 cm/sec, in comparison with electrophoretic velocities of the order of 101 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 forceinduced 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 forceinduced migration. The magnitude of the applied force could be increased by using lowfrequency AC fields; since the forceinduced 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 forceinduced migration or not. Screening could also be investigated more directly by numerical simulations that include electrostatic forces, either through explicit charges or through a PoissonBoltzmann 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 rt': 1. 1 additional laboratory and numerical experiments to further investigate this question. CHAPTER 6 CONCLUSIONS 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 endtoend 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 .iviiiii. 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. APPENDIX A DISCRETE LANGEVIN EQUATION 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. 33). Nevertheless, the fluctuationdissipation relation between the effective friction and diffusion coefficients holds within statistical errors (Fig. 31). It has been shown[21] that fluctuations in a discrete system can be different from those in continuous systems (c.f. Eqs. 215 and 216), and in order to understand why deviations from classical equipartition do not affect the longtime dynamics, we here examine a simple model with constant friction. The dynamics are described by the classical Langevin equation, dU 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 StokesEinstein relation D = T/l and equipartition m (U2) = T. The real situation is more complex [20, 33], with a timedependent 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), (A2) 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'. (A3) m Jt After n time steps, Eq. A2 has the solution n1 U(t + nAt) = (1 Q)U(t) + AU(t + kAt)(1 Q)(nlk), (A4) k0 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, (A5) tO 1 Q/2' k0 and violates the classical equipartition of energy. On the other hand, the diffusion coefficient, which can be calculated from the discrete equivalent of the GreenKubo relation SAt T D = (U(0)U(0)) + At (U(nAt) U(0)) (A6) nl 1 satisfies the StokesEinstein relation. Thus we recover the correct longtime 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. A1, U(t + At) U(t) QU(t + At) + AU(t), (A7) is often used to circumvent the stability constraint Q < 2 on the explicit solution (Eq. A2). It is straightforward to show that the StokesEinstein relation D = T/l is again satisfied, but in this case the mean kinetic energy is given by M KU2) =T (A8) The comparison of Eqs. A5 and A8 with the numerical simulations is shown in Fig. 33. It can be seen that the simplified model presented here, which neglects the timedependence of the hydrodynamic interactions, still predicts the deviations from equipartition quantitatively. APPENDIX B THE MOMENTUM CHANGE CORRELATION ON A SPHERE 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 (B1) SAt 2At J (1 ) < F(t)F(0) > dt (B2) 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) = ewtCp(w)dw, (B3) 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 (B5) 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! (B6) 1Ww2A S. + O(w3). (B7) 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)] (B10) 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 ] (B12) Substituting this function we find, < APAP > :;:.,/. TAt(R + R2) (B13) TvAtR APPENDIX C THE MOMENTUM CHANGE CORRELATION ON A PLANE WALL, CONTINUUM THEORY 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., rStokes 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) Uoevy v^ \ (Ct) It now remains to calculate the shear stress aux at the wall(y = 0). u(y =0) PW= (1 + )wt. (C2) The friction is the part that is proportional to u, w)= (1+ i) (C3) Thus the friction kernel per unit area, the power spectrum of the forces, can be written as CF(w) kT [(w) + 7(w)] (C4) CF(w) = kT 2pw. (C5) 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 TAti (C6) Again it's important to note this result should be integrated once more over the surface of interest. APPENDIX D THE MOMENTUM CHANGE CORRELATION ON PLANE WALL, DISCRETE THEORY The momentum correlation on a plane wall can also be calculated from the information on the population densities and the use of bounceback 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 (D1) APy 2(n15 17)At (D2) 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 j.ci [neq'" (cici C'2) ni* ac[p ci + c ] (D4) S S 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 + ] + + + ] (D5) 4C 4 12 24 ns = a [p + ] + [ + ] (D6) 1 4 t2 24 These give < AP>2 > < > + < n2 > (D7) Since < j .j >= pAx3kT [2] and < j* > is a third of this, and using the equation 231, and that p = 36 in our model, we find. 4 < AP2 > kT + 12kT (D8) This is the momentum change correlation in x direction per unit area, and the same expression is found for the z direction as well. REFERENCES [1] C. W. J. Beenakker and P. Mazur, "Selfdiffusion of spheres in a concentrated suspension," Ph;,:.i. A, vol. 120, pp. 388410, 1983. [2] A. J. C. Ladd and R. Verberg, "LatticeBoltzmann simulations of particlefluid suspensions," J. Stat. Phys., vol. 56, pp. 286311, 2001. [3] A. J. C. Ladd, "Shorttime motion of colloidal particles: Numerical simulation via a fluctuating latticeBoltzmann 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, LowRB ,.,. ..11 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 latticegas 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 ysphere hydrodymamic interactions and mobilities in a suspension," P,;;,'... A, vol. 115, pp. 2157, 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. 82258239, 1999. [12] B. J. Alder and T. E. Wainwright, "Decay of the velocity autocorrelation function," Phys. Rev. A, vol. 1, no. 1, pp. 1821, 1970. [13] T. Phung, J. Brady, and G. Bossis, "Stokesian dynamcis simulation of brownian suspensions," J. Fluid Mech., vol. 313, pp. 181207, 1996. [14] J. P. Boon and S. Yip, Molecular H. ,h..1 /,., ,i,:. McGraw 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, 1994. [16] N. Q. Nguyen and A. J. C. Ladd, "Lubrication corrections for latticeboltzmann simulations of particle suspensions," Phys. Rev. E, vol. 66, no. 046708, 2002. [17] 0. B. Usta, A. J. C. Ladd, and J. E. Butler, "LatticeBoltzmann simulations of the dynamics of polymer solutions in periodic and confined geometries," J. Ch. im Phys., vol. 122, pp. 0949021, 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, "LatticeBoltzmann simulations of particlefluid suspensions," J. Stat. Phys., vol. 104, pp. 11911251, 2001. [22] A. J. C. Ladd, H. Gang, J. X. Zhu, and D. A. Weitz, "Timedependent collective diffusion of colloidal particles," Phys. Rev. Lett., vol. 74, pp. 318, 1995. [23] A. J. C. Ladd, "Hydrodynamic screening in sedimenting suspensions of nonBrownian 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. 287311, 1998. [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. 705108, 1996. [27] N. S. Martys and H. C'!. i,1 "Simulation of multicomponent fluids in complex threedimensional geometries by the lattice boltzmann method," Phys. Rev. E, vol. 53, no. 1, pp. 743750, Jan 1996. [28] L. W. P. Manz, B; Gladden, "Flow and dispersion in porous media: Latticeboltzmann and nmr studies," Aiche J., vol. 45, pp. 18451854, 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. PhanThien, N. T. Yong, X. Wu, and D. Xu, \!icrochannel flow of a macromolecular suspension," Phys. Fluids, vol. 15, no. 1, pp. 1121, 2003. [33] E. H. Hauge and A. MartinLoef, "Fluctuating hydrodynamics and brownian motion," Journal of Statistical Ph'. vol. 7, no. 3, pp. 259281, 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. 11651173, 2003. [36] R. M. Jendi. ..1.: D. C. Schwartz, J. J. de Pablo, and M. D. Graham, "Shearinduced migration in flowing polymer solutions: Simulation of longchain DNA in microchannels," J. C(' ,in Phys., vol. 120, pp. 25132529, 2004. [37] M. Fixman, "Construction of the Langevin forces in the simulation of hydrodynamic interaction," Macromolecules, vol. 19, pp. 12041207, 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. 28942900, 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. 281298, 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. 299318, 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. 520528, 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, "Likecharge attraction and hydrodynamic interaction," Phys. Rev. Lett., vol. 85, no. 23, pp. 49764979, 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. 40464052, 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. 12251256, 1954. [49] J. H. Aubert and M. Tirrel, "Effective viscosity of dilute polymer solutions near confining boundaries," J. Ch in, Phys., vol. 77, pp. 553561, 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. 8588, 1986. [51] W. F. Busse, "Two decades of highpolymer physics a survey and forecast," Phys. T ./.i vol. 17, pp. 32, 1964. [52] A. B. Metzner, Y. Cohen, and C. RangelNafaile, "Inhomogeneous flows of nonNewtonian fluids: Generation of spatial concentration gradients," J. Non Newtonian Fluid Mech., vol. 5, pp. 449462, 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. 16931717, 1994. [54] K. A. Dill and B. H. Zimm, "Rheological separator for very large DNA molecules," Nucleic Acids Research, vol. 7, pp. 735749, 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. 955971, 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. 121, 97. [61] R. L. Schiek and E. S. G. Shaqfeh, "Crossstreamline migration of slender Brownian fibres in plane Poiseuille flow," J. Fluid. Mech, vol. 332, pp. 2339, 1997. [62] R. G. Larson, "The rheology of dilute solutions of flexible polymers: Progress and problems," J. Rheol., vol. 49, no. 1, pp. 170, 2005. [63] 0. B. Usta, J. E. Butler, and A. J. C. Ladd, "Flowinduced 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. 813872, 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. 45364547, 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. 2932, 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. 273287, 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. 166, 1994. 