Citation
Effects of elasticity and pressure dependent viscosity on the relative motion of binary sphere interactions in highly viscous fluid flows

Material Information

Title:
Effects of elasticity and pressure dependent viscosity on the relative motion of binary sphere interactions in highly viscous fluid flows
Creator:
Woodman, Adam H. ( author )
Place of Publication:
Denver, Colo.
Publisher:
University of Colorado Denver
Publication Date:
Language:
English
Physical Description:
1 electronic file (92 pages) : ;

Thesis/Dissertation Information

Degree:
Master's ( Master of science)
Degree Grantor:
University of Colorado Denver
Degree Divisions:
Department of Mechanical Engineering, CU Denver
Degree Disciplines:
Mechanical engineering

Subjects

Subjects / Keywords:
Frictional resistance (Hydrodynamics) ( lcsh )
Rheology ( lcsh )
Elasticity ( lcsh )
Viscosity ( lcsh )
Elasticity ( fast )
Frictional resistance (Hydrodynamics) ( fast )
Rheology ( fast )
Viscosity ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Review:
The effect of pressure-dependent fluid viscosity and particle elasticity in isothermal suspensions of spherical particles is investigated using numerical simulation software. The relative motion of two interacting particles suspended in viscous shear flows can be decoupled into four modes of motion. These modes include translation along their line of centers, translation perpendicular to their line of centers, rotation about their line of centers, and rotation about an axis perpendicular to their line of centers. Exact solutions to Stoke's equations exist for all four modes but are limited to rigid particles in a fluid with constant viscosity throughout the flow field. In reality, these particles are deformable and the viscosity of the surrounding fluid may increase significantly in high pressure regions.
Review:
In particular, high pressure regions can arise in the interstitial region between the particles when they are sufficiently close. In this research, viscous forces are assumed to dominate the inertial forces so that Stoke's equations apply. Using finite element software, the affect of pressure-dependent viscosity and particle elasticity is studied for each mode of motion at a range of particle separation distances. This data allows the quantification of how pressure-dependent viscosity and particle elasticity affect the hydrodynamic forces on binary particles.
Review:
Introducing pressure-dependent viscosity substantially increased the force on the particles when they translate along their line of centers, compared to a constant viscosity fluid. Pressure increases were not sufficient for the viscosity to increase significantly in all other modes of motion, and the numerical pressure-dependent and constant viscosity simulations produced essentially the same results. The effect of pressure-dependent viscosity is most prevalent when two particles are translating along their line of centers because the highest pressures are generated in this motion. The fluid viscosity increased significantly in this motion in the interstitial region between the particles, leading to a large increase in the force. When binary particles rotate about an axis perpendicular to their line of centers, translate about an axis perpendicular to their line of centers, or rotate about the line of centers, the fluid pressure in the interstitial region is insufficient to increase the forces or torques on the particles significantly.
Review:
Particle elasticity was found to have an insignificant affect on the hydrodynamic forces on the particles. The actual surface separation distance between particles is the determining factor for the fluid pressure and forces on the particles. While the particles deform in the case of particles translating along their line of centers, the deformation causes the surface separation distance to increase. The center-to-center distance of the particles has not changed, so it appears that the fluid pressure and hydrodynamic forces decrease more for softer particles. When the hydrodynamic forces are compared to the actual surface separation, the forces are found to be approximately equal to the forces on rigid particles. Comparing different values of the particles' Poisson's ratio, a similar affect to elasticity on the forces was observed.
Bibliography:
Includes bibliographical references.
System Details:
System requirements: Adobe Reader.
Statement of Responsibility:
by Adam H. Woodman.

Record Information

Source Institution:
University of Colorado Denver Collections
Holding Location:
Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
987252043 ( OCLC )
ocn987252043
Classification:
LD1193.E55 2016m W66 ( lcc )

Downloads

This item has the following downloads:


Full Text
EFFECTS OF ELASTICITY AND PRESSURE DEPENDENT VISCOSITY ON THE
RELATIVE MOTION OF BINARY SPHERE INTERACTIONS IN HIGHLY
VISCOUS FLUID FLOWS
by
ADAM H WOODMAN B.S., University of New Hampshire, 2013
A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Mechanical Engineering
2016


This thesis for the Master of Science degree by Adam H Woodman has been approved for the Department of Mechanical Engineering by
Marc Ingber, Chair Maryam Darbeheshti, Advisor Dana Carpenter Alan Graham
December 17, 2016
n


Woodman, Adam H (M.S. Mechanical Engineering)
Effects of Elasticity and Pressure Dependent Viscosity on the Relative Motion of Binary Sphere Interactions in Highly Viscous Fluid Flows
Thesis directed by Assistant Professor Maryam Darbeheshti
ABSTRACT
The effect of pressure-dependent fluid viscosity and particle elasticity in isothermal suspensions of spherical particles is investigated using numerical simulation software. The relative motion of two interacting particles suspended in viscous shear flows can be decoupled into four modes of motion. These modes include translation along their line of centers, translation perpendicular to their line of centers, rotation about their line of centers, and rotation about an axis perpendicular to their line of centers. Exact solutions to Stokes equations exist for all four modes but are limited to rigid particles in a fluid with constant viscosity throughout the flow field. In reality, these particles are deformable and the viscosity of the surrounding fluid may increase significantly in high pressure regions.
In particular, high pressure regions can arise in the interstitial region between the particles when they are sufficiently close. In this research, viscous forces are assumed to dominate the inertial forces so that Stokes equations apply. Using finite element software, the affect of pressure-dependent viscosity and particle elasticity is studied for each mode of motion at a range of particle separation distances. This data allows the quantification of how pres sure-dependent viscosity and particle elasticity affect the hydrodynamic forces on binary particles.
Introducing pressure-dependent viscosity substantially increased the force on the particles when they translate along their line of centers, compared to a constant viscosity
fluid. Pressure increases were not sufficient for the viscosity to increase significantly in
iii


all other modes of motion, and the numerical pressure-dependent and constant viscosity simulations produced essentially the same results. The effect of pressure-dependent viscosity is most prevalent when two particles are translating along their line of centers because the highest pressures are generated in this motion. The fluid viscosity increased significantly in this motion in the interstitial region between the particles, leading to a large increase in the force. When binary particles rotate about an axis perpendicular to their line of centers, translate about an axis perpendicular to their line of centers, or rotate about the line of centers, the fluid pressure in the interstitial region is insufficient to increase the forces or torques on the particles significantly.
Particle elasticity was found to have an insignificant affect on the hydrodynamic forces on the particles. The actual surface separation distance between particles is the determining factor for the fluid pressure and forces on the particles. While the particles deform in the case of particles translating along their line of centers, the deformation causes the surface separation distance to increase. The center-to-center distance of the particles has not changed, so it appears that the fluid pressure and hydrodynamic forces decrease more for softer particles. When the hydrodynamic forces are compared to the actual surface separation, the forces are found to be approximately equal to the forces on rigid particles. Comparing different values of the particles Poissons ratio, a similar affect to elasticity on the forces was observed.
The form and content of this abstract are approved. I recommend its publication.
Approved: Maryam Darbeheshti
IV


DEDICATION
This thesis is dedicated to my parents, who have supported everything I have pursued. They have taught me what hard work is, allowed me to find my interests, and have been my main source of encouragement and inspiration throughout my life. Nothing I have accomplished would have been possible without their endless advice, love, support, and dedication to being the best role models possible.
v


ACKNOWLEDGMENTS
This research would not have been possible without the opportunities that my advisor, Dr. Maryam Darbeheshti, has presented me with. Her knowledge, support, kindness, and positive attitude were key factors in my graduate career.
I also want to thank Dr. Ingber and Dr. Graham for their introduction to this research. Their endless knowledge on the topic made it possible for me to continue after hitting road blocks. They asked the right questions and pointed me in the direction of several sources of information that were essential to this thesis.
vi


TABLE OF CONTENTS
CHAPTER
I. INTRODUCTION......................................................... 1
1.1 Governing Equations ........................................... 2
1.2 Binary Particle Interactions................................... 4
1.3 Literature Review.............................................. 4
II. ANALYTIC SOLUTIONS.................................................. 9
2.1 Related Benchmark Solutions.................................... 9
2.2 Solutions to Four Modes of Motion............................. 10
III. NUMERICAL SIMULATIONS............................................. 16
3.1 Mesh Discretization........................................... 21
3.2 Boundary Conditions........................................... 23
IV. FLUID CONSTITUTIVE EQUATIONS...................................... 25
4.1 Viscosity Models.............................................. 25
V. RESULTS: PRESSURE DEPENDENT VISCOSITY............................. 28
5.1 Rotation About Axis Perpendicular to Line of Centers.......... 28
5.2 Translation Along Axis Perpendicular to Line of Centers... 34
5.3 Rotation About Line of Centers ............................... 39
5.4 Translation Along Line of Centers............................. 43
5.5 Combined Motions.............................................. 48
VI. RESULTS: PARTICLE ELASTICITY...................................... 53
6.1 Rotation About Axis Perpendicular to Line of Centers.......... 53
6.2 Translation Along Axis Perpendicular to Line of Centers... 57
6.3 Rotation About Line of Centers ............................... 58
6.4 Translation Along Line of Centers............................. 59
6.5 Poissons Ratio............................................... 65
VII. CONCLUSION........................................................ 70
vii


REFERENCES
72
APPENDIX
A. Mesh Convergence Studies.................................... 74
vm


FIGURES
FIGURE
1.1 Four modes of motion possible for binary particle interactions............. 4
2.1 Torque on binary spheres, one rotating about an axis perpendicular to their
line of centers and one stationary, normalized by the torque on an isolated rotating sphere.......................................................... 11
2.2 Force on binary spheres, one translating along axis perpendicular to line of
centers and one stationary, normalized by F'stokes.......................... 13
2.3 Torque on binary spheres, one rotating about line of centers and one stationary, normalized by the torque on an isolated rotating sphere................. 14
2.4 Force on binary spheres translating along their line of centers, normalized by
FStokes..................................................................... 13
3.1 Comparison of the analytic and numerical solutions for a rigid sphere approaching an infinite plane wall............................................. 19
3.2 Example of mesh for binary spheres traveling along their line of centers. . . 22
3.3 Detailed view of interstitial region mesh for binary spheres.............. 22
4.1 Comparison between Barus and Roelands models for viscosity of mineral oil. 26
4.2 Roelands pressure-viscosity model for four liquids....................... 27
5.1 Geometry used in mode 1 and 2 simulations................................. 29
5.2 Detailed view of the geometry surrounding the spheres for modes 1 and 2. . 29
5.3 Overview of the mesh used for modes 1 and 2............................... 30
5.4 Detailed view of the mesh used for modes 1 and 2. A high mesh density region is defined close to the interstitial region and the mesh grows coarser toward
the bulk fluid region....................................................... 31
5.5 Comparison of n = flip) and n = const, numerical results and n = const, analytic solution for mode 1. Fluid is mineral oil, particle velocity is Q = 1^. 33
5.6 Comparison of n = nip) and n = const, numerical results and n = const,
analytic solution for mode 2. Fluid is mineral oil, particle velocity is U = 1. 35
IX


5.7 Velocity vectors on the cross section of binary spheres for modes 1 and 2,
from left to right........................................................ 36
5.8 Mode 1 (left) and 2 (right) diagram of pressure distribution. LP denotes
low pressure, while HP denotes high pressure............................ 37
5.9 Pressure distribution on the surfaces of mode 1 spheres near the interstitial
region...................................................................... 38
5.10 Pressure distribution on the surfaces of mode 2 spheres near the interstitial
region...................................................................... 39
5.11 Mesh scheme for modes 3 and 4, with three separate regions of different
element densities........................................................... 40
5.12 Comparison of n = flip) and n = const, numerical results and n = const,
analytic solution for mode 3. The fluid is mineral oil and the rotating particle angular velocity is Q = 1%^........................................... 42
5.13 Comparison of n = nip) and n = const, numerical results and n = const, analytic solution for mode 4. Fluid is mineral oil, particle velocity is U = 2.15
m 44
5.14 Mode 4 approach velocity vs. normalized drag force for pressure-dependent
viscosity mineral oil at a separation distance of e = ^|.................. 47
5.15 Numerical results for mode 4 approach velocity vs. viscosity and pressure for
pressure-dependent viscosity mineral oil at a separation distance of e = ^|. 48
5.16 Individual motions of modes 1-3 combined with mode 4 motion............... 49
5.17 Normalized force or torque for combination of modes 1 and 2 with mode 4.
Mode 1 is normalized by the torque on an isolated rotating particle and mode
2 is normalized by the force on a particle in Stokes flow................... 52
6.1 Mesh used for the mode 1 solid mechanics module to compute the deformations of the spheres, with 143676 elements and 600459 degrees of freedom. . 54
6.2 Exaggerated deformed shape and deformation magnitudes of elastic spheres
in mode 1 motion, e = E = 0.2GPa, v = 0.4, scaled to 10% of the total geometry.............................................................. 55
6.3 Exaggerated deformed shape of elastic spheres in mode 2 motion, e =
E = 0.2GPa, v = 0.4, scaled to 10% of the total geometry. Deformations are normalized by sphere radius........................................... 57
6.4 Mesh used for a fully-coupled simulation of mode 4 with elastic spheres. . 60
x


6.5 Exaggerated deformation of mode 4 elastic particles of E=0.2 GPa, e =
in a highly viscous fluid with n = 0.5, normalized by sphere radius..... 61
6.6 Maximum fluid pressure in the interstitial region for mode 4 motion of elastic
spheres in mineral oil, at a range of particle separation distances. E = 20 GPa, v = 0.4........................................................ 62
6.7 Mode 4 hydrodynamic drag forces with elastic particles, where the separation
distance does not account for surface deformation....................... 64
6.8 Mode 4 hydrodynamic drag forces with elastic particles. Separation distance
has been adjusted to account for the surface deformation................ 64
6.9 Comparison of Poissons ratios vs. hydrodynamic force for E = 0.2 GPa,
< ................................................................. ee
6.10 Comparison of Poissons ratios of 0.25 and 0.49 for several elasticities. Left: deformation unaccounted for in separation, right: deformation accounted for. 68
6.11 Comparison of v = 0.25 and v = 0.49 for E = 0.2 GPa. Left: deformation
unaccounted for in separation, right: deformation accounted for......... 69
xi


TABLES
TABLE
3.1 Mesh convergence study results for the drag force on an isolated sphere in
Stokes flow..........................................................
3.2 Relative error of numerical simulations for a sphere approaching a plane wall,
compared to the analytic solution.....................................
4.1 Viscous properties of four different liquids...........................
5.1 Details of the numerical mesh scheme used for mode 1 simulations at each
separation distance...................................................
5.2 Analytic and numerical solutions for mode 1 rotating sphere in fj, = const.
mineral oil, normalized by the torque on an isolated sphere...........
5.3 Numerical solutions for fj, = const, and fi = fi{p) numerical solution for mode
1 rotating sphere in mineral oil, normalized by the torque on an isolated sphere.
5.4 Comparison of maximum gauge pressure and viscosity in constant and pressure dependent viscosity mineral oil, for mode 1 particle motion with Qrot =
1 f and < = .....................................................
5.5 Comparison of analytic and numerical solutions for mode 2 moving sphere in
ji = const, mineral oil, normalized by the force on a single sphere in Stokes flow..................................................................
5.6 Comparison of fj, = const, and fj, = fi{p) numerical force on mode 2 moving
sphere, normalized by the force on a single sphere in Stokes flow.....
5.7 Comparison of maximum gauge pressure and viscosity in constant and pressure-dependent viscosity mineral oil, for mode 2 motion with [7=1 and e = ra^'s.
5.8 Details of the numerical mesh scheme used for mode 3 simulations at each
separation distance...................................................
5.9 Comparison of analytical and numerical solutions for mode 3 rotating sphere in mineral oil with fj, = const., normalized by the torque on an isolated sphere.
5.10 Comparison of fj, = const, and fj, = fi{p) numerical torque on mode 3 rotating sphere in mineral oil, normalized by the torque on an isolated rotating sphere.
5.11 Comparison of maximum gauge pressure in constant and pressure-dependent
viscosity mineral oil, for mode 3 particle motion with and e = ra^'s.
18
20
26
31
32
32
33
34
35
36
41
41
42
43
xii


5.12 Details of the numerical mesh used for mode 4 simulations at each separation
distance....................................................................... 44
5.13 Mode 4 force, max. pressure, and max. viscosity for n = const, and n = flip) mineral oil at multiple approach velocities and a separation of e = . ... 45
5.14 Increases in mode 4 pressure, viscosity, and force from n = const, to n = nip)
mineral oil, at e = ra^'s and a range of approach velocities................... 46
5.15 Torque on mode 1 rotating sphere at Qi = 1^ when combined with mode 4
motion at a range of approach velocities....................................... 50
5.16 Force on mode 2 moving sphere at U2 = lgf when combined with mode 4
motion at a range of approach velocities....................................... 51
5.17 Torque on mode 3 rotating sphere at Q3 = 1^ when combined with mode 4
motion at a range of approach velocities....................................... 51
6.1 Maximum particle deformations for mode 1, corresponding to fluid pressures
obtained from rigid sphere simulation, normalized by sphere radius............. 56
6.2 Maximum particle deformations for mode 2 at several separation distances
and particle elasticities, corresponding to fluid pressures obtained from rigid sphere simulation, normalized by sphere radius.............................. 58
6.3 Maximum particle deformations for mode 3 at several separation distances
and particle elasticities, corresponding to fluid pressures obtained from rigid sphere simulation, normalized by sphere radius........................... 59
6.4 Maximum normalized mode 4 particle deformations at several separation distances and particle elasticities in a fluid with n = 0.5 Pa s. The simulation
is fully-coupled between solid mechanics and creeping flow..................... 61
6.5 Normalized forces on mode 4 elastic binary particles at a range of separation
distances and particle elasticities.............................................. 63
6.6 Percent decrease in mode 4 force on elastic particles, compared to rigid particles. 63
6.7 Normalized hydrodynamic forces at e = and a range of elasticities and
Poissons ratios................................................................ 65
6.8 Percent decrease in mode 4 force on elastic particles, compared to rigid particles. 67
6.9 Particle deformations at e = for a range of elasticities and Poissons ratios. 67
xiii


A.l Convergence study for mode 1 three-dimensional mesh at e = based on
change in force between subsequent meshes............................... 74
A.2 Convergence study for mode 2 three-dimensional mesh at e = based on
change in force between subsequent meshes............................... 75
A.3 Convergence study for mode 3 two-dimensional mesh at e = based on
change in force between subsequent meshes............................... 75
A.4 Convergence study for mode 4 two-dimensional mesh at e = based on
change in force between subsequent meshes............................... 76
A.5 Convergence study for mode 1 three-dimensional solid mechanics mesh at e = and E = 0.2 GPa based on change in deformation between subsequent meshes............................................................... 77
A.6 Convergence study for mode 3 two-dimensional solid mechanics mesh at e =
YP and E = 0.2 GPa based on change in deformation between subsequent meshes............................................................... 77
A.7 Convergence study for mode 4 two-dimensional fully-coupled simulation mesh
at e = and E = 20 GPa based on change in force between subsequent meshes. 78
xiv


CHAPTER I INTRODUCTION
Suspensions of particles in viscous fluids are important in many fields and applications such as sediment transport, chromatography, oil recovery, contaminants, and lubrication theory. The interactions between particles in isothermal, low Reynolds number flows, or Stokes flows, have been studied for a long time because of their applicability to many fields. Particles are often approximated as spheres for simplification purposes. When two spheres are close to one another in shear flows, their relative motion can be described as a combination of four different modes. The spheres can translate along their line of centers, translate perpendicular to their line of centers, rotate about their line of centers, and rotate about the axis perpendicular to their line of centers. Analytic solutions to these modes have been derived and are exact solutions to the Stokes equations. The analytic solutions assume isothermal, low Reynolds number, and incompressible flows.
The analytic solutions are limited to rigid particles, but if there are sufficient fluid pressures and the particles are elastic, they will deform. The viscosity of the fluid in these isothermal solutions is assumed to be constant throughout the flow held, but fluid viscosity is actually dependent on the local pressure. This research will investigate how particle elasticity and pressure-dependent fluid viscosity affect the hydrodynamic forces and torques on particles, assuming the hows are isothermal and incompressible.
The hnite element method presents an attractive analysis tool for problems such as this because of its ability to accurately capture the pressure gradients in the huid. For example, when two spheres are approaching each other, the huid pressure is the highest in the interstitial region between the spheres. We can dehne multiple elements across the interstitial gap, to accurately capture the large gradients in pressure and velocity. COMSOL Multiphysics software was used to solve the problems numerically and allowed for the hne-tuning of mesh discretization, boundary conditions, and solving techniques.
1


1.1 Governing Equations
The cornerstone of computational fluid dynamics is the equations that govern the fluid body. These equations describe the pressure and motion of a body of fluid at each point in the computational domain. The fluid region is governed by Stokes equations, since we assume the viscous forces dominate the inertial forces. This type of flow is known as Stokes flow or creeping flow and is characterized by very low Reynolds numbers. In particular, we assume the Reynolds number, which is a determining characteristic of the flow,
Re oc
inertial forces
(1.1)
viscous forces
is significantly less than 1, which simplifies the solution of the governing equations. The Navier-Stokes equations are commonly accepted as the governing equations for the motion of isothermal fluids. In general, the Navier-Stokes equations can be written as the conservation of momentum equation and the conservation of mass equation, which are solved together. In the case of a compressible Newtonian fluid in an isothermal flow, the momentum equation is written as
P
<9u
iv + u'Vu
(i)
^-Vp + V ^p, ^Vu + (Vu)Tj
(2) '-------------------v
(3)
2
3
fjL (V-u)I
)+v£
(4)
(1.2)
where the vector u represents fluid velocity in the three orthogonal directions, and 1.2 is three separate equations written as one, in vector notation. The continuity, or conservation of mass, equation for a compressible fluid is written as
^ + v (pu) = 0,
(1.3)
where p is the fluid density, p is the fluid pressure, t is time, u is the velocity vector, and
2


/ is the identity matrix. The individual terms in equation 1.2 correspond to the inertial forces (1), pressure forces (2), viscous forces (3), and body forces (4). In the case of the motion of binary spheres in a viscous fluid, there are a few assumptions that can be made to simplify these equations.
We assume the amount of time has been sufficient for the system to reach a steady state and a fully developed flow held. Under this steady state assumption, all terms dependent on time are neglected and
| = . (1.4)
in the momentum equation. Since the viscous forces dominate the inertial forces,
and this term can be neglected in equation 1.2. The body force is given by the gravitational potential which is combined with the pressure forces. If density variations in the fluid are negligible, we have incompressible how and p can be moved outside of the gradient in equation 1.3, resulting in
Substituting equation 1.6 into term (3) in equation 1.2, the incompressible momentum equation is written as
u Vu 0
(1.5)
V u = 0.
(1.6)
(1.7)
which is solved together with equation 1.6.
3


1.2 Binary Particle Interactions
In suspensions of solid particles, the particles have a significant affect on one another. This study will focus on how binary particles interact with one another in a viscous shear flow. Two spherical particles can interact in a combination of four different motions which can be decoupled and solved for independently [18].
Model Mode 2 Mode 3 Mode 4
Figure 1.1: Four modes of motion possible for binary particle interactions.
The motions in figure 1.1 will be identified with their respective label. Mode 1 will indicate a motion where two particles are rotating relative to one another about an axis perpendicular to their line of centers. Mode 2 will denote particles are translating relative to one another along the axis perpendicular to their line of centers. Mode 3 will be the motion where two particles are rotating relative to each other about their line of centers. Mode 4 will describe the motion of two particles translating relative to one another along their line of centers [18]. This convention will be used throughout the thesis.
1.3 Literature Review
The resistance of a spherical particle moving slowly through a viscous fluid was studied by George Stokes in the mid 1800s. Stokes Law was derived in 1851 to solve for the drag force on these particles, at constant velocities [8]. This widely accepted solution is exact only for fluid regions that extend infinitely in all directions. The presence of boundaries in close proximity to particles significantly perturbs the particle resistance
4


formulation and Stokes Law no longer applies [7]. The influence of particles on each other in suspensions became an interest to researchers because of the presence of these types of interactions in many natural and mechanical systems.
In the early 1900s, Einstein found the viscosity of slurries to increase with the volume fraction of solid particles, and developed an equation to describe the relative viscosity of dilute slurries up to a volume fraction of approximately 0.02 [12]. The equation is limited to creeping flows and incompressible Newtonian fluids. Guth and Simha extended Einsteins equation in 1936 for particle volume fractions up to approximately 0.06, by considering the bulk interaction between the solid particles [14]. Thomas [26] developed an improved equation in 1965 using empirical data to extend Einsteins work to volume fractions of approximately 0.1.
The research of dilute suspensions naturally led to the research of models to describe higher concentration systems. For suspension flows with high particle concentrations, Krieger and Dougherty developed a model in 1959 relating the viscosity to the particle packing parameter [20]. In 1971, Chong et al. [10] investigated different ratios of particle sizes and how their variation affects a slurrys bulk viscosity, Ending that the minimum viscosity is achieved with approximately 25% to 35% of the solid concentration being fine spheres, and the remainder being coarse spheres. In 1972, Batchelor and Green [5] developed a model for semi-dilute suspensions, which predicts a viscosity proportional to the square of the volume fraction. In 2012, Sohn and Koo derived an expression to describe the permeability of viscous flow through a packed bed of bidisperse sized spheres, by considering the drag force on each sphere [24],
The research of the bulk viscosity of slurries lead to a rise in popularity of researching unary and binary particle systems. In 1961, Howard Brenner introduced an expression for a sphere slowly approaching an infinite, rigid plane wall in a viscous fluid, as a modification to Stokes Law [7]. In the same time frame, Brenner formulated an expression for a spherical particle approaching a free surface, where the normal component of velocity
5


and tangential stresses at the plane surface vanish, creating a symmetry condition [7]. In 1970, Majumdar and ONeill [21] presented exact solutions to binary spheres rotating about an axis perpendicular to their line of centers and translating along the same axis. These solutions are correct to order In e, where e is the separation distance between spheres. In 1984, Jeffrey & Onishi [19] extended Majumdar & ONeills work, correct to order e In e, and also presented the solution for binary spheres rotating about their line of centers. In 1989, Ascoli, Dandy, and Leal developed a boundary integral formulation to the hydrodynamic interaction of a solid particle with a planar wall, obtaining excellent agreement with Brenners solution [1],
In 2000, Fernandez, Carrica, and Drew [13] calculated the average drag force exerted on a particle, accounting for binary particle interactions. They determined that the drag force depends on particle concentration, gradients of concentration, and differences in sphere sizes. In 2006, Challa & Van Swol [9] reported on the drag force experienced by a sphere approaching a plane wall to test the predictions of hydrodynamic theory and introduce Van der Waals forces to investigate their effect. They found that at low speeds and small separation distances, the presence of solvation forces can lead to a drag force that oscillates between positive and negative. This finding contradicts the classical hydrodynamic prediction of a divergent drag force as separation distances decrease. The result complicates the testing of hydrodynamic theory for binary particles and suspensions of particles, where static forces may be present.
In addition to the forces and couples on each individual sphere, investigation into the diffusion and migration of spheres in suspension systems has been performed. This is an important topic for many mechanical systems used in mineral processing and sedimentation. In 1970, Halow & Wills released their experimental observations of sphere migration in Couette systems and provide a method to calculate sphere trajectories when the sphere and liquid densities are equal [16]. In 1998, Tetlow, Graham, Ingber, and Subia published their findings of particle migration experiments and simulations [25]. They observed an
6


insignificant affect of particle roughness on the migration of particles. In 2008, Ingber, Feng, Graham, and Brenner presented their work on a new traction-corrected boundary element method to solve for the migration of rough spheres in nonlinear shear shows [18]. This new method was found to be very accurate for predicting the trajectories of particles. DAvino, Snijkers, Pasquino, et al. investigated the migration of a single sphere in a viscoelastic liquid in 2012, using Couette flow [11]. They found the sphere to migrate in the direction of decreasing shear rate gradient.
Reynolds established the basis for all lubrication theory problems in 1886, with the Reynolds equation [22], The equation relates the lubricant pressure to the geometry and movement of the machine components. In 1991, Hamrock [17] presented research on fluid film lubrication where two bearing surfaces are completely separated by a thin layer of fluid film. Rheology is the study of the flow and deformation of matter, under applied forces. Bair [2] provided research into linking the fields of rheology and elastohydrodynamic (EHD) lubrication. EHD lubrication problems are those where the pressure in the lubricant is high enough to result in an elastic deformation of the adjacent body. Habchi [15] developed a finite element approach to ultra-low-viscosity fluid EHD full film lubrication problems in his 2008 thesis, using a fully-coupled Newton-Raphson method. The lubricant pressures in EHD lubrication problems are extremely high, significantly affecting the viscosity. Habchi studied several pressure-viscosity and pressure-temperature-viscosity models to employ in his simulations, giving great insight into the models that are suited for different cases. He found that the shape of the fluid film regime, determined by the deformed shape of the body, significantly affects the pressure field and the fluids viscosity. Habchi also determined that the variation of pressure across the film thickness is negligible compared to the variation along the contact plane.
Research in the area of fluid constitutive equations, such as those governing the viscosity of a fluid as a function of pressure, led to many new formulations. Since the pressure of lubricants in machine components often reach values in the gigapascal range,
7


it is important to have models to describe the properties of the lubricants at these extreme pressures. A widely used equation to describe the relationship between viscosity and pressure in a fluid is the Barus equation, developed by Carl Barus in 1893 [3]. Ham-rock [17] points out that many viscosity-pressure relationships, such as the Roelands formula that was developed by Roelands in 1966 [23], were developed after the Barus equation, because it significantly over predicts experimental viscosities at higher pressures. Roelands also developed a model to describe the viscosity of a fluid based on both pressure an temperature. Both of Roelands models are empirically formulated.
In 1983, Happel & Brenner [8] bridged the gap between the first principles of classical hydrodynamics and fluid-particle dynamics and clarified how previous theoretical formulations relate to practical aspects seen in fluidization, sedimentation, and porous media. This text includes the formulation of Brenners solution to a rigid sphere approaching an infinite plane wall, as well as insight into formulations beginning with the Navier-Stokes equations, for other particle motions. Because of the slow moving nature of the problem, existing analytic solutions are limited to constant fluid viscosity and rigid particles. No solution currently exists for fluids with pressure-dependent viscosity or elastic particles. Depending on the motion of the binary particles and the type of fluid, high pressures may arise in the flow held and affect the viscosity of the fluid and the deformation of the particles.
8


CHAPTER II
ANALYTIC SOLUTIONS
Many well known analytic solutions have been developed to describe the forces and torques acting on particles in viscous fluids. A brief overview of these formulations will give insight into the formulation of binary particle interactions. All analytic solutions mentioned are limited to isothermal, low Reynolds number, incompressible flows, and assume an infinite fluid domain. These analytic solutions are important for the verification of numerical simulations, where numerical error can be introduced. We can identify if we are using the correct boundary conditions, an adequate mesh, and a correct domain size by using these analytic solutions. Once the numerical solutions agree with the analytic solutions, other parameters such as fluid viscosity and particle elasticity can be modified.
2.1 Related Benchmark Solutions
In 1851, George Stokes developed Stokes Law, which is an exact solution to a spherical particle moving through an infinite volume of viscous fluid at a constant velocity. The law is the basis for the falling ball viscometer [6] and can be expressed as
Fd =-6n:/irU, (2.1)
where Fd is the drag force on the sphere, n is the dynamic viscosity of the fluid, r is the radius of the sphere, and U is the velocity at which the particle is moving. In 1961, Howard Brenner considered the case where a spherical particle is close enough to a wall where the drag force is affected. Brenner introduced an expression describing the force on a spherical particle approaching an infinite plane wall as
F = Gir/irUX, (2.2)
where A is a modification factor to Stokes Law [7]. The wall is simulated as a no slip
boundary condition, or a boundary at which the fluid has zero velocity relative to the
9


boundary. Brenners solution is in the form of an infinite series, where A is expressed as
A = sinh a
n(n +1) 2 sinh (2n + l)ci! + (2n + 1) sinh 2a
( (2n 1)(2n + 3) _4sinh2 (n + \)a (2n + l)2 sinh2 a
- 1
(2.3)
and a as
with h being the distance from the center of sphere to the wall. The analytic solution for the torque on a sphere rotating in an infinite domain of viscous fluid is expressed as
T = 87r/rr3fl, (2.5)
where is the rotational velocity in ^ [6]. Equations 2.1 and 2.5 will frequently be used in later sections for normalizing forces and torques.
2.2 Solutions to Four Modes of Motion
Figure 1.1 shows the four possible modes of relative motion between two particles, with labels that will be used throughout the thesis. The equations for mode 1, describing the torques acting on binary spheres rotating about an axis perpendicular to their line of centers, are presented by Jeffrey & Onishi [19]. The solution is for the case where one sphere is rotating and the other is stationary. As a modification to equation 2.5, the torque on the rotating sphere with an angular velocity of is
Trotating 87T HTrVlgr (2-6)
and for the stationary sphere,
T-'stationary 87T jlTgVlg,
(2.7)
10


The modification factors, gr and gs are
and
9r
2 In e 5(1-A)
B n(A)
66 12A + 16A2 125(1 A)2
elne + 0(e)
(2.8)
9s
A2 In e 10(1 A)
-812(A)
A2 (43 + 24A + 43A2) 250(1 A)2
elne + 0(e),
(2.9)
respectively, where e is the minimum separation distance between spheres. The values for £>11 and BV2 were obtained in [19]. The A term in equations 2.8 and 2.9 relate the radius of the two spheres, where in this case these are equal,
A =
T rotating T stationary
-1.
For this case where the radius of the spheres are equal, Bn
(2.10)
.7029 and Ba = .0274.
Figure 2.1: Torque on binary spheres, one rotating about an axis perpendicular to their line of centers and one stationary, normalized by the torque on an isolated rotating sphere.
The hydrodynamic torque on the rotating and stationary spheres increases with decreased separation distance, as shown in figure 2.1. The singularity for this mode of
11


motion is log(e), where the singularity describes how the function approaches a singular solution. The torque values are normalized by the torque on a single isolated sphere, rotating in the same fluid at the same angular velocity as the rotating sphere in the binary case.
The equations to solve for the forces on binary spheres translating along an axis perpendicular to their line of centers, or mode 2, are outlined by Jeffrey & Onishi [19]. The solution is for the case where one sphere translates at a velocity of U while the other is stationary. The equations for each sphere are
F'translating -6nfirtUft
(2.11)
and
(2.12)
respectively. The modification factors for Stokes law are
4(16 + 45A + 58A2 + 45A3 + 16A4)e/ue 375(1 A)4
f0(e) (2.13)
and
(2.14)
The values for A21 and A22 were obtained by Jeffrey & Onishi, and for the case of equal
radii, are A21 = 0.9983 and A22 = 0.27 37.
12


3
-^-Translating Sphere
Figure 2.2: Force on binary spheres, one translating along axis perpendicular to line of centers and one stationary, normalized by Fstokes-
The hydrodynamic force on binary particles translating relative to one another along the axis perpendicular to their line of centers is shown in figure 2.2. The singularity for this motion is log(e), similar to mode 1. The forces are normalized with Stokes drag, from equation 2.1, with the same fluid viscosity and at the same particle velocity as the binary spheres.
In 1984, Jeffrey & Onishi published a solution for binary spheres with a relative rotating motion about their line of centers [19], or mode 3. The formulations for both the rotating and stationary spheres that are introduced as a modification to equation 2.5 can be expressed as
and
Trotating
87rg,rfahr
(2.15)
Fstationary 87T gTsVlhs, (2.16)
respectively. For the rotating sphere, the modification factor is
13


h
(2.17)
((3,(1-A) 1) dnt
(1 A)3 2(1 A)2 ^ U 6
and for the stationary sphere, the modification factor is
A3((3,1) A3dne
(1 A)3 2(1-A)2 ^
The function for ( is written as
(2.18)
CO.a) = y^(A- + a) (2.19)
k=0
Figure 2.3: Torque on binary spheres, one rotating about line of centers and one stationary, normalized by the torque on an isolated rotating sphere.
The normalized hydrodynamic torque acting on binary spheres where one rotates about their line of centers and the other is stationary is shown in figure 2.3. There is no singularity for this mode of motion and the torque reaches a constant value as the gap between the spheres vanishes.
In 1961, Brenner developed a solution for a sphere approaching a free surface [7]. A
free surface is used to express a symmetry boundary condition which reduces numerical
14


costs by reducing the degrees of freedom. A free surface is expressed as the vanishing of the normal component of velocity and tangential shear stresses. This is the solution to binary spheres traveling along their line of centers, which we have named mode 4. The formulation is expressed as
F = 6TT^rUj3,
where /5 is a modification factor to Stokes Law, expressed as
(2.20)
0
4 \ n(n + 1)
- smh a > ------------
3 ^ 2 /?,-12 /?, +3
n= 1
4 cosh2 (??, + F)a + (2n + l)2 sinh2 a 2 sinh (2n + 1)q' (2n + 1) sinh 2q-
1 .
(2.21)
Figure 2.4: Force on binary spheres translating along their line of centers, normalized by
FStokes
The normalized force of binary spheres traveling along their line of centers is shown in figure 2.4. The solution to this motion is divergent, and the singularity can be described as 4.
e
15


CHAPTER III
NUMERICAL SIMULATIONS
Numerical solvers have grown to be an extremely powerful method to solve many types of mathematical equations. Physics problems are often governed by nonlinear equations, which are often impossible to solve analytically. Numerical techniques can often be used as a substitute. COMSOL Multiphysics is a commercial software package which uses numerical methods to model and simulate many types of physics problems. The finite element method was used for the simulations in this report, where the software solves the equations with different numerical methods, depending on whether there is a combination of multiple physics or a single physics phenomenon. The fully-coupled approach combines the different physics into one matrix that is solved, compared to the segregated solver that solves one matrix for one phenomenon and applies the solutions to a separate matrix for the other phenomenon. The software also allows for the selection of an iterative solver and a direct solver. Direct solvers rely on a large amount of computer memory and the problem must be small enough to invert the matrix all at once, while iterative solvers require multiple iterations to break the problem up into smaller computations. For the simulations in this thesis, the nonlinear nature of the problems and the large amount of degrees of freedom required the use of iterative solvers.
COMSOL also has a parametric sweep feature that allows for the consecutive solving of problems with different specified parameters. This function was used to consecutively solve for a range of particle separation distances and particle elasticities. When the separation distance is very small, the pressures can be extremely high in the interstitial region and it can be difficult to solve the governing equations. To remedy this, COMSOL presents an option to use a previous solution as initial conditions for the current simulation. By default, COMSOL assumes initial conditions for fluid velocity and pressure to be zero, which is far from the solution. By using a previous solution as the
16


initial condition in the solver, we can facilitate convergence by starting with an initial solution closer to the final solution.
COMSOL allows the user to specify material properties for both fluids and solids. This can be done by using constant values or by specifying a function to govern the material property based on one or several input parameters. A pressure-dependent analytic function was used to determine the fluid viscosity in the simulations investigating pressure-dependent viscosity. In the particle elasticity simulations, the particle elasticity and Poissons ratio were varied to determine their affect on the binary particle hydrodynamics. To first obtain an understanding of the software environment, a spherical particle in Stokes flow was simulated. A sphere with a radius of 0.001 m moving through a fluid with a viscosity of 1 Pa s was used. The numerical drag force on the sphere was compared to the Stokes equation and a mesh convergence study was performed. The problem was set up with a 2-dimensional axisymmetric geometry, where the sphere is translating along the axis of revolved symmetry. This significantly reduces the number of degrees of freedom in the simulation. The cylindrical fluid region containing the sphere was defined as lOr wide and 20r tall. The fluid domain was bounded by an open boundary on the circumferential wall of the cylinder, an inflow boundary on the top of the cylinder, and an outflow boundary on the bottom of the cylinder.
17


Table 3.1: Mesh convergence study results for the drag force on an isolated sphere in
Stokes flow.
No. of Elements FDrag [A] % Change
110 .01295 N/A
169 .01325 2.23
538 .01638 19.15
2094 .01744 6.05
4658 .01853 5.90
8254 .01882 1.56
12902 .01887 0.24
Table 3.1 shows the results of a mesh convergence study performed on the simulation of Stokes law. A coarser mesh is used to begin with because the computation time can be significantly reduced if the number of degrees of freedom is minimized. Only 110 finite elements were initially used in the mesh. The mesh density was gradually increased until the difference in solution between two successive meshes was below 1%. At this point, the difference in results is so insignificant between each new simulation, that we can satisfy the mesh convergence criteria. The highest density mesh in the table results in a drag force of 0.01887 N, differing by 0.1 % from the analytic result of 0.01885 N obtained from equation 2.1. The exact solution to the problem is known in this case, but when new problems are being solved, a mesh convergence study is often necessary to guarantee convergence of the solution.
A rigid sphere in a viscous fluid approaching an infinite plane wall was also solved numerically. The results for the force on the sphere were compared to the analytic solution in equation 2.2. This problem helped gain insight into setting up the numerical parameters for the simulations to follow.
18


Figure 3.1: Comparison of the analytic and numerical solutions for a rigid sphere approaching an infinite plane wall.
The analytic and numerical solutions for the force on a sphere approaching an infinite plane wall is shown in figure 3.1. The error is insignificant such that the curves appear coincident, so this data was also tabulated.
19


Table 3.2: Relative error of numerical simulations for a sphere approaching a plane wall, compared to the analytic solution.
Separation, e % Error
5r 2.95
r 1.77
r/2 1.18
r/10 0.34
r/100 0.036
r/1000 0.005
r/2000 0.002
r/5000 0.002
r/10000 0.005
Table 3.2 shows that the numerical results agree with the analytic solution for the problem within 3% at the largest separation distance of 5r and within .005% at the smallest separation distance of raf0Ts Brenners analytic solution assumes the sphere is approaching an infinite plane within an infinite fluid domain. In numerical simulations using the finite element method, the problem must be set up using a finite domain enclosing all of the boundaries, bodies, and physics. The region size must be defined and cannot be infinite. It was found that the larger the region is, the more accurate the solution becomes, but at the cost of a significant increase in degrees of freedom and computation time. This quickly becomes a diminishing return in terms of accuracy. Therefore, a decision must be made as to how large of a region is adequate for the problem.
COMSOL allows for the use of different modules to determine the equations to be used to solve the underlying physics. For the pressure-dependent viscosity investigation with rigid particles, the creeping flow module was used, which uses Stokes equations
20


for incompressible fluids. For the problems involving elastic particles, the fluid structure interaction, or FSI module was used. This module allows for the combination of two different modules; creeping flow for the fluid and structural mechanics for the solids. The FSI module couples the two problems so the solid deformations affect the fluid pressures and velocities and vice-versa. The software combines the multiple physics by using a segregated solver, segregating each of the governing equations and solving for them individually.
3.1 Mesh Discretization
The finite element method relies on having an appropriate mesh defined, where values are calculated at individual points within a mesh. The mesh consists of elements, which are made up of nodes. For the common tetrahedral element, there are four vertices, or nodes, where the solution could be computed. These values might be deformations, temperatures, velocities, or any number of quantities. Other element types are the bar element, typically used in two-dimensional structural problems, and the three-dimensional quadrilateral element, which has 8 nodes.
There would be gaps in the solution if the values were only calculated at the nodes. Interpolation functions are used to approximate the solution between the nodes using a polynomial function. In all simulations in this thesis, second order functions were used for velocity, while linear shape functions were used for pressure. Once an element type and discretization orders have been chosen, the overall mesh geometry must be created.
21


Figure 3.2: Example of mesh for binary spheres traveling along their line of centers.
The fluid regions that have the largest gradient in velocity and pressure must have a higher mesh density, whereas the areas with small gradients across elements may have larger elements. It is possible to know before creating a mesh that a certain region will have large gradients in pressure and velocity. In figure 3.2, there are two spheres that are close together and moving along their line of centers, so we expect to see high pressures in the fluid and the mesh should be fine in this area. In general, it is best to have more than two elements of thickness in thin areas such as this interstitial region.
Figure 3.3: Detailed view of interstitial region mesh for binary spheres.
High pressures within the interstitial region fluid are likely for binary spheres moving
along their line of centers. In figure 3.3, we see that the region of most importance is
where the gradients in pressure and velocity are the largest, between the two spheres.
22


Higher density meshes with a greater number of degrees of freedom come at the expense of computation time.
3.2 Boundary Conditions
Choosing the correct boundary conditions is very important when setting up a numerical simulation. Boundary conditions are the mathematical representation of physics at interface boundaries in a simulation. If fluid is traveling adjacent to a solid body, there is a no-slip condition at the surface of the body. A no-slip condition specifies that the fluid velocity is zero relative to the solid boundary. The no-slip boundary condition is given by
u = uw,v = vw,w = ww, (3-1)
where uw,vw, and ww are the x, y, and z components of the walls velocity, respectively. This boundary condition is used at any fluid-solid interface, so the fluid has the same components of velocity as the solid, where they are in direct contact. The no-slip condition was used in all of the numerical simulations in this thesis. Another example of a boundary condition that is often used in numerical simulations is the symmetry condition. A symmetry boundary condition specifies that the wall has zero shear stresses, which is equivalent to having a symmetric boundary and effectively mirroring any physics on one side of the boundary to the other. This condition is given by
r = 0 at the symmetric boundary. (3.2)
Symmetry boundaries become important when one wants to reduce the computation time it takes to solve a problem. In the case where two spheres are approaching one another, such as the rightmost case in figure 1.1, a symmetric plane directly between the two spheres would cut the degrees of freedom to be solved in half. On larger problems, the difference becomes significant.
23


Each mode of motion that was simulated required slightly different boundary conditions. For the rotation about an axis perpendicular to the spheres line of centers, the surface of the spheres both required separate boundary conditions. The sphere that is stationary simply requires a no-slip boundary condition, as described in equation 3.1.
The outflow boundary conditions used in the simulations vary depending on the particular motion of spheres. In the simulation of binary spheres rotating about an axis perpendicular to their line of centers, a no-slip boundary was applied on the circumferential surfaces of the cylindrical container. On the flat ends of the cylinder, an open boundary condition was applied, which specifies the boundary is in contact with a large volume of fluid, and viscous stresses vanish. The same boundary conditions were applied for binary spheres translating relatively about an axis perpendicular to their line of centers. These outflow conditions were chosen because with all open boundaries, the steady state solution was that the far-held fluid velocity was nonzero, which is inaccurate. This significantly reduced the resistance on the spheres, to a value below the analytic solution. Imposing a no-slip condition on the circumferential container walls yielded very accurate results in regards to the analytic solutions.
For binary spheres rotating about their line of centers, a no-slip condition was applied on all container walls. This was done so the far-held velocity was zero, rather than a velocity generated by the rotation of the sphere, similar to mode 1. For binary spheres approaching one another along their line of centers, when the spheres are close, the majority of the force is from the high pressures in the interstitial region. The effect of the outflow boundary condition was found to be negligible and an open boundary was arbitrarily chosen.
24


CHAPTER IV
FLUID CONSTITUTIVE EQUATIONS
The two main areas of focus for this research are pressure-dependent viscosity and particle elasticity and how each affects the hydrodynamic forces and torques acting on binary particle interactions. The constitutive equations for a fluid are those that relate the physical quantities of such as viscosity, density, and pressure, to one another. Two models were studied and compared for determining the viscosity of a fluid based on its pressure. One model was chosen to use as our governing equation for a fluids viscosity, based on literature [17].
4.1 Viscosity Models
Models relating pressure to viscosity date back to the late nineteenth century. Barus developed the first widely used model for the pressure dependence of viscosity in fluids [17]. Barus stated that the relationship between pressure and viscosity can be expressed as
Inyj = ap, (4.1)
where p is the viscosity, pQ is the absolute viscosity at a gauge pressure of 0 Pa, p is the current pressure, and a is the pressure-viscosity coefficient. For over 70 years, this model was accepted as the standard, until it was proved to significantly over predict the viscosity at pressures above approximately 10s Pa. In 1966, another model was developed by Roelands to correct this [17]. Roelands modeled the viscosity at enhanced pressures as
logp + 1.20 = (logpo + 1.20) (l + , (4.2)
where z is a dimensionless viscosity-pressure index related to a. The Barus equation
is said to be accurate up to approximately 0.1 GPa, while the Roelands equation has a
25


much more accurate prediction for viscosity at pressures up to approximately 1.5 GPa [4]. To study the differences between the models, the viscosity of mineral oil was examined for each.
10
15
JO
CO
o
o
GO
"> 5
.2 10 E
n3
c
>>
Q
10
O Barns x Roelands

O
O
o
o
o
o
o
o
o
X
x
O x
O x X
O x
x
0.5 1
Gauge Pressure [GPa]
1.5
Figure 4.1: Comparison between Barus and Roelands models for viscosity of mineral oil.
In figure 4.1, it is evident that there is no limit to the viscosity of the mineral oil according to the Barus equation. We see that above approximately 0.1 GPa, the Roelands curve begins to decrease in slope and eventually level off at very high pressures. According to experimental data, this is a more accurate representation of how the viscosity of fluids behaves at high pressures [17], and we will proceed using this model. Four different fluids were chosen to plot the Roelands model, for a diverse comparison of different parameters.
Table 4.1: Viscous properties of four different liquids.
Fluid Ho [Pa s] z
Water 0.001 0.1
Polyethylene glycol 0.09 0.38
Glycerin 1.412 0.18
Mineral Oil 0.0681 0.67
26


Table 4.1 shows how the dimensionless pressure-viscosity index, z, is unique to each fluid. This parameter determines how the fluid viscosity increases with pressure.
It is known that water is considered an incompressible fluid, and in figure 4.2, we can see that the viscosity is approximately constant throughout the pressure range. The Roelands model was developed for mineral oils [17], and we can see that mineral oil has the most notable increase in viscosity at enhanced pressures out of all fluids in the plot. Polyethylene glycol and glycerin have a relatively low increase in viscosity with pressure, in comparison. Glycerin is a highly viscous fluid, as we see when comparing its viscosity at ambient pressure to any other fluid in table 4.1. It becomes evident just how significantly pressure can increase viscosity when considering the difference in viscosity between glycerin and mineral oil, in figure 4.2. Mineral oil surpasses the viscosity of glycerin at a gauge pressure of approximately 0.2 GPa.
27


CHAPTER V
RESULTS: PRESSURE DEPENDENT VISCOSITY
The viscosity of a fluid is generally dependent on both pressure and temperature. The analytic solutions for the modes 1-4 assume the fluid has a constant viscosity. However, the viscosity can vary significantly depending on the pressure held in the domain. In this chapter, we consider the pressure dependence of viscosity for isothermal systems. The fluid viscosity at each degree of freedom was calculated using equation 4.2. The hydrodynamic forces for all four modes of motion in figure 1.1 were calculated using numerical techniques, for constant and pressure-dependent viscosity fluid. Mineral oil was used as the fluid because its viscosity has the greatest dependence on pressure, compared to the other fluids in figure 4.2. The elements in the mesh were quadratic in velocity and linear in pressure.
The governing equations used in the numerical simulations assume low Reynolds number how conditions (viscous effects dominate the how and the inertial effects can be ignored). The binary particle simulations were performed at constant particle velocities, as opposed to constant loads.
5.1 Rotation About Axis Perpendicular to Line of Centers
The effect of pressure dependent viscosity where spheres rotate relative to one another about an axis perpendicular to their line of centers (mode 1) is analyzed. Numerical simulations were performed to provide a direct comparison between pressure dependent and constant viscosity huids. The simulations rehect a system of binary spheres with one rotating and one stationary, where the spheres are at four different separation distances. The simulations for modes 1 and 2 were three-dimensional, with the spheres occupying a cylindrical huid domain of radius 15rsphere and height 30rsphere- Very dense meshes were constructed in the interstitial region, with elements growing in size toward the bulk region of the huid.
28


Figure 5.1: Geometry used in mode 1 and 2 simulations.
The geometry used for the mode 1 and 2 simulations is shown in figure 5.1. The fluid region is a cylinder bisected by a symmetric plane.
Figure 5.2: Detailed view of the geometry surrounding the spheres for modes 1 and 2.
The interstitial region geometry was divided into four concentric cylindrical regions
centered along the spheres line of centers, each with different mesh densities. This
geometry is show in figure 5.2. The smallest of the cylinders had a radius of
followed by a cylinder of radius rspg0ere, rsph5ere 1 and rsph2ere. The smallest cylinder had a
very fine mesh where the highest fluid pressures and velocities were present. The mesh
29


density grew as the cylinders grew because the pressure and velocity gradients decreased far from the interstitial region.
Figure 5.3: Overview of the mesh used for modes 1 and 2.
The mesh used in the mode 1 and 2 simulations is shown in figure 5.3. In the smallest of the interstitial cylinders, triangular elements were defined with a maximum size of | and swept in the radial direction with a distribution of 40 elements per 90 degrees, producing triangular prisms. The next smallest cylinder was given a maximum element growth rate of 1.01. The next smallest was given a maximum element growth rate of 1.04. Finally, on the largest of the interstitial cylinders, a transition was made from the triangular prism elements to tetrahedral elements with a COMSOL predefined size of finer. The remaining bulk region was given a predefined element size of coarse, since a dense mesh is not needed outside of the interstitial region.
30


Figure 5.4: Detailed view of the mesh used for modes 1 and 2. A high mesh density region is defined close to the interstitial region and the mesh grows coarser toward the bulk fluid region.
A detailed view of the interstitial region of the mesh is shown in figure 5.4. The finest elements are in the smallest interstitial cylinder and grow outward toward the bulk of the fluid. The highest mesh density is in the interstitial region where the spheres are the closest since the highest gradients in pressure and velocity are present here. The swept elements appear as quadrilateral from the interior of the sphere because the sides of the triangular mesh created on the symmetry plane are swept in the radial direction.
Table 5.1: Details of the numerical mesh scheme used for mode 1 simulations at each separation distance.
e No. Elements Deg. of Freedom Computation Time [hr]
rad/10 224683 1383563 1.00
rad/102 158772 1004365 0.75
rad/103 194733 1575838 1.25
rad/104 354150 3875544 2.25
The resulting mesh statistics for the mode 1 simulations are shown in table 5.1. The
larger the degrees of freedom in a mesh, the longer the computation time becomes. The
31


smallest separation distance results in the greatest degrees of freedom and computation
time, requiring approximately 2.25 hours.
Table 5.2: Analytic and numerical solutions for mode 1 rotating sphere in fj, = const. mineral oil, normalized by the torque on an isolated sphere.
Separation, e TNumeric / Tj solated TAnalytic/ Tj solated Percent error
r/10 1.229 1.233 0.35
r/102 1.632 1.635 0.16
r/103 2.085 2.086 0.06
r/104 2.546 2.545 0.02
The torque obtained from the numerical and analytic solutions are compared in table
5.2 for constant viscosity mineral oil. The table shows the torque on the rotating sphere for each solution and the error in the numerical solutions. The error is very low at every separation distance and this indicates the mesh is adequate to introduce pressure-dependent viscosity to the fluid.
Table 5.3: Numerical solutions for fj, = const, and fi = fi{p) numerical solution for mode 1 rotating sphere in mineral oil, normalized by the torque on an isolated sphere.
Separation, e Tji=consi /TIsoiated Ti_l=i_l(^p'j /TIsoiated Percent difference
r/10 1.229 1.229 0.0
r/102 1.633 1.633 0.0
r/103 2.085 2.085 0.0
r/104 2.546 2.546 0.0
The torques on the mode 1 rotating sphere in constant and pressure dependent viscosity mineral oil are compared in table 5.3. The data in the table suggests that mode
1 produces insufficient fluid pressures to increase the fluid viscosity significantly, even
32


at very small separation distances. The percent difference column shows that introducing pressure-dependent viscosity did not affect the torque on the spheres, compared to constant viscosity fluid.
Figure 5.5: Comparison of fi = p(p) and fi = const, numerical results and fi = const, analytic solution for mode 1. Fluid is mineral oil, particle velocity is Q = 1^.
The data in tables 5.1 and 5.2 are combined in figure 5.5 to compare the analytic solution to the constant and pressure dependent numerical solutions. All three curves he on top of each other, showing that the numerical results correlate very well with the analytic solution with minimal error and that the constant and pressure-dependent viscosity fluids yield the same torques.
Table 5.4: Comparison of maximum gauge pressure and viscosity in constant and pressure dependent viscosity mineral oil, for mode 1 particle motion with VLrot = anc[ e =
radius 104 *
Viscosity Model Pgauge,max. [-P&] t^max. \P& s]
fi = const. 26550 0.0681
Roelands, fj, = p(p) 26558 0.0681
Table 5.4 contains the maximum fluid gauge pressure and viscosity in the fluid re-
33


gion for the constant and pressure-dependent viscosity simulations. The data is for the numerical simulations of binary spheres in mineral oil with a relative angular velocity of l2/^ and a separation distance of ra^s. There are no significant changes in fluid pressure or viscosity, which reinforces the data in table 5.3.
5.2 Translation Along Axis Perpendicular to Line of Centers
Binary spheres translating relatively along an axis perpendicular to their line of centers (mode 2) is considered. One sphere was given a constant velocity motion perpendicular to the line of centers, while the other sphere was stationary. We consider the force on the moving sphere to analyze the hydrodynamic effects. A three-dimensional simulation was performed, using the same geometry and mesh as the simulations in section 5.1, with very similar computation times.
Table 5.5: Comparison of analytic and numerical solutions for mode 2 moving sphere in H = const, mineral oil, normalized by the force on a single sphere in Stokes flow.
Separation, e FNumeric/ Fl solated FAnalytic/ Fjsoiatey Percent error
r/10 1.389 1.393 0.29
r/102 1.761 1.767 0.36
r/103 2.162 2.150 0.55
r/104 2.515 2.533 0.74
We see a very low relative error at every separation distance by comparing the analytic and numeric solutions to the force on the moving sphere in table 5.5. This demonstrates that the mesh being used is adequate to capture the physics of mode 2.
34


Table 5.6: Comparison of fi = const, and fi = p.(p) numerical force on mode 2 moving sphere, normalized by the force on a single sphere in Stokes flow.
Separation, e Ffi=constJ Fisolated F/j,=/j,(p) /Fjsolated Percent difference
r/10 1.389 1.389 0.0
r/102 1.761 1.761 0.0
r/103 2.162 2.162 0.0
r/104 2.515 2.515 0.0
Table 5.6 compares the numerical results for constant and pressure dependent fluid viscosity. The results show that the two cases have nearly identical results and the pressure in the fluid is not sufficient to increase the viscosity for the pressure-dependent fluid simulation.
analytic solution for mode 2. Fluid is mineral oil, particle velocity is U = 1.
Figure 5.6 shows the correlation between the mode 2 constant and pressure-dependent viscosity numerical solutions and the analytic solution. The curves are all coincident, meaning that the error in the constant viscosity numerical solution and the difference between the constant and pressure-dependent viscosity simulations are negligible.
35


Table 5.7: Comparison of maximum gauge pressure and viscosity in constant and
pressure-dependent viscosity mineral oil, for mode 2 motion with U = 1 and e =
Viscosity Model Pgauge,max. \P&] t^max. [P& s]
fi = const. 1.089 0.0681
Roelands, fj, = p(p) 3.4136 0.0681
Table 5.7 shows the maximum gauge pressure in the fluid for both the constant and pressure-dependent viscosity fluids, obtained from the numerical simulations. The results are for binary spheres in mineral oil with a relative velocity of 1 at a particle separation distance of ra^s. There is not a significant difference in fluid pressure and viscosity between the mode 2 constant and pressure-dependent viscosity fluid. The maximum pressures in table 5.7 are noticeably lower than those in table 5.4, when the spheres in both modes have the same tangential velocity at their closest point in the interstitial region.
Figure 5.7: Velocity vectors on the cross section of binary spheres for modes 1 and 2, from left to right.
The components of velocity for modes 1 and 2 are shown in figure 5.7. The upper
sphere in mode 1 is rotating about its center point, without any translation. The steady-
state solution to mode 2 was solved at the instant in time the upper sphere was traveling
36


perpendicular to the spheres line of centers.
Figure 5.8: Mode 1 (left) and 2 (right) diagram of pressure distribution. LP denotes low pressure, while HP denotes high pressure.
The distribution of pressures on the surfaces of the spheres near the interstitial region is shown in figure 5.8. The rotation in mode 1 entrains the fluid into the wedge-shape formed by the spheres surfaces, causing an increase in pressure. On the right side of the sphere, the surfaces of the spheres are separation, causing a low pressure. At the instant the moving mode 2 sphere is translating perpendicular to the spheres line of centers, it is dragging the fluid on its left side, creating a pressure drop behind it. This fluid then collides with the lower sphere, creating an increase in pressures on its surfaces. The opposite occurs on the right side of the line of centers, where the upper sphere is pushing the fluid body, causing an increase in fluid pressures. The fluid adjacent to the lower sphere is being pulled away from its surface, creating low pressures.
37


x-coordinate [m]
Figure 5.9: Pressure distribution on the surfaces of mode f spheres near the interstitial region.
The pressure distributions along the surfaces of the mode 1 spheres near the interstitial region is shown in figure 5.9. The curves are almost coincident, meaning the pressure in the interstitial region is a function of the x-coordinate only. The fluid is entrained into the wedge-shape made between the spheres, by the rotation of the upper sphere.
38


Figure 5.10: Pressure distribution on the surfaces of mode 2 spheres near the interstitial region.
The pressure distributions along the curved interstitial side of the mode 2 spheres are shown in figure 5.10. The figure shows that the fluid adjacent to the moving sphere is increasingly negative left of the line of centers, as we approach the interstitial region. At the point where the spheres are the closest, the pressure is 0, followed by an upward spike in pressure the right of the line of centers. The stationary sphere has exactly the opposite trend.
5.3 Rotation About Line of Centers
Pressure-dependent viscosity was simulated for binary spheres with a relative rotation about their line of centers. This mode of motion has no singularity and the analytic solution in figure 2.3 shows that the torque on the spheres stops increasing at small separation distances, below approximately Since the spheres are rotating relative to one another about their line of centers, the point where their surfaces are the closest has a velocity of 0 and no shearing will occur here. The simulation was two-dimensional
39


axisymmetric, with a cylindrical computational domain of radius 15rsphere and height
30 r
sphere *
Figure 5.11: Mesh scheme for modes 3 and 4, with three separate regions of different element densities.
Figure 5.11 shows the mesh scheme used for modes 3 and 4. To capture the potential gradients in fluid pressure and velocity in the interstitial region, several regions of different element densities were used. The mesh density is very high in the interstitial region and the elements grow toward the bulk region of the fluid. The smallest region in figure 5.11 has a maximum element size of |, the next smallest has a maximum element growth rate of 1.025, and the largest bulk fluid region has the predefined COMSOL mesh density of finer. With this mesh scheme, the convergence was met and the accuracy of the numerical constant viscosity solution at all separation distances was within 0.05% of the analytic constant viscosity solution.
40


Table 5.8: Details of the numerical mesh scheme used for mode 3 simulations at each
separation distance.
e No. Elements Deg. of Freedom Computation Time [s]
rad/10 7190 33158 4
rad/102 9271 43005 4
rad/103 11192 52667 4
rad/104 44862 210632 10
Table 5.8 shows details of the mesh used in figure 5.11 including the number of elements and degrees of freedom for each separation distance. The mesh becomes denser in the area between the spheres to account for the large pressure and velocity gradients in the interstitial region. The elements must grow at a low enough rate to avoid becoming stretched to an irregularly shaped triangle and reducing the mesh quality. This results in a much greater number of elements and degrees of freedom for the smaller separation distances meshes. The resulting mesh for e = has 44862 elements and 210632 degrees of freedom.
To determine the affect of pressure-dependent viscosity fluid on the hydrodynamic forces on particles, mineral oil with a constant and pressure-dependent viscosity was simulated numerically.
Table 5.9: Comparison of analytical and numerical solutions for mode 3 rotating sphere in mineral oil with n = const., normalized by the torque on an isolated sphere.
Separation, e TNumeric/Isolated TAnalytic / Tj solated Percent error
r/102 1.047 1.058 1.01
r/103 1.051 1.053 0.13
r/104 1.052 1.052 0.0
Table 5.9 shows the mesh is adequate to capture the physics for mode 3, by agreement
41


between the analytic solution and the numerical solution. Two spheres with a relative rotation about their line of centers have little to no shear stresses in the interstitial region between the spheres. The torque on the spheres is almost constant at different separation distances since the majority of the shear stress arises far from the interstitial region.
Table 5.10: Comparison of fi = const, and fi = p(p) numerical torque on mode 3 rotating sphere in mineral oil, normalized by the torque on an isolated rotating sphere.
Separation, e Tfx=const. fTl solated 7/i=/i(p) /Tlsolated Percent difference
r/102 1.047 1.047 0.0
r/103 1.051 1.051 0.0
r/104 1.052 1.052 0.0
Table 5.10 shows that the introducing pressure-dependent viscosity to the mineral oil does not result in a torque increase. This is because the velocity on the surface of the spheres is negligible in the interstitial region and the pressure and viscosity of the fluid remain constant.
1.2
1.15
1.1
1.
05^1
0.95
0.9L
10'4
-H-
10
,-3
Analytic X Numeric, ft = const. Numeric, m = m(p)
E3
10
,-2
Separation Distance/ Radius
Figure 5.12: Comparison of fi = p(p) and fi = const, numerical results and fi = const, analytic solution for mode 3. The fluid is mineral oil and the rotating particle angular velocity is Q = 1 .
42


Figure 5.12 shows that the analytic and numerical solutions to mode 3 agree very well at close separation distances and that introducing pressure-dependent viscosity does not increase the torque on the particle. This is because the fluid pressures increase negligibly and there is no change in viscosity. The slight deviation in the analytic and numerical solutions at the largest separation distances is because the analytic solution is limited to very close separation distances.
Table 5.11: Comparison of maximum gauge pressure in constant and pressure-dependent viscosity mineral oil, for mode 3 particle motion with 12 = 1^ and e =
Viscosity Model Pgauge^max. \P&] l^max. [PCL s]
H = const. 0 0.0681
Roelands, fj, = fi(p) 0 0.0681
Table 5.11 contains the maximum gauge pressure and viscosity in the fluid for the numerical simulation of binary spheres in constant and pressure-dependent viscosity mineral oil with a relative angular velocity of l2^ and separation distance of ra^'s. The numerical simulations determined the pressure in the fluid was identical in the constant and pressure-dependent viscosity fluids.
5.4 Translation Along Line of Centers
Binary spheres translating along their line of centers (mode 4) results in the highest fluid pressures out of all four modes of motion. The singularity for this motion is ^ and the increase in force is caused by the significant increase in fluid pressures between the spheres. The high pressures in the interstitial region act on the surfaces of the spheres to increase the hydrodynamic forces on them. The numerical simulation for this motion was a two-dimensional axisymmetric formulation, using a geometry and mesh similar to the simulations as discussed in section 5.3. The only difference is that mode 4 can be solved axisymmetrically and we can use symmetry between the spheres, on a plane
43


perpendicular to their line of centers.
Table 5.12: Details of the numerical mesh used for mode 4 simulations at each separation distance.
e No. Elements Deg. of Freedom Computation Time [s]
rad/10 6797 31412 4
rad/102 9949 46201 4
rad/103 13214 62226 4
rad/104 65284 304581 12
The mesh statistics for each separation distance simulated are shown in table 5.12. The closer the spheres are to one another, the denser the mesh must be in the interstitial region to capture the large gradients.
7000
6000
5000
4000
3000
2000
1000
0
I


\
\
\
\


- Analytic x Numeric,const. Numeric, |ii = jw(p)
£3 ^
i
10'4 10'3 10'2 10"1
Separation Distance / Radius
Figure 5.13: Comparison of fi = p(p) and fi = const, numerical results and fi = const, analytic solution for mode 4. Fluid is mineral oil, particle velocity is U = 2.15 .
The hydrodynamic forces increase substantially as the spheres become closer, as shown
44


in figure 5.13. The analytic solution for this motion is plotted to show that the mesh was adequately converged for the simulations, by its agreement with the numerical constant viscosity results. At this particle velocity of U = 2.15 in mineral oil with pressure-dependent viscosity, the hydrodynamic force is 39% higher at the smallest separation distance of ra^'s, compared to constant viscosity fluid.
Table 5.13: Mode 4 force, max. pressure, and max. viscosity for n = const, and n = p(p) mineral oil at multiple approach velocities and a separation of e = .
U [m/s] H = const. Roelands, n = P(P)
Pmax. [P&] F/ Fstokes Pmax. [P&] t^max. F/ Fstokes
0.25 5.108e6 5005.39 5.446e6 0.077 5112.28
0.5 1.021e7 5005.39 1.169e7 0.089 5230.53
0.75 1.532e7 5005.39 1.902e7 0.106 5363.54
1 2.043e7 5005.39 2.785e7 0.130 5515.42
1.5 3.065e7 5005.34 5.367e7 0.233 5906.57
2 4.086e7 5005.24 1.183e8 0.914 6567.58
2.1 4.291e7 5005.42 1.586e8 2.050 6800.15
2.15 4.393e7 5005.33 2.058e8 5.026 6961.23
Table 5.13 shows the difference in force on the particles in a constant viscosity and pressure-dependent viscosity mineral oil using the Roelands model, at the smallest separation distance of ra^s. The table covers a range of approach velocities up to 2.15 m/s, beyond which convergence was unobtainable. We see that the greater the particle velocity, the greater the change in fluid pressure and viscosity. This leads to larger changes in force for the pressure-dependent viscosity fluid, compared to constant viscosity.
45


Table 5.14: Increases in mode 4 pressure, viscosity, and force from n = const, to n = nip) mineral oil, at e = ra^'s and a range of approach velocities.
U [m/s] Percent increase in maximum pressure Percent increase in maximum viscosity Percent increase in force
0.25 6.6 13.8 2.1
0.5 14.5 31.8 4.5
0.75 24.2 56.3 7.2
1 36.3 91.5 10.2
1.5 75.1 241.7 18.0
2 189.4 1242.4 31.2
2.1 269.7 2910.6 35.9
2.15 367.4 7281.4 39.1
Table 5.14 contains the relative differences in pressure, viscosity, and force for mode 4 motion between pressure-dependent and constant viscosity mineral oil. The maximum viscosity for the pressure-dependent simulation is approximately 7300% greater than that of the constant viscosity simulation at the highest approach velocity simulated. The resulting higher pressures lead to an approximate 39% greater force than the constant viscosity simulations at this approach velocity of 2.15 m/s.
46


7000
6500
O
6000
tuo
ro
5500
5000
0
0.5
1.5
2.5
Approach Velocity [m/s]
Figure 5.14: Mode 4 approach velocity vs. normalized drag force for pressure-dependent viscosity mineral oil at a separation distance of e =
Figure 5.14 shows the increase in the force on mode 4 particles, normalized by the drag force on a single particle in Stokes flow. The fluid pressure and viscosity are dependent on the approach velocity, which also results in changes in hydrodynamic force on the particles.
47


i xlO8
2
as
Q.
as
L_
in
as
0
Figure 5.15: Numerical results for mode 4 approach velocity vs. viscosity and pressure for pressure-dependent viscosity mineral oil at a separation distance of e = ^.
The increase in fluid pressure and viscosity based on the mode 4 approach velocity for pressure-dependent mineral oil is shown in figure 5.15. The particles are at a separation distance of e = ra^4S and we can see that as the velocity increases, the pressure and viscosity increase nonlinearly.
5.5 Combined Motions
Mode 4 motion generates sufficient fluid pressures in the interstitial region of the spheres to increase the fluids viscosity significantly, while modes 1-3 generate insufficient pressures to increase the viscosity. The motion of particles in suspensions is realistically a combination of the four modes, so we will investigate whether the combination of mode 4 with modes 1-3 affects the hydrodynamic forces on these modes. Combining mode 4 with the other motions imposes the high pressures and the change in fluid viscosity of mode 4 on the other modes, to determine their affect on each mode.
48


Mode 1,4 Mode 2,4 Mode 3,4
! i
i |
i !
Figure 5.16: Individual motions of modes 1-3 combined with mode 4 motion.
We will consider each combined motion in figure 5.16. Each combined motion requires a different geometry since mode 3 and 4 are solvable in a two-dimensional axisynnnetric configuration, while the others require three dimensions. We can make use of symmetry when combining modes 1 and 4 on a plane bisecting the spheres, perpendicular to the axis of rotation. Combining modes 2 and 4 also requires a three-dimensional geometry, but we can again make use of symmetry on a plane bisecting the spheres. Finally, modes 3 and 4 were solvable with a two-dimensional axisynnnetric geometry, since COMSOL allows an angular velocity about the axis of symmetry. Since the highest increases in viscosity were observed at e = the comparisons were performed at this separation distance, also in mineral oil.
The same three-dimensional mesh described in section 5.1 was used for combining modes 1 and 4 as well as modes 2 and 4. For combining modes 3 and 4, the two-dimensional axisynnnetric geometry and mesh from section 5.3 was used. The computation times for all combined motions were very similar to those in these previous sections, with their respective meshes. Modes 1 and 4 were combined by specifying a moving boundary condition on each sphere toward one another along the line of centers, while also specifying a rotational boundary condition on one of the spheres about the axis perpendicular to the line of centers. For mode 2, the same approach velocities were given
to mode 4 and a moving boundary condition was given to mode 2. Mode 3 is similar
49


to mode 1 in that the spheres are approaching one another, and one sphere is rotating. The results from each combination of modes are compared to the results from the individual modes 1 3 in the absence of mode 4. A range of mode 4 approach velocities are considered since vastly different pressures and changes in viscosity result for each. A single mode 1 angular velocity of fh = 1 rad/s was used, which can be compared with the results from section 5.1.
Table 5.15: Torque on mode 1 rotating sphere at fli = 1^ when combined with mode 4 motion at a range of approach velocities.
Mode 4 Velocity, U4 [m/s\ Tm ode 1 Percent difference from 6/4 = 0
TIsolated
0 (Mode 1 Only) 2.545 N/A
0.5 2.559 0.56
1 2.572 1.07
1.5 2.597 2.08
2 2.691 5.75
The torque on a mode 1 rotating sphere at fh = l1^ for different mode 4 approach velocities is shown in table 5.15. An approach velocity of 0 is included for mode 1 motion in the absence of mode 4 motion. The percent difference column shows how the torque changed from only having mode 1 motion, to the combined motion at the specified approach velocity. We can see the mode 1 torque increasing as the mode 4 approach velocity increases, with a 5% increase in torque at the largest approach velocity of 2.
50


Table 5.16: Force on mode 2 moving sphere at U2 = 1 when combined with mode 4
motion at a range of approach velocities
Mode 4 Velocity, U4 [m/s] FMode 2 Percent difference from U4 = 0
^Stokes
0 (Mode 2 Only) 2.515 N/A
0.5 2.535 0.83
1 2.570 2.21
1.5 2.621 4.22
2 2.757 9.62
The mode 2 force when combined with different mode 4 approach velocities is shown in table 5.16. As in the previous table, an approach velocity of 0 means only mode 2 motion exists. The table shows that an increase in mode 4 approach velocity increases the mode 2 force. At the largest approach velocity of 2, the mode 2 force increases by 9.6% compared to mode 2 spheres in the absence of mode 4 motion.
Table 5.17: Torque on mode 3 rotating sphere at fl3 = 1%^ when combined with mode 4 motion at a range of approach velocities
Mode 4 Velocity, U4 [m/s] TMode 3 Percent difference from 6/4 = 0
TIsolated
0 (Mode 3 Only) 1.052 N/A
0.5 1.052 0.0
1 1.052 0.0
1.5 1.052 0.0
2 1.052 0.0
Table 5.17 shows that the change in mode 3 torque is negligible when mode 4 motion is
combined with mode 3 motion. Although extreme pressures are present in the interstitial
51


region between the spheres, the velocity related to the mode 3 motion is very low on this area of the sphere and there are no shear stresses to influence. The fluid pressure far from the interstitial region is essentially unaffected by the particles approaching and we can say mode 4 motion has no effect on the torque on mode 3 particles.
2.8
X Mode 1 Mode 2
01
=J
D"
CL)
O
2.7
[]
T3
01
in
E
2.6

2.5i

X
X

__I_
0.5
B
1.5
Mode 4 Approach Velocity [m/s]
Figure 5.17: Normalized force or torque for combination of modes 1 and 2 with mode 4. Mode 1 is normalized by the torque on an isolated rotating particle and mode 2 is normalized by the force on a particle in Stokes flow.
The normalized force or torque in tables 5.15 and 5.16 are combined in figure 5.17 by comparing the approach velocity to the normalized force or torque on the rotating or translating sphere. The plot shows that mode 2 has a greater increase in force compared to the increase in torque for mode 1. Mode 3 is not included in the figure because there is no increase in torque due to the addition of the mode 4 motion.
52


CHAPTER VI
RESULTS: PARTICLE ELASTICITY
The work presented so far has considered rigid particles. The particles in suspension systems may have a sufficiently low elasticity and the fluid may be at high enough pressures to result in particle deformation. This section will investigate how particle elasticity affects the individual modes of motion in binary interactions. In particular, we will consider how the deformation affects the hydrodynamic forces and torques on the particles. The simulations are based on Stokes equations for the fluid, coupled with linear elasticity equations for the calculation of elastic deformation.
The analysis will consider four different Youngs moduli of 200 GPa, 20 GPa, 2 GPa, and 0.2 GPa, comparable to those of the common materials, steel, concrete, ABS plastic, and soft rubbers, respectively. A Poissons ratio of 0.4 is used when comparing different elasticities. The last section in this chapter will consider a range of Poissons ratio values for each elasticity. Investigating such a wide range of particle elasticities will help give insight into whether there is a difference between large and small deformations in terms of the hydrodynamic drag force and torque. The hydrodynamic forces and torques were normalized with the hydrodynamic forces and torques for an isolated sphere in the same viscous fluid as the elasticity simulations, allowing us to disregard the units and the type of fluid.
6.1 Rotation About Axis Perpendicular to Line of Centers
Mode 1 requires a three-dimensional geometry but it is possible to describe this
motion with a symmetric model, which reduces the degrees of freedom significantly. The
results in the previous sections regarding pressure dependent viscosity relied on only the
creeping flow module. The spheres were assumed to be rigid and there was no need
to simulate their stresses and strains. Introducing elasticity to the particles requires
the combination of creeping flow and solid mechanics physics. COMSOL has a module
named fluid-structure interaction which allows for the two-way coupling of both of
53


these physics. Applying the same boundary conditions as mode 1 for the rigid spheres and assigning a Youngs modulus, E, to each sphere, we attempted to run the simulation. COMSOL was unable to converge on a result in this three-dimensional simulation, and the problem was decoupled.
If the coupling between physics is removed, we can obtain the fluid pressures adjacent to the spheres surfaces using the creeping flow module. Next, these pressures can be mapped via an interpolation function onto the surfaces of two elastic spheres in a separate solid mechanics simulation. The solid mechanics simulation does not explicitly have the rotating boundary condition applied to the spheres, but rather a pressure held boundary load, as if they were rotating. Once these pressures are mapped to the surfaces of the spheres and the simulation is performed, we have insight into whether the pressures are high enough to significantly deform these spheres for a given elasticity. The geometry and mesh used in the creeping flow portion of the simulation were the same as in the viscosity simulation for mode 1, in section 5.1. The solid mechanics module only requires the meshing of the two spheres since the pressure fields are prescribed from the fluid simulation.
Figure 6.1: Mesh used for the mode 1 solid mechanics module to compute the deformations of the spheres, with 143676 elements and 600459 degrees of freedom.
As shown in figure 6.1, a smaller spherical region of radius rsp^re was created in
54


the center of each sphere and a fixed constraint was defined here. In solid mechanics simulations where external forces are applied, all deformations would be infinite if some degrees of freedom are not fixed in space. Symmetry was applied to the plane bisecting the spheres to reduce the degrees of freedom, similar to the creeping flow simulation. The resulting mesh consists of 143676 elements, 600459 degrees of freedom, and required approximately 9 minutes to solve. The same mesh was used for every separation distance since the pressures in the interstitial region are already determined in a separate simulation and are being mapped onto the elastic spheres surfaces.
Figure 6.2: Exaggerated deformed shape and deformation magnitudes of elastic spheres in mode 1 motion, e = E = 0.2GPa, v = 0.4, scaled to 10% of the total geometry.
The resulting deformation from mode 1 is shown in figure 6.2, where we see a compression left of the line of centers and an expansion to the right of the line of centers. The deformed shape is exaggerated because actual deformation is not visible to the eye. In the simulation, the lower sphere is stationary and the upper sphere is rotating counterclockwise. The upper sphere is dragging fluid into the gap and compressing into the wedge shape formed by the spheres curvature. This results in raised pressures and a compression of the surfaces to the left of the line of centers. The spheres surfaces to the right of the line of centers are separating and causing a vacuum which contracts the surfaces of the sphere inward.
55


The stationary sphere experiences a larger deformation magnitude because it only has a large rightward force in the interstitial region. There is no force on the opposite side of the sphere, counteracting this rightward force. In contrast, the upper sphere has a leftward force in the interstitial region, but also has a rightward force on its opposite side, from the rotational resistance. We can see the upper sphere stays relatively centered in its original position, while all points in the lower sphere have a net rightward shift, which causes the large deformation magnitude. This can be summarized be the upper sphere having more balanced shear stresses on its surface than the lower sphere.
Table 6.1: Maximum particle deformations for mode 1, corresponding to fluid pressures obtained from rigid sphere simulation, normalized by sphere radius.
Separation, e Max. Deformation / Particle Radius
E = 200 GPa E = 20 GPa E = 2 GPa E = 0.2 GPa
r/10 2.42e-12 2.42e-ll 2.42e-10 2.42e-9
r/100 1.34e-ll 1.34e-10 1.34e-9 1.34e-8
r/1000 1.03e-10 1.03e-9 1.03e-8 1.03e-7
r/10000 1.13e-9 1.13e-8 1.13e-7 1.13e-6
The maximum deformation of the rotating sphere is shown in table 6.1, where we see that even at the smallest separation distance and lowest elastic modulus, a normalized deformation of only 1.13e-6 is experienced. This value is so insignificant that we dont expect a change in hydrodynamic forces compared to rigid spheres. Since the two different physics are decoupled, the pressures obtained are for rigid spheres which we expect to be larger than the pressures on elastic spheres. If there are high fluid pressures in the interstitial region, the gap will grow as the surfaces deform and the pressures will decrease. The coupling of these solutions would likely show a decrease in deformation, compared to the results in table 6.1.
56


6.2 Translation Along Axis Perpendicular to Line of Centers
Translation perpendicular to the line of centers, or mode 2, introduces a shearing motion to the surfaces of the sphere. The upper sphere was given a translation velocity of hm/s along the axis perpendicular to the line of centers, with the lower sphere stationary. The same geometry and mesh was used as in the mode 1 elasticity simulations in section 6.1 and the simulations were performed with the same decoupled method.
Figure 6.3: Exaggerated deformed shape of elastic spheres in mode 2 motion, e = pyr, E = 0.2GPa, v = 0.4, scaled to 10% of the total geometry. Deformations are normalized by sphere radius.
The deformed shape of the spheres are shown in figure 6.3, where the upper sphere is moving toward the right and the lower sphere is stationary. The black circles represent the original shape of the spheres. From the plot, we can see that the upper sphere has a greater deformation in the direction opposite of motion than the lower, stationary sphere. This is because the Stokes drag on the moving sphere is added to the shear stress from the pressure generated in the interstitial region. There is also an upward compression of the surface of the moving sphere in the interstitial region, from the higher pressures in the fluid. The lower sphere has a small amount of deformation, caused by the shearing effect from the moving sphere, but there is no Stokes drag on this sphere so the deformation is lower.
57


Table 6.2: Maximum particle deformations for mode 2 at several separation distances
and particle elasticities, corresponding to fluid pressures obtained from rigid sphere sim-
ulation, normalized by sphere radius.
Separation, e Max. Deformation / Particle Radius
E = 200 GPa E = 20 GPa E = 2 GPa E = 0.2 GPa
r/10 3.75e-12 3.75e-ll 3.75e-10 3.75e-9
r/100 4.08e-12 4.08e-ll 4.08e-10 4.08e-9
r/1000 4.22e-12 4.22e-ll 4.22e-10 4.22e-9
r/10000 4.30e-12 4.30e-ll 4.30e-10 4.30e-9
The maximum deformation of the particles with different elasticities are shown in table 6.2 at several separation distances. At all separation distances, as the elasticity is decreased by an order of magnitude, the deformation also decreases by an order of magnitude. The maximum deformation was found to be independent of the separation distance for this motion. The surface deformations are so insignificant that we can expect the hydrodynamic forces to be approximately equal to those of rigid spheres for this mode.
6.3 Rotation About Line of Centers
Binary spheres rotating about their line of centers is a motion that we can utilize axisymmetric geometry for, since COMSOL allows a velocity component in the radial, or cf>, direction. This motion was also solved with a decoupled approach, where the creeping flow module was used to simulate rigid spheres rotating in mineral oil. The pressures from this simulation were used in a separate solid mechanics simulation of elastic spheres. The geometry and mesh for the creeping flow simulations used to obtain fluid pressures were identical to the mesh for mode 3 pressure dependent viscosity simulations. The two-dimensional solid mechanics mesh is composed of 6522 elements with 26628 degrees of freedom. The pressure dependent viscosity simulations for mode 3 in chapter
5 determined that insufficient pressures are generated from this motion to raise the fluid
58


viscosity. The motion has no singularity and we can expect to see similar results in the elasticity simulations.
Table 6.3: Maximum particle deformations for mode 3 at several separation distances and particle elasticities, corresponding to fluid pressures obtained from rigid sphere simulation, normalized by sphere radius.
Separation, e Max. Deformation / Particle Radius
E = 200 GPa E = 20 GPa E = 2 GPa E = 0.2 GPa
r/10 0 0 0 0
r/100 0 0 0 0
r/1000 0 0 0 0
r/10000 0 0 0 0
Table 6.3 shows that elastic spheres in mode 3 motion experience no deformation. This is because the gauge pressure at every point on the spheres surface is negligible for this motion. Since there is no deformation, we can expect hydrodynamic forces that are equal to those of rigid spheres. The motion for rotation about the particles line of centers has no singularity, and the only deformation would be in direction tangent to the surface, which would not cause an increase in the hydrodynamic torque for a steady state solution.
6.4 Translation Along Line of Centers
Spheres translating along their line of centers can be solved with a two-dimensional axisymmetric geometry. A geometry similar to the mode 4 pressure dependent viscosity simulations in section 5.4 was used, except the problem is fully coupled and includes a domain for the sphere, rather than only its boundary. The other elasticity simulations required the decoupling of the creeping flow and solid mechanics physics, but having the two-dimensional geometry and a simple motion allowed for a successful solution of the
59


fully coupled problem. In this case, the results from the creeping flow simulation are used as inputs to the solid mechanics module, which passes its results back to the creeping flow module for another iteration. This process is repeated until the solution reaches a steady state, at which time the results reflect a fully-coupled simulation.
Figure 6.4: Mesh used for a fully-coupled simulation of mode 4 with elastic spheres.
A separate mesh is constructed for the creeping flow and solid mechanics regions, as seen in figure 6.4. First the fluid mesh was constructed identical to the mesh in the mode 4 pressure dependent viscosity simulations in section 5.4. Next, the mesh for the elastic sphere was constructed. The adjacent meshes share nodes along their boundaries, so the bottom of the sphere where it meets the fluid has a dense mesh and the elements farther from the bottom grow in size. The greatest amount of stress and deformation is where the fluid pressures are high in the interstitial region, so a dense mesh for the sphere is necessary here. There were 66231 triangular elements in the mesh for the sphere, in addition to the elements in the fluid region. Simulations for mode 4 with elastic particles were performed for spheres of a material with a Poissons ratio of 0.4 at a range of particle elasticities and separation distances.
60


Figure 6.5: Exaggerated deformation of mode 4 elastic particles of E=0.2 GPa, e = ^ in a highly viscous fluid with /i = 0.5, normalized by sphere radius.
The area of the highest deformation is shown in figure 6.5, in the hemisphere closest to the interstitial region. The area of highest pressures in the fluid acts on the adjacent sphere boundary, causing straining on this area of the sphere. The black line in the figure shows the original curvature of the sphere.
Table 6.4: Maximum normalized mode 4 particle deformations at several separation distances and particle elasticities in a fluid with p, = 0.5 Pa s. The simulation is fully-coupled between solid mechanics and creeping flow.
Separation, e M ax .Deformati on / Particle Radius
E = 200 GPa E = 20 GPa E = 2 GPa E = 0.2 GPa
r/10 4.56e-9 4.56e-8 4.56e-7 4.56e-6
r/100 4.32e-8 4.31e-7 4.31e-6 4.28e-5
r/1000 6.41e-7 6.33e-6 5.73e-5 3.58e-4
r/10000 l.lle-5 5.73e-5 N/A N/A
The maximum deformations for elastic particles in a highly viscous fluid of p. = 0.5 Pa s are shown in table 6.4. The absent values for the lower elasticity particles are a result of large deformations on the particles surface. The numerical solver was unable to converge for the cases when deformations were very large. The largest normalized
61


deformation of 3.58e-4 is observed at e = with a particle elasticity of E = 0.2 GPa. At the larger separation distances, the deformation increases by an order of magnitude when the elasticity decreases by an order of magnitude, but this trend abates as the separation distance decreases.
03
Q.
(D
=3
CO
03
<13
0.
=3
03
(3
X
03
-4 -3 -2 -1 0
10 10 10 10 10
Separation Distance / Radius
Figure 6.6: Maximum fluid pressure in the interstitial region for mode 4 motion of elastic spheres in mineral oil, at a range of particle separation distances. E = 20 GPa, v = 0.4.
As the spheres become closer together, the fluid pressure between them increases rapidly, which is directly related to the force on the particle and the particles deformation. An order of magnitude separation distance decrease does not result in an order of magnitude pressure increase at small separation distances, as we see in figure 6.6. The figure shows the maximum pressure in the interstitial region for particles of E = 20 GPa at a range of separation distances. The relationship between pressure and both particle deformation and force is nonlinear at small separation distances.
x 10
62


Table 6.5: Normalized forces on mode 4 elastic binary particles at a range of separation distances and particle elasticities.
Separation, e FElastic/ FStokes
E = 200 GPa E = 20 GPa E = 2 GPa E = 0.2 GPa Rigid
r/10 7.406 7.404 7.403 7.403 7.413
r/100 53.019 53.016 52.983 52.657 53.422
r/1000 503.460 499.290 464.414 324.559 504.434
r/10000 4370.819 2792.957 N/A N/A 5004.90
We can see in table 6.5 that the force on elastic particles decreases significantly from rigid particles. This is most evident at the lowest elasticity value of E = 0.2 GPa, which we can compare with the force on rigid spheres.
Table 6.6: Percent decrease in mode 4 force on elastic particles, compared to rigid particles.
Separation, e Percent decrease from rigid particle force
E = 200GPa E = 20 GPa E=2GPa E = 0.2 GPa
r/10 0.10 0.13 0.13 0.14
r/100 0.76 0.76 0.82 1.43
r/1000 0.19 1.02 7.93 35.66
r/10000 14.51 44.20 N/A N/A
Table 6.6 shows that the smallest separation distances have the greatest decrease in drag force. This is because the spheres have the largest amount of surface deformation at these close distances as a result of the extreme fluid pressures.
63


Figure 6.7: Mode 4 hydrodynamic drag forces with elastic particles, where the separation distance does not account for surface deformation.
It appears from figure 6.7 that the particles with lower elasticity have a significantly lower drag force at small separation distances, when the pressure is the highest. This is because as the surfaces of the spheres deform, the gap grows and the pressure between the spheres is relieved. The change in separation distance from the particle deformation is unaccounted for in this figure.
Figure 6.8: Mode 4 hydrodynamic drag forces with elastic particles. Separation distance has been adjusted to account for the surface deformation.
We can see from figure 6.8 that if the true surface separation distance is plotted
64


against the drag force, the results essentially collapse back onto the curve for rigid spheres. This is done by adjusting the separation distance to account for the surface deformation of both spheres. This leads to the conclusion that the drag force is a function of the separation distance of the spheres only, and not of the particle elasticity.
6.5 Poissons Ratio
When a material is compressed along an axis, it tends to expand in the other two orthogonal directions. The Poissons ratio of an isotropic, linear elastic material describes the relationship between transverse and axial strain, by how much the material expands relative to the deformation in the axis of compression or tension. Materials that demonstrate a significant transverse displacement when compressed axially, such as rubbers, have a Poissons ratio close to 0.5 and are classified as nearly incompressible. Materials that experience a negligible transverse displacement compared to the axial displacement, such as cork, have a Poissons ratio close to zero. Studying the effects of the modulus of elasticity is directly related to Poissons ratio and we will expand upon the investigation for mode 4 where the particle deformations were the largest of the modes.
Table 6.7: Normalized hydrodynamic forces at e = and a range of elasticities and Poissons ratios.
Poissons ratio, v F /FStokes
E = 200 GPa E = 20 GPa E = 2 GPa E = 0.2 GPa
0.25 503.418 498.887 461.539 317.917
0.3 503.428 498.993 462.281 319.583
0.35 503.439 499.121 463.226 321.758
0.4 503.460 499.290 464.414 324.559
0.45 503.481 499.502 465.963 328.379
0.49 503.513 499.778 468.096 333.992
Table 6.7 shows how the different Poissons ratio values affect the hydrodynamic
65


forces at each modulus of elasticity. We can see that as the Poissons ratio increases, the drag force on the particles increases. The difference in drag force between v = 0.25 and 0.49 for an elasticity of 0.2 GPa is approximately 5%, which is a small difference when comparing to the effect of elastic modulus on particle drag force before accounting for the surface deformation. Although this is by no means a direct comparison, we can say that the Poissons ratio of the material has little to no effect on the drag force, at least at the separation distances e = A smaller separation distance was not simulated because the solver would not converge for the combination of separations below e = ^ and elasticities below E = 20 GPa.
Poissons Ratio, v
Figure 6.9: Comparison of Poissons ratios vs. hydrodynamic force for E = 0.2 GPa,
_ rad
t 103 *
A range of Poissons ratios are plotted against the normalized force on particles with E = 0.2 GPa in mode 4 motion at e = in figure 6.9. The plot shows that the greater the Poissons ratio, the greater the corresponding drag force is.
66


Table 6.8: Percent decrease in mode 4 force on elastic particles, compared to rigid par-
ticles.
Poissons ratio, v Percent decrease from rigid particle force
E = 200 GPa E = 20 GPa E = 2 GPa E = 0.2 GPa
0.25 0.20 1.10 8.50 36.98
0.3 0.20 1.08 8.36 36.65
0.35 0.20 1.05 8.17 36.21
0.4 0.19 1.02 7.93 35.66
0.45 0.19 0.98 7.63 34.90
0.49 0.18 0.92 7.20 33.79
The percent decreases in drag force compared to a rigid particle are shown in table 6.8 for a range of particle elasticity and Poissons ratio values. For each value of E, the decrease in force is less for the higher Poissons ratio materials. This implies that the more compressible a material is, the more the drag force will be affected by the particles change in shape.
Table 6.9: Particle deformations at e = for a range of elasticities and Poissons ratios.
Poissons ratio, v Maximum surface deformation / radius
E = 200 GPa E = 20 GPa E = 2 GPa E = 0.2 GPa
0.25 6.55e-7 6.47e-6 5.85e-5 3.63e-4
0.3 6.41e-7 6.33e-6 5.73e-5 3.58e-4
0.35 6.23e-7 6.15e-6 5.58e-5 3.51e-4
0.4 6.00e-7 5.94e-6 5.40e-5 3.42e-4
0.45 5.72e-7 5.66e-6 5.17e-5 3.31e-4
0.49 5.37e-7 5.32e-6 4.88e-5 3.17e-4
67


Table 6.9 contains the data we can use to see if the Poissons ratio of the particles has the same effect as the elasticity. The maximum deformation for each combination of E and v are shown in the table, where we see the greatest deformation for the lowest E and lowest v. This combination of material properties correspond to the softest, most compressible materials.
Deformation unaccounted for Deformation accounted for
Figure 6.10: Comparison of Poissons ratios of 0.25 and 0.49 for several elasticities. Left: deformation unaccounted for in separation, right: deformation accounted for.
The plots for v = 0.25 and v = 0.49 are shown for elasticities of 20 GPa, 2 GPa, and 0.2 GPa, in figure 6.10. On the left plot, we can see that the force appears lower for the elastic materials, which is what was seen in section 6.4. It also appears that the lower the Poissons ratio, the lower the force. The particle deformation is accounted for in the separation distance on the right plot and we see that the forces for every combination of E and v agree with the force on a rigid particle. Particles with a low Poissons ratio of 0.25 have a larger deformation than those with a Poissons ratio of 0.49 and the point has a larger rightward shift than the higher Poissons ratio of 0.49 on the right plot.
68


Deformation unaccounted for Deformation accounted for
Figure 6.11: Comparison of v = 0.25 and v = 0.49 for E = 0.2 GPa. Left: deformation unaccounted for in separation, right: deformation accounted for.
The forces on particles with Poissons ratios of 0.25 and 0.49 are shown in figure 6.11 for a particle elasticity of 0.2 GPa. The left plot shows that as the separation distance decreases, the lower Poissons ratio particles have a slightly lower drag force. The less compressible the material is, for instance a material with a Poissons ratio of 0.49, the higher the drag force appears to be. This is because the deformation of the particles surface is greater since the material is more compressible, giving more space for the fluid in the interstitial region to relieve itself of high pressures. In the right plot, we adjust the separation distance by the particles surface deformation and the curves essentially collapse onto one another, which is similar to the result in the elasticity simulations in the previous sections. This behavior is consistent with the observation that the drag force is a function of the actual separation distance between the particles surfaces, and not of their Poissons ratio.
69


CHAPTER VII CONCLUSION
Pressure-dependent fluid viscosity and particle elasticity both have unique affects on the hydrodynamic forces and torques on binary spheres in suspension systems. Fluid viscosity can increase significantly at enhanced pressures, which can affect the motion of binary spheres in viscous fluid. We saw that binary spheres translating along their line of centers have the largest increase in hydrodynamic force due to pressure-dependent viscosity, compared to the other modes of motion. An approximate 39% increase in hydrodynamic force was observed for particles in this motion at a velocity of 2.15 in mineral oil. The majority of the force is from the extreme fluid pressures in the interstitial region between the spheres. The other three modes of motion did not significantly increase the hydrodynamic force or torque when introducing pressure-dependent viscosity because the pressures generated are insufficient to increase the fluids viscosity, significantly.
The actual motion of binary particles in suspensions is a combination of the four individual motions, so mode 4 was combined with the other three modes to determine if the presence of large pressures affected the hydrodynamic force or torque in modes 1-3. The mode 4 particle velocity was 2 in these simulations. The increase in torque on mode 1 particles rotating at 1 ^ was 5.75% and the force increase on mode 2 particles translating at 1 was 9.62%. Mode 3 particles rotating at 1 ^ experienced no torque increase when combined with mode 4 motion because the mode 3 components of velocity are negligible in the interstitial region where the high pressures are present. The pressure-dependent viscosity simulations show that if particle velocities are sufficient, the hydrodynamic forces imposed can be significantly perturbed by a fluid with pressure-dependent viscosity, compared to constant viscosity.
Particle elasticity was simulated for binary spherical particles in viscous fluids with
each of the four modes of motion. Particles in mode 4 motion were found to have the
70


largest deformations compared to the other modes. A significant decrease in drag force was observed for mode 4 elastic particles, compared to rigid particles, at comparable center-to-center particle distances. The largest force decrease of 44% in mode 4 was observed at a separation distance of and a particle elasticity of 20 GPa. The simulations for particle elasticities below 20 GPa did not converge due to large mesh deformations. Although a 44% decrease from the force on rigid particles appears to be significant, a shift of the separation distance by the amount of particle deformation shows otherwise. After measuring the surface deformation of the particles and accounting for this in the separation distance, the gap size increases and these points shift onto the original curve for a rigid particle. These results determine that particle elasticity does not have a direct affect on the hydrodynamic force because elastic and rigid particles incur the same forces when the elastic particle separation is adjusted to account for deformation.
The simulations used to investigate the affect of particle elasticity assumed a Poissons ratio value of 0.4. Because of its relationship with the modulus of elasticity, the Poissons ratio of the particles was also studied. A range of Poissons ratios from v = 0.25 through v = 0.49 were considered for each elastic modulus value. The forces on particles did not change significantly at different values of v for each value of E. Similar to the results for particle elasticity, if the separation distance is adjusted to account for the surface deformation, the forces essentially fall back on the curve for a rigid particle. We can conclude that the Poissons ratio of a material does not have a significant affect on the hydrodynamic forces or torques on binary particles in viscous fluids.
71


REFERENCES
[1] E.P. Ascoli, D.S. Dandy, L.G. Leal, Low Reynolds Number Hydrodynamic Interaction of a Solid Particle With a Planar Wall. J. Numer. Meth. Fluids 9, 651-688.
[2] S. Bair, High-Pressure Rheology for Quantitative Elastohydrodynamics, Elsevier, Amsterdam, The Netherlands, 2007.
[3] C. Barus, Isothermals, Isopiestics and Isometrics Relative to Viscosity, Am. J. Science, Third Series, Vol. XLV, No. 266, 1893.
[4] A.W. Batchelor, G.W. Stachowiak, Engineering Tribology, second edition, Butterworth-Heinemann, Oxford, UK, 2001.
[5] G.K. Batchelor, J.T. Green, The determination of the bulk stress in a suspension of spherical particles to order e?, J. Fluid Mech., 56 (1972), 401-427.
[6] R.B. Bird, E.N. Lightfoot, W.E. Stewart, Transport Phenomena, second edition, John Wiley & Sons, New York, N.Y., 2001.
[7] H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface, (1961). Chem. Eng. Sci. 16, 242-251.
[8] H. Brenner, J. Happel, Low Reynolds number hydrodynamics, Springer, The Netherlands, 1983.
[9] S.R. Challa, F.V. Swol, Molecular simulations of lubrication and solvation forces, Phys. Rev. E 73 (2006), 016306.
[10] J.S. Chong, E.B. Christiansen, A.D. Baer, Rheology of Concentrated Suspensions, J. Applied Polymer Sci., 15 (1971), 2997-2021.
[11] G. DAvino, F. Snijkers, R. Pasquino, et al., Migration of a sphere suspended in viscoelastic liquids in Couette flow: experiments and simulations, Rheologica Acta (2012) 51:215.
[12] A. Einstein, Eine neue Bestimmung der Molekldimensionen, (1906). Annalen der Physik. 19 (2): 289.
[13] G.E. Fernandez, P.M. Carrica, D.A. Drew (2000) Binary Interaction Force on Spheres I: Low Reynolds Number Flow, Chem. Eng. Comm.
[14] E. Guth, R. Simha, Untersuchungen ber die Viskositt von Suspensionen und Lsun-gen. 3. ber die Viskositt von Kugelsuspensionen, (1936). Kolloid Z. 74 (3): 266.
[15] W. Habchi. A Full-System Finite Element Approach to Elastohydrodynamic Lubrication Problems : Application to Ultra-Low-Viscosity Fluids. PhD thesis, INSA de Lyon, Villeurbanne, Prance; 2008.
72


[16] J.S. Halow, G.B. Wills, Experimental Observations of Sphere Migration in Cou-ette Systems, (1970). Ind. Eng. Chem. Fundamen. 9 (4): 603-607
[17] B.J. Hamrock, Fundamentals of Fluid Film Lubrication, Nasa Reference Publication 1255 (1991), 55-59.
[18] M.S. Ingber, S. Feng, A.L. Graham, H. Brenner, The analysis of self-diffusion and migration of rough spheres in nonlinear shear flow using a traction-corrected boundary element method, J. Fluid Mech. 598 (2008), 267-292.
[19] D.J. Jeffrey, Y. Onishi, The forces and couples acting on two nearly touching spheres in low-Reynolds-number flow, Z. Angevj. Math. Phys. 35 (1984), 634-641.
[20] I.M. Krieger, T.J. Dougherty, A mechanism for non-Newtonian flow in suspensions of rigid spheres., Trans. Soc. Rheol. Ill (1959), 137-152.
[21] S.R. Majumdar, M.E. ONeill, Asymmetrical slow viscous fluid motions caused by the rotation or translation of two spheres. Part I. The determination of exact solutions for any values of the ratio of radii and separation parameters, Z. Angevj. Math. Phys. 21 (1970), 164-179.
[22] O. Reynolds, On the theory of lubrication and its application to Mr. Beauchamp Towers experiments, including an experimental determination of the viscosity of olive oil, Phil. Trans. Royal Soc., 177 (1886), 157-234.
[23] C.J.A. Roelands, Correlational Aspects of the Viscosity-Temperature-Pressure Relationship of Lubricating Oils, Ph.D. Thesis, University of Technology, Delft, 1966.
[24] H. Sohn, S. Koo, Permeability of Viscous Flow Through Packed Bed of Bidisperse Hard Spheres, (2012). Korean Chemical Engineering Research. 50 (1): 66-71.
[25] N. Tetlow, A.L. Graham, M.S. Ingber, S.R. Subia, A.L. Mondy, S.A. Altobelli, 1998. Particle migration in a Couette apparatus: Experiment and modeling. J. Rheol. 42, 307-327.
[26] D.G. Thomas, Transport characteristics of suspension: VIII. A note on the viscosity of Newtonian suspensions of uniform spherical particles, (1965). J. Colloid Sci. 20 (3): 267
73


APPENDIX
A. Mesh Convergence Studies
Convergence studies are important in numerical simulations to guarantee the mesh construction is adequate for the numerical simulation to produce reliable results for the physics in the simulation. A convergence study does not guarantee that the results are correct, but rather that convergence is being obtained and the results between successive simulations of varying (but similar) meshes are not changing more than a predetermined threshold. An example is for the deformation of a beam, where an extremely coarse mesh is first defined. A simulation using this mesh will yield an initial result which we are not sure is converged or not. The mesh is then made slightly denser and the new result might be much different than the previous result.
This iterative process is continued until the mesh is fine enough where the difference between solutions is below approximately 1%. At this point, we have determined the mesh is converged below the threshold. A convergence study is provided for each simulation performed in this thesis, at the smallest separation distance of e = ^|.
Viscosity Simulations
Mode 1 (Rotation about axis perpendicular to line of centers)
Table A.l: Convergence study for mode 1 three-dimensional mesh at e = based on change in force between subsequent meshes.
No. elem. Deg. of freedom Torque [N-m] Percent T change
264996 2890202 4.104 N/A
298614 3267810 4.249 3.53
315246 3449502 4.326 1.81
354150 3875544 4.357 0.71
74


Increasing the total number of elements in the mesh allowed convergence to be obtained at 354150 elements, with only a 0.71% change from the previous mesh.
Mode 2 (Translation along axis perpendicular to line of centers)
Table A.2: Convergence study for mode 2 three-dimensional mesh at e = ^ based on change in force between subsequent meshes.
No. elem. Deg. of freedom Force [N] Percent F change
264996 2890202 2.371 N/A
298614 3267810 2.449 3.18
315246 3449502 2.498 1.96
354150 3875544 2.515 0.68
The mode 2 mesh convergence is similar to the mode 2 results. Convergence was reached at the same amount of elements for both modes.
Mode 3 (Rotation about line of centers)
Table A.3: Convergence study for mode 3 two-dimensional mesh at e = ^ based on change in force between subsequent meshes.
No. elem. Deg. of freedom Torque [N-m] Percent T change
4586 32522 1.941 N/A
16551 112643 1.801 6.96
36962 201664 1.793 0.55
44862 210632 1.791 0.11
Mode 4 (Translation along line of centers)
In mode 4, the majority of the force is from the high fluid pressures in the interstitial region. The bulk fluid mesh can be relatively coarse, with a high density mesh in the interstitial region.
75


Table A.4: Convergence study for mode 4 two-dimensional mesh at e = based on
change in force between subsequent meshes.
No. elem. Deg. of freedom Force [N] Percent F change
30195 141010 13956 N/A
58665 274553 13837 0.85
62143 290344 13814 0.16
65284 304581 13814 0.00
A large difference is not seen when increasing the number of elements in the mode 4 meshed because the interstitial elements are guaranteed to be fine enough to fit in the gap. If this condition was satisfied, the mesh was satisfactory since most of the force came form the interstitial region pressures.
Elasticity Simulations Mode 1 and 2
The geometry and mesh was identical for mode 1 and 2 elasticity simulations. Once the simulation was performed for the fluid portion of the problem, the resulting pressures were mapped onto elastic spheres in another solid mechanics simulation. The mesh for the fluid portion was identical to the mesh for the mode 1 and 2 constant viscosity numerical simulations. The same solid mechanics mesh was used in modes 1 and 2, based on Mode 1 convergence results.
76


Table A.5: Convergence study for mode 1 three-dimensional solid mechanics mesh at
e = and E = 0.2 GPa based on change in deformation between subsequent meshes.
No. elem. Deg. of freedom DcfyVorm. Percent def. change
26442 112537 1.03e-6 N/A
31468 135426 1.09e-6 5.82
76998 331168 1.12e-6 2.75
143676 600459 1.13e-6 0.89
The density of the mesh in the interstitial region is important to determine the correct displacement for modes 1 and 2. The pressures in this area are large, so mapping the pressures onto a denser mesh in this area is important. A relatively uniform mesh was used for the whole sphere, so the amount of elements is large to guarantee a high enough density in the interstitial region.
Mode 3
Table A.6: Convergence study for mode 3 two-dimensional solid mechanics mesh at e = and E = 0.2 GPa based on change in deformation between subsequent meshes.
No. elem. Deg. of freedom Dcf/Vorm. Percent def. change
976 4188 0 N/A
1480 6204 0 0
6522 26628 0 0
Mode 3 experiences no deformation since fluid pressures from the creeping flow simulations are very low, so the density of the mesh was not critical.
Mode 4
Mode 4 used a fully coupled approach where the mesh of the solid mechanics and creeping flow modules were adjacent in the same simulation.
77


Table A.7: Convergence study for mode 4 two-dimensional fully-coupled simulation mesh
at e = Yp and E = 20 GPa based on change in force between subsequent meshes.
No. solid elem. Total deg. of freedom F Percent F change
FStokes
41248 235938 2652.2 N/A
58912 338744 2720.6 2.59
62985 359014 2768.10 1.75
66231 380184 2792.95 0.89
The density of the fluid mesh is limited by the density of the fluid mesh which was constructed first. The solid mesh begins with the size of elements from the adjacent fluid mesh and expands them upward away from the interstitial region. The mesh is not able to grow too coarse because of this, so the forces for all meshes are not far from the converged mesh.
78


Full Text

PAGE 1

EFFECTSOFELASTICITYANDPRESSUREDEPENDENTVISCOSITYONTHE RELATIVEMOTIONOFBINARYSPHEREINTERACTIONSINHIGHLY VISCOUSFLUIDFLOWS by ADAMHWOODMAN B.S.,UniversityofNewHampshire,2013 Athesissubmittedtothe FacultyoftheGraduateSchoolofthe UniversityofColoradoinpartialfulllment oftherequirementsforthedegreeof MasterofScience MechanicalEngineering 2016

PAGE 2

ThisthesisfortheMasterofSciencedegreeby AdamHWoodman hasbeenapprovedforthe DepartmentofMechanicalEngineering by MarcIngber,Chair MaryamDarbeheshti,Advisor DanaCarpenter AlanGraham December17,2016 ii

PAGE 3

Woodman,AdamHM.S.MechanicalEngineering EectsofElasticityandPressureDependentViscosityontheRelativeMotionofBinary SphereInteractionsinHighlyViscousFluidFlows ThesisdirectedbyAssistantProfessorMaryamDarbeheshti ABSTRACT Theeectofpressure-dependentuidviscosityandparticleelasticityinisothermal suspensionsofsphericalparticlesisinvestigatedusingnumericalsimulationsoftware. Therelativemotionoftwointeractingparticlessuspendedinviscousshearowscan bedecoupledintofourmodesofmotion.Thesemodesincludetranslationalongtheir lineofcenters,translationperpendiculartotheirlineofcenters,rotationabouttheir lineofcenters,androtationaboutanaxisperpendiculartotheirlineofcenters.Exact solutionstoStoke'sequationsexistforallfourmodesbutarelimitedtorigidparticles inauidwithconstantviscositythroughouttheoweld.Inreality,theseparticlesare deformableandtheviscosityofthesurroundinguidmayincreasesignicantlyinhigh pressureregions. Inparticular,highpressureregionscanariseintheinterstitialregionbetweenthe particleswhentheyaresucientlyclose.Inthisresearch,viscousforcesareassumedto dominatetheinertialforcessothatStoke'sequationsapply.Usingniteelementsoftware, theaectofpressure-dependentviscosityandparticleelasticityisstudiedforeachmode ofmotionatarangeofparticleseparationdistances.Thisdataallowsthequantication ofhowpressure-dependentviscosityandparticleelasticityaectthehydrodynamicforces onbinaryparticles. Introducingpressure-dependentviscositysubstantiallyincreasedtheforceontheparticleswhentheytranslatealongtheirlineofcenters,comparedtoaconstantviscosity uid.Pressureincreaseswerenotsucientfortheviscositytoincreasesignicantlyin iii

PAGE 4

allothermodesofmotion,andthenumericalpressure-dependentandconstantviscositysimulationsproducedessentiallythesameresults.Theeectofpressure-dependent viscosityismostprevalentwhentwoparticlesaretranslatingalongtheirlineofcenters becausethehighestpressuresaregeneratedinthismotion.Theuidviscosityincreased signicantlyinthismotionintheinterstitialregionbetweentheparticles,leadingtoa largeincreaseintheforce.Whenbinaryparticlesrotateaboutanaxisperpendicular totheirlineofcenters,translateaboutanaxisperpendiculartotheirlineofcenters,or rotateaboutthelineofcenters,theuidpressureintheinterstitialregionisinsucient toincreasetheforcesortorquesontheparticlessignicantly. Particleelasticitywasfoundtohaveaninsignicantaectonthehydrodynamic forcesontheparticles.Theactualsurfaceseparationdistancebetweenparticlesisthe determiningfactorfortheuidpressureandforcesontheparticles.Whiletheparticles deforminthecaseofparticlestranslatingalongtheirlineofcenters,thedeformation causesthesurfaceseparationdistancetoincrease.Thecenter-to-centerdistanceofthe particleshasnotchanged,soitappearsthattheuidpressureandhydrodynamicforces decreasemoreforsofterparticles.Whenthehydrodynamicforcesarecomparedtothe actualsurfaceseparation,theforcesarefoundtobeapproximatelyequaltotheforces onrigidparticles.Comparingdierentvaluesoftheparticles'Poisson'sratio,asimilar aecttoelasticityontheforceswasobserved. Theformandcontentofthisabstractareapproved.Irecommenditspublication. Approved:MaryamDarbeheshti iv

PAGE 5

DEDICATION Thisthesisisdedicatedtomyparents,whohavesupportedeverythingIhavepursued. Theyhavetaughtmewhathardworkis,allowedmetondmyinterests,andhavebeen mymainsourceofencouragementandinspirationthroughoutmylife.NothingIhave accomplishedwouldhavebeenpossiblewithouttheirendlessadvice,love,support,and dedicationtobeingthebestrolemodelspossible. v

PAGE 6

ACKNOWLEDGMENTS Thisresearchwouldnothavebeenpossiblewithouttheopportunitiesthatmyadvisor,Dr.MaryamDarbeheshti,haspresentedmewith.Herknowledge,support,kindness, andpositiveattitudewerekeyfactorsinmygraduatecareer. IalsowanttothankDr.IngberandDr.Grahamfortheirintroductiontothis research.Theirendlessknowledgeonthetopicmadeitpossibleformetocontinueafter hittingroadblocks.Theyaskedtherightquestionsandpointedmeinthedirectionof severalsourcesofinformationthatwereessentialtothisthesis. vi

PAGE 7

TABLEOFCONTENTS CHAPTER I.INTRODUCTION.................................1 1.1GoverningEquations............................2 1.2BinaryParticleInteractions.........................4 1.3LiteratureReview..............................4 II.ANALYTICSOLUTIONS.............................9 2.1RelatedBenchmarkSolutions........................9 2.2SolutionstoFourModesofMotion.....................10 III.NUMERICALSIMULATIONS..........................16 3.1MeshDiscretization.............................21 3.2BoundaryConditions............................23 IV.FLUIDCONSTITUTIVEEQUATIONS.....................25 4.1ViscosityModels...............................25 V.RESULTS:PRESSUREDEPENDENTVISCOSITY..............28 5.1RotationAboutAxisPerpendiculartoLineofCenters..........28 5.2TranslationAlongAxisPerpendiculartoLineofCenters.........34 5.3RotationAboutLineofCenters......................39 5.4TranslationAlongLineofCenters.....................43 5.5CombinedMotions..............................48 VI.RESULTS:PARTICLEELASTICITY......................53 6.1RotationAboutAxisPerpendiculartoLineofCenters..........53 6.2TranslationAlongAxisPerpendiculartoLineofCenters.........57 6.3RotationAboutLineofCenters......................58 6.4TranslationAlongLineofCenters.....................59 6.5Poisson'sRatio................................65 VII.CONCLUSION...................................70 vii

PAGE 8

REFERENCES.....................................72 APPENDIX A.MeshConvergenceStudies.............................74 viii

PAGE 9

FIGURES FIGURE 1.1Fourmodesofmotionpossibleforbinaryparticleinteractions.........4 2.1Torqueonbinaryspheres,onerotatingaboutanaxisperpendiculartotheir lineofcentersandonestationary,normalizedbythetorqueonanisolated rotatingsphere...................................11 2.2Forceonbinaryspheres,onetranslatingalongaxisperpendiculartolineof centersandonestationary,normalizedby F Stokes ................13 2.3Torqueonbinaryspheres,onerotatingaboutlineofcentersandonestationary,normalizedbythetorqueonanisolatedrotatingsphere..........14 2.4Forceonbinaryspherestranslatingalongtheirlineofcenters,normalizedby F Stokes .......................................15 3.1Comparisonoftheanalyticandnumericalsolutionsforarigidsphereapproachinganinniteplanewall..........................19 3.2Exampleofmeshforbinaryspherestravelingalongtheirlineofcenters...22 3.3Detailedviewofinterstitialregionmeshforbinaryspheres...........22 4.1ComparisonbetweenBarusandRoelandsmodelsforviscosityofmineraloil.26 4.2Roeland'spressure-viscositymodelforfourliquids...............27 5.1Geometryusedinmode1and2simulations...................29 5.2Detailedviewofthegeometrysurroundingthespheresformodes1and2..29 5.3Overviewofthemeshusedformodes1and2..................30 5.4Detailedviewofthemeshusedformodes1and2.Ahighmeshdensityregion isdenedclosetotheinterstitialregionandthemeshgrowscoarsertoward thebulkuidregion................................31 5.5Comparisonof = p and =const.numericalresultsand =const. analyticsolutionformode1.Fluidismineraloil,particlevelocityis=1 rad s .33 5.6Comparisonof = p and =const.numericalresultsand =const. analyticsolutionformode2.Fluidismineraloil,particlevelocityis U =1 m s .35 ix

PAGE 10

5.7Velocityvectorsonthecrosssectionofbinaryspheresformodes1and2, fromlefttoright..................................36 5.8Mode1leftand2rightdiagramofpressuredistribution.'LP'denotes lowpressure,while'HP'denoteshighpressure.................37 5.9Pressuredistributiononthesurfacesofmode1spheresneartheinterstitial region........................................38 5.10Pressuredistributiononthesurfacesofmode2spheresneartheinterstitial region........................................39 5.11Meshschemeformodes3and4,withthreeseparateregionsofdierent elementdensities..................................40 5.12Comparisonof = p and =const.numericalresultsand =const. analyticsolutionformode3.Theuidismineraloilandtherotatingparticle angularvelocityis=1 rad s ............................42 5.13Comparisonof = p and =const.numericalresultsand =const. analyticsolutionformode4.Fluidismineraloil,particlevelocityisU=2.15 m s ..........................................44 5.14Mode4approachvelocityvs.normalizeddragforceforpressure-dependent viscositymineraloilataseparationdistanceof = rad 10 4 .............47 5.15Numericalresultsformode4approachvelocityvs.viscosityandpressurefor pressure-dependentviscositymineraloilataseparationdistanceof = rad 10 4 .48 5.16Individualmotionsofmodes1-3combinedwithmode4motion.......49 5.17Normalizedforceortorqueforcombinationofmodes1and2withmode4. Mode1isnormalizedbythetorqueonanisolatedrotatingparticleandmode 2isnormalizedbytheforceonaparticleinStokesow............52 6.1Meshusedforthemode1solidmechanicsmoduletocomputethedeformationsofthespheres,with143676elementsand600459degreesoffreedom..54 6.2Exaggerateddeformedshapeanddeformationmagnitudesofelasticspheres inmode1motion. = rad 10 E =0 : 2 GPa =0 : 4,scaledto10%ofthetotal geometry......................................55 6.3Exaggerateddeformedshapeofelasticspheresinmode2motion. = rad 10 2 E =0 : 2 GPa =0 : 4,scaledto10%ofthetotalgeometry.Deformationsare normalizedbysphereradius............................57 6.4Meshusedforafully-coupledsimulationofmode4withelasticspheres...60 x

PAGE 11

6.5Exaggerateddeformationofmode4elasticparticlesofE=0.2GPa, = rad 10 3 inahighlyviscousuidwith =0 : 5,normalizedbysphereradius......61 6.6Maximumuidpressureintheinterstitialregionformode4motionofelastic spheresinmineraloil,atarangeofparticleseparationdistances.E=20 GPa, =0.4....................................62 6.7Mode4hydrodynamicdragforceswithelasticparticles,wheretheseparation distancedoesnotaccountforsurfacedeformation...............64 6.8Mode4hydrodynamicdragforceswithelasticparticles.Separationdistance hasbeenadjustedtoaccountforthesurfacedeformation...........64 6.9ComparisonofPoisson'sratiosvs.hydrodynamicforceforE=0.2GPa, = rad 10 3 .......................................66 6.10ComparisonofPoisson'sratiosof0.25and0.49forseveralelasticities.Left: deformationunaccountedforinseparation,right:deformationaccountedfor.68 6.11Comparisonof =0.25and =0.49forE=0.2GPa.Left:deformation unaccountedforinseparation,right:deformationaccountedfor........69 xi

PAGE 12

TABLES TABLE 3.1Meshconvergencestudyresultsforthedragforceonanisolatedspherein Stoke'sow.....................................18 3.2Relativeerrorofnumericalsimulationsforasphereapproachingaplanewall, comparedtotheanalyticsolution........................20 4.1Viscouspropertiesoffourdierentliquids....................26 5.1Detailsofthenumericalmeshschemeusedformode1simulationsateach separationdistance................................31 5.2Analyticandnumericalsolutionsformode1rotatingspherein = const: mineraloil,normalizedbythetorqueonanisolatedsphere..........32 5.3Numericalsolutionsfor = const: and = p numericalsolutionformode 1rotatingsphereinmineraloil,normalizedbythetorqueonanisolatedsphere.32 5.4Comparisonofmaximumgaugepressureandviscosityinconstantandpressuredependentviscositymineraloil,formode1particlemotionwith rot = 1 rad s and = radius 10 4 .................................33 5.5Comparisonofanalyticandnumericalsolutionsformode2movingspherein = const: mineraloil,normalizedbytheforceonasinglesphereinStoke's ow.........................................34 5.6Comparisonof = const: and = p numericalforceonmode2moving sphere,normalizedbytheforceonasinglesphereinStokesow.......35 5.7Comparisonofmaximumgaugepressureandviscosityinconstantandpressuredependentviscositymineraloil,formode2motionwith U =1 m s and = radius 10 4 .36 5.8Detailsofthenumericalmeshschemeusedformode3simulationsateach separationdistance................................41 5.9Comparisonofanalyticalandnumericalsolutionsformode3rotatingsphere inmineraloilwith = const: ,normalizedbythetorqueonanisolatedsphere.41 5.10Comparisonof = const: and = p numericaltorqueonmode3rotating sphereinmineraloil,normalizedbythetorqueonanisolatedrotatingsphere.42 5.11Comparisonofmaximumgaugepressureinconstantandpressure-dependent viscositymineraloil,formode3particlemotionwith=1 rad s and = radius 10 4 .43 xii

PAGE 13

5.12Detailsofthenumericalmeshusedformode4simulationsateachseparation distance.......................................44 5.13Mode4force,max.pressure,andmax.viscosityfor = const: and = p mineraloilatmultipleapproachvelocitiesandaseparationof = rad 10 4 ....45 5.14Increasesinmode4pressure,viscosity,andforcefrom =const.to = p mineraloil,at = radius 10 4 andarangeofapproachvelocities..........46 5.15Torqueonmode1rotatingsphereat 1 =1 rad s whencombinedwithmode4 motionatarangeofapproachvelocities.....................50 5.16Forceonmode2movingsphereat U 2 =1 m s whencombinedwithmode4 motionatarangeofapproachvelocities....................51 5.17Torqueonmode3rotatingsphereat 3 =1 rad s whencombinedwithmode4 motionatarangeofapproachvelocities....................51 6.1Maximumparticledeformationsformode1,correspondingtouidpressures obtainedfromrigidspheresimulation,normalizedbysphereradius......56 6.2Maximumparticledeformationsformode2atseveralseparationdistances andparticleelasticities,correspondingtouidpressuresobtainedfromrigid spheresimulation,normalizedbysphereradius.................58 6.3Maximumparticledeformationsformode3atseveralseparationdistances andparticleelasticities,correspondingtouidpressuresobtainedfromrigid spheresimulation,normalizedbysphereradius.................59 6.4Maximumnormalizedmode4particledeformationsatseveralseparationdistancesandparticleelasticitiesinauidwith =0 : 5 Pa s .Thesimulation isfully-coupledbetweensolidmechanicsandcreepingow...........61 6.5Normalizedforcesonmode4elasticbinaryparticlesatarangeofseparation distancesandparticleelasticities.........................63 6.6Percentdecreaseinmode4forceonelasticparticles,comparedtorigidparticles.63 6.7Normalizedhydrodynamicforcesat = rad 10 3 andarangeofelasticitiesand Poisson'sratio's..................................65 6.8Percentdecreaseinmode4forceonelasticparticles,comparedtorigidparticles.67 6.9Particledeformationsat = rad 1000 forarangeofelasticitiesandPoisson'sratio's.67 xiii

PAGE 14

A.1Convergencestudyformode1three-dimensionalmeshat = rad 10 4 basedon changeinforcebetweensubsequentmeshes...................74 A.2Convergencestudyformode2three-dimensionalmeshat = rad 10 4 basedon changeinforcebetweensubsequentmeshes...................75 A.3Convergencestudyformode3two-dimensionalmeshat = rad 10 4 basedon changeinforcebetweensubsequentmeshes...................75 A.4Convergencestudyformode4two-dimensionalmeshat = rad 10 4 basedon changeinforcebetweensubsequentmeshes...................76 A.5Convergencestudyformode1three-dimensionalsolidmechanicsmeshat = rad 10 4 andE=0.2GPabasedonchangeindeformationbetweensubsequent meshes.......................................77 A.6Convergencestudyformode3two-dimensionalsolidmechanicsmeshat = rad 10 4 andE=0.2GPabasedonchangeindeformationbetweensubsequent meshes.......................................77 A.7Convergencestudyformode4two-dimensionalfully-coupledsimulationmesh at = rad 10 4 andE=20GPabasedonchangeinforcebetweensubsequentmeshes.78 xiv

PAGE 15

Suspensionsofparticlesinviscousruidsareimportantinmanyfeldsandapplicationssuchassedimenttransport,chromatography,oilrecovery,contaminants,andlubricationtheory.Theinteractionsbetweenparticlesinisothermal,lowReynoldsnumberrows,orStoke'srows,havebeenstudiedforalongtimebecauseoftheirapplicabilitytomanyfelds.Particlesareoftenapproximatedasspheresforsimplifcationpurposes.Whentwospheresareclosetooneanotherinshearrows,theirrelativemotioncanbedescribedasacombinationoffourdierentmodes.Thespherescantranslatealongtheirlineofcenters,translateperpendiculartotheirlineofcenters,rotateabouttheirlineofcenters,androtateabouttheaxisperpendiculartotheirlineofcenters.AnalyticsolutionstothesemodeshavebeenderivedandareexactsolutionstotheStokesequations.Theanalyticsolutionsassumeisothermal,lowReynoldsnumber,andincompressiblerows.Theanalyticsolutionsarelimitedtorigidparticles,butiftherearesucientruidpressuresandtheparticlesareelastic,theywilldeform.Theviscosityoftheruidintheseisothermalsolutionsisassumedtobeconstantthroughouttherowfeld,butruidviscosityisactuallydependentonthelocalpressure.Thisresearchwillinvestigatehowparticleelasticityandpressure-dependentruidviscosityaectthehydrodynamicforcesandtorquesonparticles,assumingtherowsareisothermalandincompressible.Thefniteelementmethodpresentsanattractiveanalysistoolforproblemssuchasthisbecauseofitsabilitytoaccuratelycapturethepressuregradientsintheruid.Forexample,whentwospheresareapproachingeachother,theruidpressureisthehighestintheinterstitialregionbetweenthespheres.Wecandefnemultipleelementsacrosstheinterstitialgap,toaccuratelycapturethelargegradientsinpressureandvelocity.COMSOLMultiphysicssoftwarewasusedtosolvetheproblemsnumericallyandallowedforthefne-tuningofmeshdiscretization,boundaryconditions,andsolvingtechniques.1 CHAPTERIINTRODUCTION

PAGE 16

1.1GoverningEquations Thecornerstoneofcomputationaluiddynamicsistheequationsthatgovernthe uidbody.Theseequationsdescribethepressureandmotionofabodyofuidateach pointinthecomputationaldomain.TheuidregionisgovernedbyStoke'sequations, sinceweassumetheviscousforcesdominatetheinertialforces.Thistypeofowisknown asStoke'soworcreepingowandischaracterizedbyverylowReynoldsnumbers.In particular,weassumetheReynoldsnumber,whichisadeterminingcharacteristicofthe ow, Re / inertialforces viscousforces ; .1 issignicantlylessthan1,whichsimpliesthesolutionofthegoverningequations.The Navier-Stokesequationsarecommonlyacceptedasthegoverningequationsforthemotionofisothermaluids.Ingeneral,theNavier-Stokesequationscanbewrittenasthe conservationofmomentumequationandtheconservationofmassequation,whichare solvedtogether.InthecaseofacompressibleNewtonianuidinanisothermalow,the momentumequationiswrittenas @ u @t + u r u | {z } = r p | {z } + r r u + r u T )]TJ/F15 11.9552 Tf 13.151 8.088 Td [(2 3 r u I | {z } + F |{z} ; .2 wherethevector u representsuidvelocityinthethreeorthogonaldirections,and1.2is threeseparateequationswrittenasone,invectornotation.Thecontinuity,orconservationofmass,equationforacompressibleuidiswrittenas @ @t + r u =0 ; .3 where istheuiddensity, p istheuidpressure, t istime, u isthevelocityvector,and 2

PAGE 17

I istheidentitymatrix.Theindividualtermsinequation1.2correspondtotheinertial forces,pressureforces,viscousforces,andbodyforces.Inthecaseofthe motionofbinaryspheresinaviscousuid,thereareafewassumptionsthatcanbemade tosimplifytheseequations. Weassumetheamountoftimehasbeensucientforthesystemtoreachasteady stateandafullydevelopedoweld.Underthissteadystateassumption,allterms dependentontimeareneglectedand @ u @t =0 ; .4 inthemomentumequation.Sincetheviscousforcesdominatetheinertialforces, u r u 0 ; .5 andthistermcanbeneglectedinequation1.2.Thebodyforceisgivenbythegravitationalpotentialwhichiscombinedwiththepressureforces.Ifdensityvariationsin theuidarenegligible,wehaveincompressibleowand canbemovedoutsideofthe gradientinequation1.3,resultingin r u =0 : .6 Substitutingequation1.6intoterminequation1.2,theincompressiblemomentum equationiswrittenas 0= r p + r r u + r u T ; .7 whichissolvedtogetherwithequation1.6. 3

PAGE 18

1.2BinaryParticleInteractions Insuspensionsofsolidparticles,theparticleshaveasignicantaectononeanother. Thisstudywillfocusonhowbinaryparticlesinteractwithoneanotherinaviscousshear ow.Twosphericalparticlescaninteractinacombinationoffourdierentmotionswhich canbedecoupledandsolvedforindependently[18]. Figure1.1:Fourmodesofmotionpossibleforbinaryparticleinteractions. Themotionsingure1.1willbeidentiedwiththeirrespectivelabel.Mode1will indicateamotionwheretwoparticlesarerotatingrelativetooneanotheraboutanaxis perpendiculartotheirlineofcenters.Mode2willdenoteparticlesaretranslatingrelative tooneanotheralongtheaxisperpendiculartotheirlineofcenters.Mode3willbethe motionwheretwoparticlesarerotatingrelativetoeachotherabouttheirlineofcenters. Mode4willdescribethemotionoftwoparticlestranslatingrelativetooneanotheralong theirlineofcenters[18].Thisconventionwillbeusedthroughoutthethesis. 1.3LiteratureReview Theresistanceofasphericalparticlemovingslowlythroughaviscousuidwas studiedbyGeorgeStokesinthemid1800's.Stoke'sLawwasderivedin1851tosolvefor thedragforceontheseparticles,atconstantvelocities[8].Thiswidelyacceptedsolution isexactonlyforuidregionsthatextendinnitelyinalldirections.Thepresenceof boundariesincloseproximitytoparticlessignicantlyperturbstheparticleresistance 4

PAGE 19

formulationandStoke'sLawnolongerapplies[7].Theinuenceofparticlesoneach otherinsuspensionsbecameaninteresttoresearchersbecauseofthepresenceofthese typesofinteractionsinmanynaturalandmechanicalsystems. Intheearly1900's,Einsteinfoundtheviscosityofslurriestoincreasewiththevolume fractionofsolidparticles,anddevelopedanequationtodescribetherelativeviscosity ofdiluteslurriesuptoavolumefractionofapproximately0.02[12].Theequationis limitedtocreepingowsandincompressibleNewtonianuids.GuthandSimhaextended Einstein'sequationin1936forparticlevolumefractionsuptoapproximately0.06,by consideringthebulkinteractionbetweenthesolidparticles[14].Thomas[26]developed animprovedequationin1965usingempiricaldatatoextendEinstein'sworktovolume fractionsofapproximately0.1. Theresearchofdilutesuspensionsnaturallyledtotheresearchofmodelstodescribe higherconcentrationsystems.Forsuspensionowswithhighparticleconcentrations, KriegerandDoughertydevelopedamodelin1959relatingtheviscositytotheparticle packingparameter[20].In1971,Chongetal.[10]investigateddierentratiosofparticle sizesandhowtheirvariationaectsaslurry'sbulkviscosity,ndingthattheminimum viscosityisachievedwithapproximately25%to35%ofthesolidconcentrationbeing nespheres,andtheremainderbeingcoarsespheres.In1972,BatchelorandGreen[5] developedamodelforsemi-dilutesuspensions,whichpredictsaviscosityproportional tothesquareofthevolumefraction.In2012,SohnandKooderivedanexpressionto describethepermeabilityofviscousowthroughapackedbedofbidispersesizedspheres, byconsideringthedragforceoneachsphere[24]. Theresearchofthebulkviscosityofslurriesleadtoariseinpopularityofresearching unaryandbinaryparticlesystems.In1961,HowardBrennerintroducedanexpression forasphereslowlyapproachinganinnite,rigidplanewallinaviscousuid,asamodicationtoStoke'sLaw[7].Inthesametimeframe,Brennerformulatedanexpressionfor asphericalparticleapproachingafreesurface,wherethenormalcomponentofvelocity 5

PAGE 20

andtangentialstressesattheplanesurfacevanish,creatingasymmetrycondition[7]. In1970,MajumdarandO'Neill[21]presentedexactsolutionstobinaryspheresrotating aboutanaxisperpendiculartotheirlineofcentersandtranslatingalongthesameaxis. Thesesolutionsarecorrecttoorder ln ,where istheseparationdistancebetween spheres.In1984,Jerey&Onishi[19]extendedMajumdar&O'Neill'swork,correctto order ln ,andalsopresentedthesolutionforbinaryspheresrotatingabouttheirline ofcenters.In1989,Ascoli,Dandy,andLealdevelopedaboundaryintegralformulation tothehydrodynamicinteractionofasolidparticlewithaplanarwall,obtainingexcellent agreementwithBrenner'ssolution[1]. In2000,Fernandez,Carrica,andDrew[13]calculatedtheaveragedragforceexerted onaparticle,accountingforbinaryparticleinteractions.Theydeterminedthatthedrag forcedependsonparticleconcentration,gradientsofconcentration,anddierencesin spheresizes.In2006,Challa&VanSwol[9]reportedonthedragforceexperienced byasphereapproachingaplanewalltotestthepredictionsofhydrodynamictheory andintroduceVanderWaalsforcestoinvestigatetheireect.Theyfoundthatat lowspeedsandsmallseparationdistances,thepresenceofsolvationforcescanlead toadragforcethatoscillatesbetweenpositiveandnegative.Thisndingcontradicts theclassicalhydrodynamicpredictionofadivergentdragforceasseparationdistances decrease.Theresultcomplicatesthetestingofhydrodynamictheoryforbinaryparticles andsuspensionsofparticles,wherestaticforcesmaybepresent. Inadditiontotheforcesandcouplesoneachindividualsphere,investigationintothe diusionandmigrationofspheresinsuspensionsystemshasbeenperformed.Thisisan importanttopicformanymechanicalsystemsusedinmineralprocessingandsedimentation.In1970,Halow&Willsreleasedtheirexperimentalobservationsofspheremigration inCouettesystemsandprovideamethodtocalculatespheretrajectorieswhenthesphere andliquiddensitiesareequal[16].In1998,Tetlow,Graham,Ingber,andSubiapublished theirndingsofparticlemigrationexperimentsandsimulations[25].Theyobservedan 6

PAGE 21

insignicantaectofparticleroughnessonthemigrationofparticles.In2008,Ingber, Feng,Graham,andBrennerpresentedtheirworkonanewtraction-correctedboundary elementmethodtosolveforthemigrationofroughspheresinnonlinearshearshows[18]. Thisnewmethodwasfoundtobeveryaccurateforpredictingthetrajectoriesofparticles.D'Avino,Snijkers,Pasquino,etal.investigatedthemigrationofasinglespherein aviscoelasticliquidin2012,usingCouetteow[11].Theyfoundthespheretomigrate inthedirectionofdecreasingshearrategradient. Reynoldsestablishedthebasisforalllubricationtheoryproblemsin1886,withthe Reynoldsequation[22].Theequationrelatesthelubricantpressuretothegeometry andmovementofthemachinecomponents.In1991,Hamrock[17]presentedresearch onuidlmlubricationwheretwobearingsurfacesarecompletelyseparatedbyathin layerofuidlm.Rheologyisthestudyoftheowanddeformationofmatter,underappliedforces.Bair[2]providedresearchintolinkingtheeldsofrheologyand elastohydrodynamicEHDlubrication.EHDlubricationproblemsarethosewherethe pressureinthelubricantishighenoughtoresultinanelasticdeformationoftheadjacentbody.Habchi[15]developedaniteelementapproachtoultra-low-viscosityuid EHDfulllmlubricationproblemsinhis2008thesis,usingafully-coupledNewtonRaphsonmethod.ThelubricantpressuresinEHDlubricationproblemsareextremely high,signicantlyaectingtheviscosity.Habchistudiedseveralpressure-viscosityand pressure-temperature-viscositymodelstoemployinhissimulations,givinggreatinsight intothemodelsthataresuitedfordierentcases.Hefoundthattheshapeoftheuid lmregime,determinedbythedeformedshapeofthebody,signicantlyaectsthepressureeldandtheuid'sviscosity.Habchialsodeterminedthatthevariationofpressure acrossthelmthicknessisnegligiblecomparedtothevariationalongthecontactplane. Researchintheareaofuidconstitutiveequations,suchasthosegoverningthe viscosityofauidasafunctionofpressure,ledtomanynewformulations.Sincethe pressureoflubricantsinmachinecomponentsoftenreachvaluesinthegigapascalrange, 7

PAGE 22

itisimportanttohavemodelstodescribethepropertiesofthelubricantsattheseextremepressures.Awidelyusedequationtodescribetherelationshipbetweenviscosity andpressureinauidistheBarusequation,developedbyCarlBarusin1893[3].Hamrock[17]pointsoutthatmanyviscosity-pressurerelationships,suchastheRoelands formulathatwasdevelopedbyRoelandsin1966[23],weredevelopedaftertheBarus equation,becauseitsignicantlyoverpredictsexperimentalviscositiesathigherpressures.Roelandsalsodevelopedamodeltodescribetheviscosityofauidbasedonboth pressureantemperature.BothofRoelands'modelsareempiricallyformulated. In1983,Happel&Brenner[8]bridgedthegapbetweentherstprinciplesofclassical hydrodynamicsanduid-particledynamicsandclariedhowprevioustheoreticalformulationsrelatetopracticalaspectsseeninuidization,sedimentation,andporousmedia. ThistextincludestheformulationofBrenner'ssolutiontoarigidsphereapproachingan inniteplanewall,aswellasinsightintoformulationsbeginningwiththeNavier-Stoke's equations,forotherparticlemotions.Becauseoftheslowmovingnatureoftheproblem, existinganalyticsolutionsarelimitedtoconstantuidviscosityandrigidparticles.No solutioncurrentlyexistsforuidswithpressure-dependentviscosityorelasticparticles. Dependingonthemotionofthebinaryparticlesandthetypeofuid,highpressures mayariseintheoweldandaecttheviscosityoftheuidandthedeformationofthe particles. 8

PAGE 23

2.1RelatedBenchmarkSolutionsIn1851,GeorgeStokesdevelopedStoke'sLaw,whichisanexactsolutiontoaspher-icalparticlemovingthroughaninfnitevolumeofviscousruidataconstantvelocity.Thelawisthebasisforthefallingballviscometer[6]andcanbeexpressedasFd=6rU;(2.1)whereFdisthedragforceonthesphere,isthedynamicviscosityoftheruid,ristheradiusofthesphere,andUisthevelocityatwhichtheparticleismoving.In1961,HowardBrennerconsideredthecasewhereasphericalparticleiscloseenoughtoawallwherethedragforceisaected.BrennerintroducedanexpressiondescribingtheforceonasphericalparticleapproachinganinfniteplanewallasF=6rU;(2.2)whereisamodifcationfactortoStoke'sLaw[7].Thewallissimulatedasanoslipboundarycondition,oraboundaryatwhichtheruidhaszerovelocityrelativetothe9 CHAPTERIIANALYTICSOLUTIONS Manywellknownanalyticsolutionshavebeendevelopedtodescribetheforcesandtorquesactingonparticlesinviscousruids.Abriefoverviewoftheseformulationswillgiveinsightintotheformulationofbinaryparticleinteractions.Allanalyticsolutionsmentionedarelimitedtoisothermal,lowReynoldsnumber,incompressiblerows,andas-sumeaninfniteruiddomain.Theseanalyticsolutionsareimportantfortheverifcationofnumericalsimulations,wherenumericalerrorcanbeintroduced.Wecanidentifyifweareusingthecorrectboundaryconditions,anadequatemesh,andacorrectdomainsizebyusingtheseanalyticsolutions.Oncethenumericalsolutionsagreewiththeanalyticsolutions,otherparameterssuchasruidviscosityandparticleelasticitycanbemodifed.

PAGE 24

boundary.Brenner'ssolutionisintheformofaninniteseries,where isexpressedas = 4 3 sinh 1 X n =1 n n +1 n )]TJ/F15 11.9552 Tf 11.955 0 Td [(1 n +3 2sinh n +1 + n +1sinh2 4sinh 2 n + 1 2 )]TJ/F15 11.9552 Tf 11.956 0 Td [( n +1 2 sinh 2 )]TJ/F15 11.9552 Tf 11.955 0 Td [(1 ; .3 and as =ln h r + s h r 2 )]TJ/F15 11.9552 Tf 11.955 0 Td [(1 ; .4 with h beingthedistancefromthecenterofspheretothewall.Theanalyticsolution forthetorqueonasphererotatinginaninnitedomainofviscousuidisexpressedas T = )]TJ/F15 11.9552 Tf 9.299 0 Td [(8 r 3 ; .5 whereistherotationalvelocityin rad s [6].Equations2.1and2.5willfrequentlybe usedinlatersectionsfornormalizingforcesandtorques. 2.2SolutionstoFourModesofMotion Figure1.1showsthefourpossiblemodesofrelativemotionbetweentwoparticles, withlabelsthatwillbeusedthroughoutthethesis.Theequationsformode1,describing thetorquesactingonbinaryspheresrotatingaboutanaxisperpendiculartotheirline ofcenters,arepresentedbyJerey&Onishi[19].Thesolutionisforthecasewhere onesphereisrotatingandtheotherisstationary.Asamodicationtoequation2.5,the torqueontherotatingspherewithanangularvelocityofis T rotating = )]TJ/F15 11.9552 Tf 9.298 0 Td [(8 r 3 r g r .6 andforthestationarysphere, T stationary = )]TJ/F15 11.9552 Tf 9.299 0 Td [(8 r 3 s g s : .7 10

PAGE 25

Themodicationfactors, g r and g s are g r = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2ln 5 )]TJ/F19 11.9552 Tf 11.955 0 Td [( + B 11 )]TJ/F15 11.9552 Tf 13.151 8.087 Td [(66 )]TJ/F15 11.9552 Tf 11.955 0 Td [(12 +16 2 125 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 2 ln +0 .8 and g s = 2 ln 10 )]TJ/F19 11.9552 Tf 11.955 0 Td [( + B 12 + 2 +24 +43 2 250 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 2 ln +0 ; .9 respectively,where istheminimumseparationdistancebetweenspheres.Thevalues for B 11 and B 12 wereobtainedin[19].The terminequations2.8and2.9relatethe radiusofthetwospheres,whereinthiscasetheseareequal, = )]TJ/F19 11.9552 Tf 14.779 8.088 Td [(r rotating r stationary = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 : .10 Forthiscasewheretheradiusofthespheresareequal, B 11 = : 7029and B 12 = : 0274. Figure2.1:Torqueonbinaryspheres,onerotatingaboutanaxisperpendiculartotheir lineofcentersandonestationary,normalizedbythetorqueonanisolatedrotatingsphere. Thehydrodynamictorqueontherotatingandstationaryspheresincreaseswith decreasedseparationdistance,asshowningure2.1.Thesingularityforthismodeof 11

PAGE 26

motionis log ,wherethesingularitydescribeshowthefunctionapproachesasingular solution.Thetorquevaluesarenormalizedbythetorqueonasingleisolatedsphere, rotatinginthesameuidatthesameangularvelocityastherotatingsphereinthe binarycase. Theequationstosolvefortheforcesonbinaryspherestranslatingalonganaxis perpendiculartotheirlineofcenters,ormode2,areoutlinedbyJerey&Onishi[19]. Thesolutionisforthecasewhereonespheretranslatesatavelocityof U whiletheother isstationary.Theequationsforeachsphereare F translating = )]TJ/F15 11.9552 Tf 9.298 0 Td [(6 r t Uf t .11 and F stationary = )]TJ/F15 11.9552 Tf 9.298 0 Td [(6 r s Uf s ; .12 respectively.ThemodicationfactorsforStokeslaware f t = )]TJ/F15 11.9552 Tf 9.298 0 Td [(4 )]TJ/F19 11.9552 Tf 11.955 0 Td [( +2 2 ln 15 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 3 + A 21 )]TJ/F15 11.9552 Tf 11.444 8.088 Td [(4+45 +58 2 +45 3 +16 4 ln 375 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 4 +0 .13 and f s = )]TJ/F15 11.9552 Tf 9.298 0 Td [(4 )]TJ/F19 11.9552 Tf 11.955 0 Td [( +2 2 ln 15 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 3 + A 22 )]TJ/F15 11.9552 Tf 13.15 8.088 Td [(4 +45 +58 2 +45 3 +16 4 ln 375 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 4 +0 : .14 Thevaluesfor A 21 and A 22 wereobtainedbyJerey&Onishi,andforthecaseofequal radii,are A 21 =0 : 9983and A 22 = )]TJ/F15 11.9552 Tf 9.299 0 Td [(0 : 2737. 12

PAGE 27

Figure2.2:Forceonbinaryspheres,onetranslatingalongaxisperpendiculartolineof centersandonestationary,normalizedby F Stokes Thehydrodynamicforceonbinaryparticlestranslatingrelativetooneanotheralong theaxisperpendiculartotheirlineofcentersisshowningure2.2.Thesingularityfor thismotionis log ,similartomode1.TheforcesarenormalizedwithStoke'sdrag, fromequation2.1,withthesameuidviscosityandatthesameparticlevelocityasthe binaryspheres. In1984,Jerey&Onishipublishedasolutionforbinarysphereswitharelative rotatingmotionabouttheirlineofcenters[19],ormode3.Theformulationsforboth therotatingandstationaryspheresthatareintroducedasamodicationtoequation2.5 canbeexpressedas T rotating = )]TJ/F15 11.9552 Tf 9.299 0 Td [(8 r 3 r h r .15 and T stationary = )]TJ/F15 11.9552 Tf 9.298 0 Td [(8 r 3 s h s ; .16 respectively.Fortherotatingsphere,themodicationfactoris 13

PAGE 28

h r = ; )]TJ/F19 11.9552 Tf 11.955 0 Td [( )]TJ/F15 11.9552 Tf 7.085 -4.339 Td [(1 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 3 )]TJ/F19 11.9552 Tf 26.547 8.087 Td [(ln 2 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 2 +0 ; .17 andforthestationarysphere,themodicationfactoris h s = )]TJ/F19 11.9552 Tf 9.298 0 Td [( 3 ; 1 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 3 )]TJ/F19 11.9552 Tf 20.767 8.087 Td [( 3 ln 2 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 2 +0 : .18 Thefunctionfor iswrittenas z;a = 1 X k =0 k + a )]TJ/F20 7.9701 Tf 6.587 0 Td [(z : .19 Figure2.3:Torqueonbinaryspheres,onerotatingaboutlineofcentersandonestationary,normalizedbythetorqueonanisolatedrotatingsphere. Thenormalizedhydrodynamictorqueactingonbinarysphereswhereonerotates abouttheirlineofcentersandtheotherisstationaryisshowningure2.3.Thereisno singularityforthismodeofmotionandthetorquereachesaconstantvalueasthegap betweenthespheresvanishes. In1961,Brennerdevelopedasolutionforasphereapproachingafreesurface[7].A freesurfaceisusedtoexpressasymmetryboundaryconditionwhichreducesnumerical 14

PAGE 29

costsbyreducingthedegreesoffreedom.Afreesurfaceisexpressedasthevanishingof thenormalcomponentofvelocityandtangentialshearstresses.Thisisthesolutionto binaryspherestravelingalongtheirlineofcenters,whichwehavenamedmode4.The formulationisexpressedas F = )]TJ/F15 11.9552 Tf 9.298 0 Td [(6 rU; .20 where isamodicationfactortoStoke'sLaw,expressedas = 4 3 sinh 1 X n =1 n n +1 n )]TJ/F15 11.9552 Tf 11.955 0 Td [(1 n +3 4cosh 2 n + 1 2 + n +1 2 sinh 2 2sinh n +1 )]TJ/F15 11.9552 Tf 11.955 0 Td [( n +1sinh2 )]TJ/F15 11.9552 Tf 11.955 0 Td [(1 : .21 Figure2.4:Forceonbinaryspherestranslatingalongtheirlineofcenters,normalizedby F Stokes Thenormalizedforceofbinaryspherestravelingalongtheirlineofcentersisshownin gure2.4.Thesolutiontothismotionisdivergent,andthesingularitycanbedescribed as 1 15

PAGE 30

Numericalsolvershavegrowntobeanextremelypowerfulmethodtosolvemanytypesofmathematicalequations.Physicsproblemsareoftengovernedbynonlinearequations,whichareoftenimpossibletosolveanalytically.Numericaltechniquescanoftenbeusedasasubstitute.COMSOLMultiphysicsisacommercialsoftwarepackagewhichusesnumericalmethodstomodelandsimulatemanytypesofphysicsproblems.Thefniteelementmethodwasusedforthesimulationsinthisreport,wherethesoft-waresolvestheequationswithdierentnumericalmethods,dependingonwhetherthereisacombinationofmultiplephysicsorasinglephysicsphenomenon.Thefully-coupledapproachcombinesthedierentphysicsintoonematrixthatissolved,comparedtothesegregatedsolverthatsolvesonematrixforonephenomenonandappliesthesolutionstoaseparatematrixfortheotherphenomenon.Thesoftwarealsoallowsfortheselectionofaniterativesolverandadirectsolver.Directsolversrelyonalargeamountofcom-putermemoryandtheproblemmustbesmallenoughtoinvertthematrixallatonce,whileiterativesolversrequiremultipleiterationstobreaktheproblemupintosmallercomputations.Forthesimulationsinthisthesis,thenonlinearnatureoftheproblemsandthelargeamountofdegreesoffreedomrequiredtheuseofiterativesolvers.COMSOLalsohasa\parametricsweep"featurethatallowsfortheconsecutivesolvingofproblemswithdierentspecifedparameters.Thisfunctionwasusedtocon-secutivelysolveforarangeofparticleseparationdistancesandparticleelasticities.Whentheseparationdistanceisverysmall,thepressurescanbeextremelyhighintheinter-stitialregionanditcanbediculttosolvethegoverningequations.Toremedythis,COMSOLpresentsanoptiontouseaprevioussolutionasinitialconditionsforthecur-rentsimulation.Bydefault,COMSOLassumesinitialconditionsforruidvelocityandpressuretobezero,whichisfarfromthesolution.Byusingaprevioussolutionasthe16 CHAPTERIIINUMERICAL SIMULATIONS

PAGE 31

initialconditioninthesolver,wecanfacilitateconvergencebystartingwithaninitial solutionclosertothenalsolution. COMSOLallowstheusertospecifymaterialpropertiesforbothuidsandsolids. Thiscanbedonebyusingconstantvaluesorbyspecifyingafunctiontogovernthe materialpropertybasedononeorseveralinputparameters.Apressure-dependentanalyticfunctionwasusedtodeterminetheuidviscosityinthesimulationsinvestigating pressure-dependentviscosity.Intheparticleelasticitysimulations,theparticleelasticity andPoisson'sratiowerevariedtodeterminetheiraectonthebinaryparticlehydrodynamics.Torstobtainanunderstandingofthesoftwareenvironment,asphericalparticle inStoke'sowwassimulated.Aspherewitharadiusof0.001 m movingthroughauid withaviscosityof1 Pas wasused.ThenumericaldragforceonthespherewascomparedtotheStoke'sequationandameshconvergencestudywasperformed.Theproblem wassetupwitha2-dimensionalaxisymmetricgeometry,wherethesphereistranslating alongtheaxisofrevolvedsymmetry.Thissignicantlyreducesthenumberofdegreesof freedominthesimulation.Thecylindricaluidregioncontainingthespherewasdened as10 r wideand20 r tall.Theuiddomainwasboundedbyanopenboundaryonthe circumferentialwallofthecylinder,aninowboundaryonthetopofthecylinder,and anoutowboundaryonthebottomofthecylinder. 17

PAGE 32

Table3.1:Meshconvergencestudyresultsforthedragforceonanisolatedspherein Stoke'sow. No.ofElements F Drag [ N ] %Change 110 .01295 N/A 169 .01325 2.23 538 .01638 19.15 2094 .01744 6.05 4658 .01853 5.90 8254 .01882 1.56 12902 .01887 0.24 Table3.1showstheresultsofameshconvergencestudyperformedonthesimulation ofStoke'slaw.Acoarsermeshisusedtobeginwithbecausethecomputationtimecan besignicantlyreducedifthenumberofdegreesoffreedomisminimized.Only110nite elementswereinitiallyusedinthemesh.Themeshdensitywasgraduallyincreaseduntil thedierenceinsolutionbetweentwosuccessivemesheswasbelow1%.Atthispoint,the dierenceinresultsissoinsignicantbetweeneachnewsimulation,thatwecansatisfy themeshconvergencecriteria.Thehighestdensitymeshinthetableresultsinadrag forceof0.01887N,dieringby0.1%fromtheanalyticresultof0.01885Nobtained fromequation2.1.Theexactsolutiontotheproblemisknowninthiscase,butwhen newproblemsarebeingsolved,ameshconvergencestudyisoftennecessarytoguarantee convergenceofthesolution. Arigidsphereinaviscousuidapproachinganinniteplanewallwasalsosolved numerically.Theresultsfortheforceonthespherewerecomparedtotheanalytic solutioninequation2.2.Thisproblemhelpedgaininsightintosettingupthenumerical parametersforthesimulationstofollow. 18

PAGE 33

Figure3.1:Comparisonoftheanalyticandnumericalsolutionsforarigidsphereapproachinganinniteplanewall. Theanalyticandnumericalsolutionsfortheforceonasphereapproachinganinnite planewallisshowningure3.1.Theerrorisinsignicantsuchthatthecurvesappear coincident,sothisdatawasalsotabulated. 19

PAGE 34

Table3.2:Relativeerrorofnumericalsimulationsforasphereapproachingaplanewall, comparedtotheanalyticsolution. Separation, %Error 5r 2.95 r 1.77 r/2 1.18 r/10 0.34 r/100 0.036 r/1000 0.005 r/2000 0.002 r/5000 0.002 r/10000 0.005 Table3.2showsthatthenumericalresultsagreewiththeanalyticsolutionforthe problemwithin3%atthelargestseparationdistanceof5 r andwithin.005%atthe smallestseparationdistanceof radius 10 4 .Brenner'sanalyticsolutionassumesthesphereis approachinganinniteplanewithinaninniteuiddomain.Innumericalsimulations usingtheniteelementmethod,theproblemmustbesetupusinganitedomain enclosingalloftheboundaries,bodies,andphysics.Theregionsizemustbedened andcannotbeinnite.Itwasfoundthatthelargertheregionis,themoreaccurate thesolutionbecomes,butatthecostofasignicantincreaseindegreesoffreedomand computationtime.Thisquicklybecomesadiminishingreturnintermsofaccuracy. Therefore,adecisionmustbemadeastohowlargeofaregionisadequateforthe problem. COMSOLallowsfortheuseofdierentmodules"todeterminetheequationstobe usedtosolvetheunderlyingphysics.Forthepressure-dependentviscosityinvestigation withrigidparticles,thecreepingow"modulewasused,whichusesStokesequations 20

PAGE 35

forincompressibleuids.Fortheproblemsinvolvingelasticparticles,theuidstructure interaction",orFSImodulewasused.Thismoduleallowsforthecombinationoftwo dierentmodules;creepingowfortheuidandstructuralmechanicsforthesolids. TheFSImodulecouplesthetwoproblemssothesoliddeformationsaecttheuid pressuresandvelocitiesandvice-versa.Thesoftwarecombinesthemultiplephysicsby usingasegregatedsolver,segregatingeachofthegoverningequationsandsolvingfor themindividually. 3.1MeshDiscretization Theniteelementmethodreliesonhavinganappropriatemeshdened,wherevaluesarecalculatedatindividualpointswithinamesh.Themeshconsistsofelements, whicharemadeupofnodes.Forthecommontetrahedralelement,therearefourvertices, ornodes,wherethesolutioncouldbecomputed.Thesevaluesmightbedeformations, temperatures,velocities,oranynumberofquantities.Otherelementtypesarethebarelement,typicallyusedintwo-dimensionalstructuralproblems,andthethree-dimensional quadrilateralelement,whichhas8nodes. Therewouldbegapsinthesolutionifthevalueswereonlycalculatedatthenodes. Interpolationfunctionsareusedtoapproximatethesolutionbetweenthenodesusinga polynomialfunction.Inallsimulationsinthisthesis,secondorderfunctionswereused forvelocity,whilelinearshapefunctionswereusedforpressure.Onceanelementtype anddiscretizationordershavebeenchosen,theoverallmeshgeometrymustbecreated. 21

PAGE 36

Figure3.2:Exampleofmeshforbinaryspherestravelingalongtheirlineofcenters. Theuidregionsthathavethelargestgradientinvelocityandpressuremusthave ahighermeshdensity,whereastheareaswithsmallgradientsacrosselementsmayhave largerelements.Itispossibletoknowbeforecreatingameshthatacertainregionwill havelargegradientsinpressureandvelocity.Ingure3.2,therearetwospheresthatare closetogetherandmovingalongtheirlineofcenters,soweexpecttoseehighpressures intheuidandthemeshshouldbeneinthisarea.Ingeneral,itisbesttohavemore thantwoelementsofthicknessinthinareassuchasthisinterstitialregion. Figure3.3:Detailedviewofinterstitialregionmeshforbinaryspheres. Highpressureswithintheinterstitialregionuidarelikelyforbinaryspheresmoving alongtheirlineofcenters.Ingure3.3,weseethattheregionofmostimportanceis wherethegradientsinpressureandvelocityarethelargest,betweenthetwospheres. 22

PAGE 37

Higherdensitymesheswithagreaternumberofdegreesoffreedomcomeattheexpense ofcomputationtime. 3.2BoundaryConditions Choosingthecorrectboundaryconditionsisveryimportantwhensettingupanumericalsimulation.Boundaryconditionsarethemathematicalrepresentationofphysics atinterfaceboundariesinasimulation.Ifuidistravelingadjacenttoasolidbody, thereisano-slipconditionatthesurfaceofthebody.Ano-slipconditionspeciesthat theuidvelocityiszerorelativetothesolidboundary.Theno-slipboundarycondition isgivenby u = u w ;v = v w ;w = w w ; .1 where u w ;v w ; and w w arethex,y,andzcomponentsofthewall'svelocity,respectively.Thisboundaryconditionisusedatanyuid-solidinterface,sotheuidhasthe samecomponentsofvelocityasthesolid,wheretheyareindirectcontact.Theno-slip conditionwasusedinallofthenumericalsimulationsinthisthesis.Anotherexampleof aboundaryconditionthatisoftenusedinnumericalsimulationsisthesymmetrycondition.Asymmetryboundaryconditionspeciesthatthewallhaszeroshearstresses, whichisequivalenttohavingasymmetricboundaryandeectivelymirroringanyphysics ononesideoftheboundarytotheother.Thisconditionisgivenby =0atthesymmetricboundary..2 Symmetryboundariesbecomeimportantwhenonewantstoreducethecomputation timeittakestosolveaproblem.Inthecasewheretwospheresareapproachingone another,suchastherightmostcaseingure1.1,asymmetricplanedirectlybetweenthe twosphereswouldcutthedegreesoffreedomtobesolvedinhalf.Onlargerproblems, thedierencebecomessignicant. 23

PAGE 38

Eachmodeofmotionthatwassimulatedrequiredslightlydierentboundaryconditions.Fortherotationaboutanaxisperpendiculartothespheres'lineofcenters,the surfaceofthespheresbothrequiredseparateboundaryconditions.Thespherethatis stationarysimplyrequiresano-slipboundarycondition,asdescribedinequation3.1. Theoutowboundaryconditionsusedinthesimulationsvarydependingontheparticularmotionofspheres.Inthesimulationofbinaryspheresrotatingaboutanaxis perpendiculartotheirlineofcenters,ano-slipboundarywasappliedonthecircumferentialsurfacesofthecylindricalcontainer.Ontheatendsofthecylinder,anopen boundaryconditionwasapplied,whichspeciestheboundaryisincontactwithalarge volumeofuid,andviscousstressesvanish.Thesameboundaryconditionswereapplied forbinaryspherestranslatingrelativelyaboutanaxisperpendiculartotheirlineofcenters.Theseoutowconditionswerechosenbecausewithallopenboundaries,thesteady statesolutionwasthatthefar-elduidvelocitywasnonzero,whichisinaccurate.This signicantlyreducedtheresistanceonthespheres,toavaluebelowtheanalyticsolution. Imposingano-slipconditiononthecircumferentialcontainerwallsyieldedveryaccurate resultsinregardstotheanalyticsolutions. Forbinaryspheresrotatingabouttheirlineofcenters,ano-slipconditionwasapplied onallcontainerwalls.Thiswasdonesothefar-eldvelocitywaszero,ratherthana velocitygeneratedbytherotationofthesphere,similartomode1.Forbinaryspheres approachingoneanotheralongtheirlineofcenters,whenthespheresareclose,the majorityoftheforceisfromthehighpressuresintheinterstitialregion.Theeectof theoutowboundaryconditionwasfoundtobenegligibleandanopenboundarywas arbitrarilychosen. 24

PAGE 39

Thetwomainareasoffocusforthisresearcharepressure-dependentviscosityandparticleelasticityandhoweachaectsthehydrodynamicforcesandtorquesactingonbinaryparticleinteractions.Theconstitutiveequationsforaruidarethosethatrelatethephysicalquantitiesofsuchasviscosity,density,andpressure,tooneanother.Twomodelswerestudiedandcomparedfordeterminingtheviscosityofaruidbasedonitspressure.Onemodelwaschosentouseasourgoverningequationforaruidsviscosity,basedonliterature[17].4.1ViscosityModelsModelsrelatingpressuretoviscositydatebacktothelatenineteenthcentury.Barusdevelopedthefrstwidelyusedmodelforthepressuredependenceofviscosityinruids[17].Barusstatedthattherelationshipbetweenpressureandviscositycanbeexpressedasln o=p;(4.1)whereistheviscosity,oistheabsoluteviscosityatagaugepressureof0Pa,pisthecurrentpressure,andisthepressure-viscositycoecient.Forover70years,thismodelwasacceptedasthestandard,untilitwasprovedtosignifcantlyoverpredicttheviscosityatpressuresaboveapproximately108Pa.In1966,anothermodelwasdevelopedbyRoelandstocorrectthis[17].Roelandsmodeledtheviscosityatenhancedpressuresaslog+1:20=(logo+1:20)1+p 2000z;(4.2)wherezisadimensionlessviscosity-pressureindexrelatedto.TheBarusequationissaidtobeaccurateuptoapproximately0.1GPa,whiletheRoelandsequationhasa25 CHAPTERIVFLUIDCONSTITUTIVEEQUATIONS

PAGE 40

muchmoreaccuratepredictionforviscosityatpressuresuptoapproximately1.5GPa[4]. Tostudythedierencesbetweenthemodels,theviscosityofmineraloilwasexamined foreach. Figure4.1:ComparisonbetweenBarusandRoelandsmodelsforviscosityofmineraloil. Ingure4.1,itisevidentthatthereisnolimittotheviscosityofthemineraloil accordingtotheBarusequation.Weseethataboveapproximately0 : 1GPa,theRoelands curvebeginstodecreaseinslopeandeventuallyleveloatveryhighpressures.According toexperimentaldata,thisisamoreaccuraterepresentationofhowtheviscosityofuids behavesathighpressures[17],andwewillproceedusingthismodel.Fourdierentuids werechosentoplottheRoelandsmodel,foradiversecomparisonofdierentparameters. Table4.1:Viscouspropertiesoffourdierentliquids. Fluid o [Pas] z Water 0.001 0.1 Polyethyleneglycol 0.09 0.38 Glycerin 1.412 0.18 MineralOil 0.0681 0.67 26

PAGE 41

Table4.1showshowthedimensionlesspressure-viscosityindex, z ,isuniquetoeach uid.Thisparameterdetermineshowtheuidviscosityincreaseswithpressure. Figure4.2:Roeland'spressure-viscositymodelforfourliquids. Itisknownthatwaterisconsideredanincompressibleuid,andingure4.2,we canseethattheviscosityisapproximatelyconstantthroughoutthepressurerange.The Roelandsmodelwasdevelopedformineraloils[17],andwecanseethatmineraloil hasthemostnotableincreaseinviscosityatenhancedpressuresoutofalluidsinthe plot.Polyethyleneglycolandglycerinhavearelativelylowincreaseinviscositywith pressure,incomparison.Glycerinisahighlyviscousuid,asweseewhencomparingits viscosityatambientpressuretoanyotheruidintable4.1.Itbecomesevidentjusthow signicantlypressurecanincreaseviscositywhenconsideringthedierenceinviscosity betweenglycerinandmineraloil,ingure4.2.Mineraloilsurpassestheviscosityof glycerinatagaugepressureofapproximately0.2 GPa 27

PAGE 42

Theviscosityofaruidisgenerallydependentonbothpressureandtemperature.Theanalyticsolutionsforthemodes1-4assumetheruidhasaconstantviscosity.However,theviscositycanvarysignifcantlydependingonthepressurefeldinthedomain.Inthischapter,weconsiderthepressuredependenceofviscosityforisothermalsystems.Theruidviscosityateachdegreeoffreedomwascalculatedusingequation4.2.Thehydrodynamicforcesforallfourmodesofmotioninfgure1.1werecalculatedusingnumericaltechniques,forconstantandpressure-dependentviscosityruid.Mineraloilwasusedastheruidbecauseitsviscosityhasthegreatestdependenceonpressure,comparedtotheotherruidsinfgure4.2.Theelementsinthemeshwerequadraticinvelocityandlinearinpressure.ThegoverningequationsusedinthenumericalsimulationsassumelowReynoldsnumberrowconditions(viscouseectsdominatetherowandtheinertialeectscanbeignored).Thebinaryparticlesimulationswereperformedatconstantparticlevelocities,asopposedtoconstantloads.5.1RotationAboutAxisPerpendiculartoLineofCentersTheeectofpressuredependentviscositywherespheresrotaterelativetooneanotheraboutanaxisperpendiculartotheirlineofcenters(mode1)isanalyzed.Numericalsimulationswereperformedtoprovideadirectcomparisonbetweenpressuredependentandconstantviscosityruids.Thesimulationsrerectasystemofbinarysphereswithonerotatingandonestationary,wherethespheresareatfourdierentseparationdistances.Thesimulationsformodes1and2werethree-dimensional,withthespheresoccupyingacylindricalruiddomainofradius15rsphereandheight30rsphere.Verydensemesheswereconstructedintheinterstitialregion,withelementsgrowinginsizetowardthebulkregionoftheruid.28 CHAPTERVRESULTS:PRESSUREDEPENDENTVISCOSITY

PAGE 43

Figure5.1:Geometryusedinmode1and2simulations. Thegeometryusedforthemode1and2simulationsisshowningure5.1.Theuid regionisacylinderbisectedbyasymmetricplane. Figure5.2:Detailedviewofthegeometrysurroundingthespheresformodes1and2. Theinterstitialregiongeometrywasdividedintofourconcentriccylindricalregions centeredalongthespheres'lineofcenters,eachwithdierentmeshdensities.This geometryisshowingure5.2.Thesmallestofthecylindershadaradiusof r sphere 200 followedbyacylinderofradius r sphere 50 r sphere 5 ,and r sphere 2 .Thesmallestcylinderhada verynemeshwherethehighestuidpressuresandvelocitieswerepresent.Themesh 29

PAGE 44

densitygrewasthecylindersgrewbecausethepressureandvelocitygradientsdecreased farfromtheinterstitialregion. Figure5.3:Overviewofthemeshusedformodes1and2. Themeshusedinthemode1and2simulationsisshowningure5.3.Inthe smallestoftheinterstitialcylinders,triangularelementsweredenedwithamaximum sizeof 3 andsweptintheradialdirectionwithadistributionof40elementsper90 degrees,producingtriangularprisms.Thenextsmallestcylinderwasgivenamaximum elementgrowthrateof1.01.Thenextsmallestwasgivenamaximumelementgrowth rateof1.04.Finally,onthelargestoftheinterstitialcylinders,atransitionwasmade fromthetriangularprismelementstotetrahedralelementswithaCOMSOLpredened sizeof`ner'.Theremainingbulkregionwasgivenapredenedelementsizeof`coarse', sinceadensemeshisnotneededoutsideoftheinterstitialregion. 30

PAGE 45

Figure5.4:Detailedviewofthemeshusedformodes1and2.Ahighmeshdensity regionisdenedclosetotheinterstitialregionandthemeshgrowscoarsertowardthe bulkuidregion. Adetailedviewoftheinterstitialregionofthemeshisshowningure5.4.Thenest elementsareinthesmallestinterstitialcylinderandgrowoutwardtowardthebulkof theuid.Thehighestmeshdensityisintheinterstitialregionwherethespheresarethe closestsincethehighestgradientsinpressureandvelocityarepresenthere.Theswept elementsappearasquadrilateralfromtheinteriorofthespherebecausethesidesofthe triangularmeshcreatedonthesymmetryplanearesweptintheradialdirection. Table5.1:Detailsofthenumericalmeshschemeusedformode1simulationsateach separationdistance. No.Elements Deg.ofFreedom ComputationTime[hr] rad= 10 224683 1383563 1.00 rad= 10 2 158772 1004365 0.75 rad= 10 3 194733 1575838 1.25 rad= 10 4 354150 3875544 2.25 Theresultingmeshstatisticsforthemode1simulationsareshownintable5.1.The largerthedegreesoffreedominamesh,thelongerthecomputationtimebecomes.The 31

PAGE 46

smallestseparationdistanceresultsinthegreatestdegreesoffreedomandcomputation time,requiringapproximately2.25hours. Table5.2:Analyticandnumericalsolutionsformode1rotatingspherein = const: mineraloil,normalizedbythetorqueonanisolatedsphere. Separation, T Numeric =T Isolated T Analytic =T Isolated Percenterror r= 10 1.229 1.233 0.35 r= 10 2 1.632 1.635 0.16 r= 10 3 2.085 2.086 0.06 r= 10 4 2.546 2.545 0.02 Thetorqueobtainedfromthenumericalandanalyticsolutionsarecomparedintable 5.2forconstantviscositymineraloil.Thetableshowsthetorqueontherotatingsphere foreachsolutionandtheerrorinthenumericalsolutions.Theerrorisverylowat everyseparationdistanceandthisindicatesthemeshisadequatetointroducepressuredependentviscositytotheuid. Table5.3:Numericalsolutionsfor = const: and = p numericalsolutionformode 1rotatingsphereinmineraloil,normalizedbythetorqueonanisolatedsphere. Separation, T = const: =T Isolated T = p =T Isolated Percentdierence r= 10 1.229 1.229 0.0 r= 10 2 1.633 1.633 0.0 r= 10 3 2.085 2.085 0.0 r= 10 4 2.546 2.546 0.0 Thetorquesonthemode1rotatingsphereinconstantandpressuredependentviscositymineraloilarecomparedintable5.3.Thedatainthetablesuggeststhatmode 1producesinsucientuidpressurestoincreasetheuidviscositysignicantly,even 32

PAGE 47

atverysmallseparationdistances.Thepercentdierencecolumnshowsthatintroducingpressure-dependentviscositydidnotaectthetorqueonthespheres,comparedto constantviscosityuid. Figure5.5:Comparisonof = p and =const.numericalresultsand =const. analyticsolutionformode1.Fluidismineraloil,particlevelocityis=1 rad s Thedataintables5.1and5.2arecombinedingure5.5tocomparetheanalytic solutiontotheconstantandpressuredependentnumericalsolutions.Allthreecurves lieontopofeachother,showingthatthenumericalresultscorrelateverywellwith theanalyticsolutionwithminimalerrorandthattheconstantandpressure-dependent viscosityuidsyieldthesametorques. Table5.4:Comparisonofmaximumgaugepressureandviscosityinconstantandpressure dependentviscositymineraloil,formode1particlemotionwith rot =1 rad s and = radius 10 4 ViscosityModel p gauge;max: [ Pa ] max: [ Pa s ] = const: 26550 0.0681 Roelands; = p 26558 0.0681 Table5.4containsthemaximumuidgaugepressureandviscosityintheuidre33

PAGE 48

gionfortheconstantandpressure-dependentviscositysimulations.Thedataisforthe numericalsimulationsofbinaryspheresinmineraloilwitharelativeangularvelocityof 1 rad s andaseparationdistanceof radius 10 4 .Therearenosignicantchangesinuidpressure orviscosity,whichreinforcesthedataintable5.3. 5.2TranslationAlongAxisPerpendiculartoLineofCenters Binaryspherestranslatingrelativelyalonganaxisperpendiculartotheirlineof centersmode2isconsidered.Onespherewasgivenaconstantvelocitymotionperpendiculartothelineofcenters,whiletheotherspherewasstationary.Weconsidertheforce onthemovingspheretoanalyzethehydrodynamiceects.Athree-dimensionalsimulationwasperformed,usingthesamegeometryandmeshasthesimulationsinsection 5.1,withverysimilarcomputationtimes. Table5.5:Comparisonofanalyticandnumericalsolutionsformode2movingspherein = const: mineraloil,normalizedbytheforceonasinglesphereinStoke'sow. Separation, F Numeric =F Isolated F Analytic =F Isolated Percenterror r= 10 1.389 1.393 0.29 r= 10 2 1.761 1.767 0.36 r= 10 3 2.162 2.150 0.55 r= 10 4 2.515 2.533 0.74 Weseeaverylowrelativeerrorateveryseparationdistancebycomparingtheanalytic andnumericsolutionstotheforceonthemovingsphereintable5.5.Thisdemonstrates thatthemeshbeingusedisadequatetocapturethephysicsofmode2. 34

PAGE 49

Table5.6:Comparisonof = const: and = p numericalforceonmode2moving sphere,normalizedbytheforceonasinglesphereinStokesow. Separation, F = const: =F Isolated F = p =F Isolated Percentdierence r= 10 1.389 1.389 0.0 r= 10 2 1.761 1.761 0.0 r= 10 3 2.162 2.162 0.0 r= 10 4 2.515 2.515 0.0 Table5.6comparesthenumericalresultsforconstantandpressuredependentuid viscosity.Theresultsshowthatthetwocaseshavenearlyidenticalresultsandthe pressureintheuidisnotsucienttoincreasetheviscosityforthepressure-dependent uidsimulation. Figure5.6:Comparisonof = p and =const.numericalresultsand =const. analyticsolutionformode2.Fluidismineraloil,particlevelocityis U =1 m s Figure5.6showsthecorrelationbetweenthemode2constantandpressure-dependent viscositynumericalsolutionsandtheanalyticsolution.Thecurvesareallcoincident, meaningthattheerrorintheconstantviscositynumericalsolutionandthedierence betweentheconstantandpressure-dependentviscositysimulationsarenegligible. 35

PAGE 50

Table5.7:Comparisonofmaximumgaugepressureandviscosityinconstantand pressure-dependentviscositymineraloil,formode2motionwith U =1 m s and = radius 10 4 ViscosityModel p gauge;max: [ Pa ] max: [ Pa s ] = const: 1.089 0.0681 Roelands; = p 3.4136 0.0681 Table5.7showsthemaximumgaugepressureintheuidforboththeconstantand pressure-dependentviscosityuids,obtainedfromthenumericalsimulations.Theresults areforbinaryspheresinmineraloilwitharelativevelocityof1 m s ataparticleseparation distanceof radius 10 4 .Thereisnotasignicantdierenceinuidpressureandviscosity betweenthemode2constantandpressure-dependentviscosityuid.Themaximum pressuresintable5.7arenoticeablylowerthanthoseintable5.4,whenthespheresin bothmodeshavethesametangentialvelocityattheirclosestpointintheinterstitial region. Figure5.7:Velocityvectorsonthecrosssectionofbinaryspheresformodes1and2, fromlefttoright. Thecomponentsofvelocityformodes1and2areshowningure5.7.Theupper sphereinmode1isrotatingaboutitscenterpoint,withoutanytranslation.Thesteadystatesolutiontomode2wassolvedattheinstantintimetheupperspherewastraveling 36

PAGE 51

perpendiculartothespheres'lineofcenters. Figure5.8:Mode1leftand2rightdiagramofpressuredistribution.'LP'denotes lowpressure,while'HP'denoteshighpressure. Thedistributionofpressuresonthesurfacesofthespheresneartheinterstitialregion isshowningure5.8.Therotationinmode1entrainstheuidintothewedge-shape formedbythespheres'surfaces,causinganincreaseinpressure.Ontherightsideofthe sphere,thesurfacesofthespheresareseparation,causingalowpressure.Attheinstant themovingmode2sphereistranslatingperpendiculartothespheres'lineofcenters, itisdraggingtheuidonitsleftside,creatingapressuredropbehindit.Thisuid thencollideswiththelowersphere,creatinganincreaseinpressuresonitssurfaces.The oppositeoccursontherightsideofthelineofcenters,wheretheuppersphereispushing theuidbody,causinganincreaseinuidpressures.Theuidadjacenttothelower sphereisbeingpulledawayfromitssurface,creatinglowpressures. 37

PAGE 52

Figure5.9:Pressuredistributiononthesurfacesofmode1spheresneartheinterstitial region. Thepressuredistributionsalongthesurfacesofthemode1spheresneartheinterstitialregionisshowningure5.9.Thecurvesarealmostcoincident,meaningthepressure intheinterstitialregionisafunctionofthex-coordinateonly.Theuidisentrainedinto thewedge-shapemadebetweenthespheres,bytherotationoftheuppersphere. 38

PAGE 53

Figure5.10:Pressuredistributiononthesurfacesofmode2spheresneartheinterstitial region. Thepressuredistributionsalongthecurvedinterstitialsideofthemode2spheres areshowningure5.10.Thegureshowsthattheuidadjacenttothemovingsphere isincreasinglynegativeleftofthelineofcenters,asweapproachtheinterstitialregion. Atthepointwherethesphere'saretheclosest,thepressureis0,followedbyanupward spikeinpressuretherightofthelineofcenters.Thestationaryspherehasexactlythe oppositetrend. 5.3RotationAboutLineofCenters Pressure-dependentviscositywassimulatedforbinarysphereswitharelativerotation abouttheirlineofcenters.Thismodeofmotionhasnosingularityandtheanalytic solutioningure2.3showsthatthetorqueonthespheresstopsincreasingatsmall separationdistances,belowapproximately r 100 .Sincethespheresarerotatingrelative tooneanotherabouttheirlineofcenters,thepointwheretheirsurfacesaretheclosest hasavelocityof0andnoshearingwilloccurhere.Thesimulationwastwo-dimensional 39

PAGE 54

axisymmetric,withacylindricalcomputationaldomainofradius15 r sphere andheight 30 r sphere Figure5.11:Meshschemeformodes3and4,withthreeseparateregionsofdierent elementdensities. Figure5.11showsthemeshschemeusedformodes3and4.Tocapturethepotential gradientsinuidpressureandvelocityintheinterstitialregion,severalregionsofdierent elementdensitieswereused.Themeshdensityisveryhighintheinterstitialregionand theelementsgrowtowardthebulkregionoftheuid.Thesmallestregioningure5.11 hasamaximumelementsizeof 4 ,thenextsmallesthasamaximumelementgrowthrate of1.025,andthelargestbulkuidregionhasthepredenedCOMSOLmeshdensity of`ner'.Withthismeshscheme,theconvergencewasmetandtheaccuracyofthe numericalconstantviscositysolutionatallseparationdistanceswaswithin0.05%ofthe analyticconstantviscositysolution. 40

PAGE 55

Table5.8:Detailsofthenumericalmeshschemeusedformode3simulationsateach separationdistance. No.Elements Deg.ofFreedom ComputationTime[s] rad= 10 7190 33158 4 rad= 10 2 9271 43005 4 rad= 10 3 11192 52667 4 rad= 10 4 44862 210632 10 Table5.8showsdetailsofthemeshusedingure5.11includingthenumberofelements anddegreesoffreedomforeachseparationdistance.Themeshbecomesdenserinthe areabetweenthespherestoaccountforthelargepressureandvelocitygradientsinthe interstitialregion.Theelementsmustgrowatalowenoughratetoavoidbecoming stretchedtoanirregularlyshapedtriangleandreducingthemeshquality.Thisresults inamuchgreaternumberofelementsanddegreesoffreedomforthesmallerseparation distancesmeshes.Theresultingmeshfor = rad 10 4 has44862elementsand210632degrees offreedom. Todeterminetheaectofpressure-dependentviscosityuidonthehydrodynamic forcesonparticles,mineraloilwithaconstantandpressure-dependentviscositywas simulatednumerically. Table5.9:Comparisonofanalyticalandnumericalsolutionsformode3rotatingsphere inmineraloilwith = const: ,normalizedbythetorqueonanisolatedsphere. Separation, T Numeric =T Isolated T Analytic =T Isolated Percenterror r= 10 2 1.047 1.058 1.01 r= 10 3 1.051 1.053 0.13 r= 10 4 1.052 1.052 0.0 Table5.9showsthemeshisadequatetocapturethephysicsformode3,byagreement 41

PAGE 56

betweentheanalyticsolutionandthenumericalsolution.Twosphereswitharelative rotationabouttheirlineofcentershavelittletonoshearstressesintheinterstitialregion betweenthespheres.Thetorqueonthespheresisalmostconstantatdierentseparation distancessincethemajorityoftheshearstressarisesfarfromtheinterstitialregion. Table5.10:Comparisonof = const: and = p numericaltorqueonmode3rotating sphereinmineraloil,normalizedbythetorqueonanisolatedrotatingsphere. Separation, T = const: =T Isolated T = p =T Isolated Percentdierence r= 10 2 1.047 1.047 0.0 r= 10 3 1.051 1.051 0.0 r= 10 4 1.052 1.052 0.0 Table5.10showsthattheintroducingpressure-dependentviscositytothemineral oildoesnotresultinatorqueincrease.Thisisbecausethevelocityonthesurfaceofthe spheresisnegligibleintheinterstitialregionandthepressureandviscosityoftheuid remainconstant. Figure5.12:Comparisonof = p and =const.numericalresultsand =const. analyticsolutionformode3.Theuidismineraloilandtherotatingparticleangular velocityis=1 rad s 42

PAGE 57

Figure5.12showsthattheanalyticandnumericalsolutionstomode3agreeverywell atcloseseparationdistancesandthatintroducingpressure-dependentviscositydoesnot increasethetorqueontheparticle.Thisisbecausetheuidpressuresincreasenegligibly andthereisnochangeinviscosity.Theslightdeviationintheanalyticandnumerical solutionsatthelargestseparationdistancesisbecausetheanalyticsolutionislimitedto verycloseseparationdistances. Table5.11:Comparisonofmaximumgaugepressureinconstantandpressure-dependent viscositymineraloil,formode3particlemotionwith=1 rad s and = radius 10 4 ViscosityModel p gauge;max: [ Pa ] max: [ Pa s ] = const: 0 0.0681 Roelands; = p 0 0.0681 Table5.11containsthemaximumgaugepressureandviscosityintheuidforthe numericalsimulationofbinaryspheresinconstantandpressure-dependentviscositymineraloilwitharelativeangularvelocityof1 rad s andseparationdistanceof radius 10 4 .The numericalsimulationsdeterminedthepressureintheuidwasidenticalintheconstant andpressure-dependentviscosityuids. 5.4TranslationAlongLineofCenters Binaryspherestranslatingalongtheirlineofcentersmode4resultsinthehighest uidpressuresoutofallfourmodesofmotion.Thesingularityforthismotionis 1 and theincreaseinforceiscausedbythesignicantincreaseinuidpressuresbetweenthe spheres.Thehighpressuresintheinterstitialregionactonthesurfacesofthespheres toincreasethehydrodynamicforcesonthem.Thenumericalsimulationforthismotion wasatwo-dimensionalaxisymmetricformulation,usingageometryandmeshsimilar tothesimulationsasdiscussedinsection5.3.Theonlydierenceisthatmode4can besolvedaxisymmetricallyandwecanusesymmetrybetweenthespheres,onaplane 43

PAGE 58

perpendiculartotheirlineofcenters. Table5.12:Detailsofthenumericalmeshusedformode4simulationsateachseparation distance. No.Elements Deg.ofFreedom ComputationTime[s] rad= 10 6797 31412 4 rad= 10 2 9949 46201 4 rad= 10 3 13214 62226 4 rad= 10 4 65284 304581 12 Themeshstatisticsforeachseparationdistancesimulatedareshownintable5.12. Thecloserthespheresaretooneanother,thedenserthemeshmustbeintheinterstitial regiontocapturethelargegradients. Figure5.13:Comparisonof = p and =const.numericalresultsand =const. analyticsolutionformode4.Fluidismineraloil,particlevelocityisU=2.15 m s Thehydrodynamicforcesincreasesubstantiallyasthespheresbecomecloser,asshown 44

PAGE 59

ingure5.13.Theanalyticsolutionforthismotionisplottedtoshowthatthemeshwas adequatelyconvergedforthesimulations,byitsagreementwiththenumericalconstant viscosityresults.Atthisparticlevelocityof U =2 : 15 m s inmineraloilwithpressuredependentviscosity,thehydrodynamicforceis39%higheratthesmallestseparation distanceof radius 10 4 ,comparedtoconstantviscosityuid. Table5.13:Mode4force,max.pressure,andmax.viscosityfor = const: and = p mineraloilatmultipleapproachvelocitiesandaseparationof = rad 10 4 U[m/s] =const. Roelands, = p p max: [Pa] F=F Stokes p max: [Pa] max: F=F Stokes 0.25 5.108e6 5005.39 5.446e6 0.077 5112.28 0.5 1.021e7 5005.39 1.169e7 0.089 5230.53 0.75 1.532e7 5005.39 1.902e7 0.106 5363.54 1 2.043e7 5005.39 2.785e7 0.130 5515.42 1.5 3.065e7 5005.34 5.367e7 0.233 5906.57 2 4.086e7 5005.24 1.183e8 0.914 6567.58 2.1 4.291e7 5005.42 1.586e8 2.050 6800.15 2.15 4.393e7 5005.33 2.058e8 5.026 6961.23 Table5.13showsthedierenceinforceontheparticlesinaconstantviscosityand pressure-dependentviscositymineraloilusingtheRoelandsmodel,atthesmallestseparationdistanceof radius 10 4 .Thetablecoversarangeofapproachvelocitiesupto2.15 m=s beyondwhichconvergencewasunobtainable.Weseethatthegreatertheparticlevelocity,thegreaterthechangeinuidpressureandviscosity.Thisleadstolargerchangesin forceforthepressure-dependentviscosityuid,comparedtoconstantviscosity. 45

PAGE 60

Table5.14:Increasesinmode4pressure,viscosity,andforcefrom =const.to = p mineraloil,at = radius 10 4 andarangeofapproachvelocities. U[m/s] Percentincreasein maximumpressure Percentincreasein maximumviscosity Percentincrease inforce 0.25 6.6 13.8 2.1 0.5 14.5 31.8 4.5 0.75 24.2 56.3 7.2 1 36.3 91.5 10.2 1.5 75.1 241.7 18.0 2 189.4 1242.4 31.2 2.1 269.7 2910.6 35.9 2.15 367.4 7281.4 39.1 Table5.14containstherelativedierencesinpressure,viscosity,andforceformode 4motionbetweenpressure-dependentandconstantviscositymineraloil.Themaximum viscosityforthepressure-dependentsimulationisapproximately7300%greaterthan thatoftheconstantviscositysimulationatthehighestapproachvelocitysimulated.The resultinghigherpressuresleadtoanapproximate39%greaterforcethantheconstant viscositysimulationsatthisapproachvelocityof2.15 m=s 46

PAGE 61

Figure5.14:Mode4approachvelocityvs.normalizeddragforceforpressure-dependent viscositymineraloilataseparationdistanceof = rad 10 4 Figure5.14showstheincreaseintheforceonmode4particles,normalizedbythedrag forceonasingleparticleinStokesow.Theuidpressureandviscosityaredependent ontheapproachvelocity,whichalsoresultsinchangesinhydrodynamicforceonthe particles. 47

PAGE 62

Figure5.15:Numericalresultsformode4approachvelocityvs.viscosityandpressure forpressure-dependentviscositymineraloilataseparationdistanceof = rad 10 4 Theincreaseinuidpressureandviscositybasedonthemode4approachvelocityfor pressure-dependentmineraloilisshowningure5.15.Theparticlesareataseparation distanceof = radius 10 4 andwecanseethatasthevelocityincreases,thepressureand viscosityincreasenonlinearly. 5.5CombinedMotions Mode4motiongeneratessucientuidpressuresintheinterstitialregionofthe spherestoincreasetheuid'sviscositysignicantly,whilemodes1-3generateinsucient pressurestoincreasetheviscosity.Themotionofparticlesinsuspensionsisrealistically acombinationofthefourmodes,sowewillinvestigatewhetherthecombinationofmode 4withmodes1-3aectsthehydrodynamicforcesonthesemodes.Combiningmode 4withtheothermotionsimposesthehighpressuresandthechangeinuidviscosityof mode4ontheothermodes,todeterminetheiraectoneachmode. 48

PAGE 63

Figure5.16:Individualmotionsofmodes1-3combinedwithmode4motion. Wewillconsidereachcombinedmotioningure5.16.Eachcombinedmotionrequires adierentgeometrysincemode3and4aresolvableinatwo-dimensionalaxisymmetric conguration,whiletheothersrequirethreedimensions.Wecanmakeuseofsymmetry whencombiningmodes1and4onaplanebisectingthespheres,perpendiculartothe axisofrotation.Combiningmodes2and4alsorequiresathree-dimensionalgeometry, butwecanagainmakeuseofsymmetryonaplanebisectingthespheres.Finally,modes 3and4weresolvablewithatwo-dimensionalaxisymmetricgeometry,sinceCOMSOL allowsanangularvelocityabouttheaxisofsymmetry.Sincethehighestincreasesin viscositywereobservedat = rad 10 4 ,thecomparisonswereperformedatthisseparation distance,alsoinmineraloil. Thesamethree-dimensionalmeshdescribedinsection5.1wasusedforcombining modes1and4aswellasmodes2and4.Forcombiningmodes3and4,thetwodimensionalaxisymmetricgeometryandmeshfromsection5.3wasused.Thecomputationtimesforallcombinedmotionswereverysimilartothoseintheseprevioussections, withtheirrespectivemeshes.Modes1and4werecombinedbyspecifyingamoving boundaryconditiononeachspheretowardoneanotheralongthelineofcenters,while alsospecifyingarotationalboundaryconditionononeofthespheresabouttheaxis perpendiculartothelineofcenters.Formode2,thesameapproachvelocitiesweregiven tomode4andamovingboundaryconditionwasgiventomode2.Mode3issimilar 49

PAGE 64

tomode1inthatthespheresareapproachingoneanother,andonesphereisrotating. Theresultsfromeachcombinationofmodesarecomparedtotheresultsfromtheindividualmodes1-3intheabsenceofmode4.Arangeofmode4approachvelocities areconsideredsincevastlydierentpressuresandchangesinviscosityresultforeach.A singlemode1angularvelocityof 1 =1 rad=s wasused,whichcanbecomparedwith theresultsfromsection5.1. Table5.15:Torqueonmode1rotatingsphereat 1 =1 rad s whencombinedwithmode 4motionatarangeofapproachvelocities. Mode4Velocity, U 4 [ m=s ] T Mode 1 T Isolated Percentdierence from U 4 =0 0Mode1Only 2.545 N/A 0.5 2.559 0.56 1 2.572 1.07 1.5 2.597 2.08 2 2.691 5.75 Thetorqueonamode1rotatingsphereat 1 =1 rad s fordierentmode4approach velocitiesisshownintable5.15.Anapproachvelocityof0 m s isincludedformode1 motionintheabsenceofmode4motion.Thepercentdierencecolumnshowshowthe torquechangedfromonlyhavingmode1motion,tothecombinedmotionatthespecied approachvelocity.Wecanseethemode1torqueincreasingasthemode4approach velocityincreases,witha5%increaseintorqueatthelargestapproachvelocityof2 m s 50

PAGE 65

Table5.16:Forceonmode2movingsphereat U 2 =1 m s whencombinedwithmode4 motionatarangeofapproachvelocities Mode4Velocity, U 4 [ m=s ] F Mode 2 F Stokes Percentdierence from U 4 =0 0Mode2Only 2.515 N/A 0.5 2.535 0.83 1 2.570 2.21 1.5 2.621 4.22 2 2.757 9.62 Themode2forcewhencombinedwithdierentmode4approachvelocitiesisshown intable5.16.Asintheprevioustable,anapproachvelocityof0 m s meansonlymode2 motionexists.Thetableshowsthatanincreaseinmode4approachvelocityincreases themode2force.Atthelargestapproachvelocityof2 m s ,themode2forceincreasesby 9.6%comparedtomode2spheresintheabsenceofmode4motion. Table5.17:Torqueonmode3rotatingsphereat 3 =1 rad s whencombinedwithmode 4motionatarangeofapproachvelocities Mode4Velocity, U 4 [ m=s ] T Mode 3 T Isolated Percentdierence from U 4 =0 0Mode3Only 1.052 N/A 0.5 1.052 0.0 1 1.052 0.0 1.5 1.052 0.0 2 1.052 0.0 Table5.17showsthatthechangeinmode3torqueisnegligiblewhenmode4motionis combinedwithmode3motion.Althoughextremepressuresarepresentintheinterstitial 51

PAGE 66

regionbetweenthespheres,thevelocityrelatedtothemode3motionisverylowonthis areaofthesphereandtherearenoshearstressestoinuence.Theuidpressurefar fromtheinterstitialregionisessentiallyunaectedbytheparticlesapproachingandwe cansaymode4motionhasnoeectonthetorqueonmode3particles. Figure5.17:Normalizedforceortorqueforcombinationofmodes1and2withmode 4.Mode1isnormalizedbythetorqueonanisolatedrotatingparticleandmode2is normalizedbytheforceonaparticleinStokesow. Thenormalizedforceortorqueintables5.15and5.16arecombinedingure5.17 bycomparingtheapproachvelocitytothenormalizedforceortorqueontherotatingor translatingsphere.Theplotshowsthatmode2hasagreaterincreaseinforcecompared totheincreaseintorqueformode1.Mode3isnotincludedinthegurebecausethere isnoincreaseintorqueduetotheadditionofthemode4motion. 52

PAGE 67

Theworkpresentedsofarhasconsideredrigidparticles.Theparticlesinsuspen-sionsystemsmayhaveasucientlylowelasticityandtheruidmaybeathighenoughpressurestoresultinparticledeformation.Thissectionwillinvestigatehowparticleelasticityaectstheindividualmodesofmotioninbinaryinteractions.Inparticular,wewillconsiderhowthedeformationaectsthehydrodynamicforcesandtorquesontheparticles.ThesimulationsarebasedonStokesequationsfortheruid,coupledwithlinearelasticityequationsforthecalculationofelasticdeformation.TheanalysiswillconsiderfourdierentYoung'smoduliof200GPa,20GPa,2GPa,and0.2GPa,comparabletothoseofthecommonmaterials,steel,concrete,ABSplastic,andsoftrubbers,respectively.APoisson'sratioof0.4isusedwhencomparingdierentelasticities.ThelastsectioninthischapterwillconsiderarangeofPoisson'sratiovaluesforeachelasticity.Investigatingsuchawiderangeofparticleelasticitieswillhelpgiveinsightintowhetherthereisadierencebetweenlargeandsmalldeformationsintermsofthehydrodynamicdragforceandtorque.Thehydrodynamicforcesandtorqueswerenormalizedwiththehydrodynamicforcesandtorquesforanisolatedsphereinthesameviscousruidastheelasticitysimulations,allowingustodisregardtheunitsandthetypeofruid.6.1RotationAboutAxisPerpendiculartoLineofCentersMode1requiresathree-dimensionalgeometrybutitispossibletodescribethismotionwithasymmetricmodel,whichreducesthedegreesoffreedomsignifcantly.Theresultsintheprevioussectionsregardingpressuredependentviscosityreliedononlythecreepingrowmodule.Thesphereswereassumedtoberigidandtherewasnoneedtosimulatetheirstressesandstrains.Introducingelasticitytotheparticlesrequiresthecombinationofcreepingrowandsolidmechanicsphysics.COMSOLhasamodulenamed"ruid-structureinteraction"whichallowsforthetwo-waycouplingofbothof53 CHAPTERVIRESULTS:PARTICLEELASTICITY

PAGE 68

thesephysics.Applyingthesameboundaryconditionsasmode1fortherigidspheres andassigningaYoung'smodulus, E ,toeachsphere,weattemptedtorunthesimulation. COMSOLwasunabletoconvergeonaresultinthisthree-dimensionalsimulation,and theproblemwasdecoupled. Ifthecouplingbetweenphysicsisremoved,wecanobtaintheuidpressuresadjacent tothespheres'surfacesusingthecreepingowmodule.Next,thesepressurescanbe mappedviaaninterpolationfunctionontothesurfacesoftwoelasticspheresinaseparate solidmechanicssimulation.Thesolidmechanicssimulationdoesnotexplicitlyhavethe rotatingboundaryconditionappliedtothespheres,butratherapressureeldboundary load,asiftheywererotating.Oncethesepressuresaremappedtothesurfacesofthe spheresandthesimulationisperformed,wehaveinsightintowhetherthepressuresare highenoughtosignicantlydeformthesespheresforagivenelasticity.Thegeometry andmeshusedinthecreepingowportionofthesimulationwerethesameasinthe viscositysimulationformode1,insection5.1.Thesolidmechanicsmoduleonlyrequires themeshingofthetwospheressincethepressureeldsareprescribedfromtheuid simulation. Figure6.1:Meshusedforthemode1solidmechanicsmoduletocomputethedeformationsofthespheres,with143676elementsand600459degreesoffreedom. Asshowningure6.1,asmallersphericalregionofradius r sphere 10 wascreatedin 54

PAGE 69

thecenterofeachsphereandaxedconstraintwasdenedhere.Insolidmechanics simulationswhereexternalforcesareapplied,alldeformationswouldbeinniteifsome degreesoffreedomarenotxedinspace.Symmetrywasappliedtotheplanebisecting thespherestoreducethedegreesoffreedom,similartothecreepingowsimulation. Theresultingmeshconsistsof143676elements,600459degreesoffreedom,andrequired approximately9minutestosolve.Thesamemeshwasusedforeveryseparationdistancesincethepressuresintheinterstitialregionarealreadydeterminedinaseparate simulationandarebeingmappedontotheelasticspheres'surfaces. Figure6.2:Exaggerateddeformedshapeanddeformationmagnitudesofelasticspheres inmode1motion. = rad 10 E =0 : 2 GPa =0 : 4,scaledto10%ofthetotalgeometry. Theresultingdeformationfrommode1isshowningure6.2,whereweseeacompressionleftofthelineofcentersandanexpansiontotherightofthelineofcenters. Thedeformedshapeisexaggeratedbecauseactualdeformationisnotvisibletotheeye. Inthesimulation,thelowersphereisstationaryandtheuppersphereisrotatingcounterclockwise.Theuppersphereisdragginguidintothegapandcompressingintothe wedgeshapeformedbythespheres'curvature.Thisresultsinraisedpressuresanda compressionofthesurfacestotheleftofthelineofcenters.Thespheres'surfacesto therightofthelineofcentersareseparatingandcausingavacuumwhichcontractsthe surfacesofthesphereinward. 55

PAGE 70

Thestationarysphereexperiencesalargerdeformationmagnitudebecauseitonly hasalargerightwardforceintheinterstitialregion.Thereisnoforceontheopposite sideofthesphere,counteractingthisrightwardforce.Incontrast,theupperspherehas aleftwardforceintheinterstitialregion,butalsohasarightwardforceonitsopposite side,fromtherotationalresistance.Wecanseetheupperspherestaysrelativelycentered initsoriginalposition,whileallpointsinthelowerspherehaveanetrightwardshift, whichcausesthelargedeformationmagnitude.Thiscanbesummarizedbetheupper spherehavingmorebalancedshearstressesonitssurfacethanthelowersphere. Table6.1:Maximumparticledeformationsformode1,correspondingtouidpressures obtainedfromrigidspheresimulation,normalizedbysphereradius. Separation, Max:Deformation=ParticleRadius E =200GPa E =20GPa E =2GPa E =0 : 2GPa r= 10 2.42e-12 2.42e-11 2.42e-10 2.42e-9 r= 100 1.34e-11 1.34e-10 1.34e-9 1.34e-8 r= 1000 1.03e-10 1.03e-9 1.03e-8 1.03e-7 r= 10000 1.13e-9 1.13e-8 1.13e-7 1.13e-6 Themaximumdeformationoftherotatingsphereisshownintable6.1,wherewesee thatevenatthesmallestseparationdistanceandlowestelasticmodulus,anormalized deformationofonly1.13e-6isexperienced.Thisvalueissoinsignicantthatwedon't expectachangeinhydrodynamicforcescomparedtorigidspheres.Sincethetwodierent physicsaredecoupled,thepressuresobtainedareforrigidsphereswhichweexpectto belargerthanthepressuresonelasticspheres.Iftherearehighuidpressuresinthe interstitialregion,thegapwillgrowasthesurfacesdeformandthepressureswilldecrease. Thecouplingofthesesolutionswouldlikelyshowadecreaseindeformation,compared totheresultsintable6.1. 56

PAGE 71

6.2TranslationAlongAxisPerpendiculartoLineofCenters Translationperpendiculartothelineofcenters,ormode2,introducesashearing motiontothesurfacesofthesphere.Theupperspherewasgivenatranslationvelocityof 1 m=s alongtheaxisperpendiculartothelineofcenters,withthelowerspherestationary. Thesamegeometryandmeshwasusedasinthemode1elasticitysimulationsinsection 6.1andthesimulationswereperformedwiththesamedecoupledmethod. Figure6.3:Exaggerateddeformedshapeofelasticspheresinmode2motion. = rad 10 2 E =0 : 2 GPa =0 : 4,scaledto10%ofthetotalgeometry.Deformationsarenormalized bysphereradius. Thedeformedshapeofthespheresareshowningure6.3,wheretheuppersphereis movingtowardtherightandthelowersphereisstationary.Theblackcirclesrepresent theoriginalshapeofthespheres.Fromtheplot,wecanseethattheupperspherehasa greaterdeformationinthedirectionoppositeofmotionthanthelower,stationarysphere. ThisisbecausetheStoke'sdragonthemovingsphereisaddedtotheshearstressfrom thepressuregeneratedintheinterstitialregion.Thereisalsoanupwardcompressionof thesurfaceofthemovingsphereintheinterstitialregion,fromthehigherpressuresinthe uid.Thelowerspherehasasmallamountofdeformation,causedbytheshearingeect fromthemovingsphere,butthereisnoStoke'sdragonthisspheresothedeformation islower. 57

PAGE 72

Table6.2:Maximumparticledeformationsformode2atseveralseparationdistances andparticleelasticities,correspondingtouidpressuresobtainedfromrigidspheresimulation,normalizedbysphereradius. Separation, Max:Deformation=ParticleRadius E =200GPa E =20GPa E =2GPa E =0 : 2GPa r= 10 3.75e-12 3.75e-11 3.75e-10 3.75e-9 r= 100 4.08e-12 4.08e-11 4.08e-10 4.08e-9 r= 1000 4.22e-12 4.22e-11 4.22e-10 4.22e-9 r= 10000 4.30e-12 4.30e-11 4.30e-10 4.30e-9 Themaximumdeformationoftheparticleswithdierentelasticitiesareshownin table6.2atseveralseparationdistances.Atallseparationdistances,astheelasticity isdecreasedbyanorderofmagnitude,thedeformationalsodecreasesbyanorderof magnitude.Themaximumdeformationwasfoundtobeindependentoftheseparation distanceforthismotion.Thesurfacedeformationsaresoinsignicantthatwecanexpect thehydrodynamicforcestobeapproximatelyequaltothoseofrigidspheresforthismode. 6.3RotationAboutLineofCenters Binaryspheresrotatingabouttheirlineofcentersisamotionthatwecanutilize axisymmetricgeometryfor,sinceCOMSOLallowsavelocitycomponentintheradial, or ,direction.Thismotionwasalsosolvedwithadecoupledapproach,wherethe creepingowmodulewasusedtosimulaterigidspheresrotatinginmineraloil.The pressuresfromthissimulationwereusedinaseparatesolidmechanicssimulationofelastic spheres.Thegeometryandmeshforthecreepingowsimulationsusedtoobtainuid pressureswereidenticaltothemeshformode3pressuredependentviscositysimulations. Thetwo-dimensionalsolidmechanicsmeshiscomposedof6522elementswith26628 degreesoffreedom.Thepressuredependentviscositysimulationsformode3inchapter 5determinedthatinsucientpressuresaregeneratedfromthismotiontoraisetheuid 58

PAGE 73

viscosity.Themotionhasnosingularityandwecanexpecttoseesimilarresultsinthe elasticitysimulations. Table6.3:Maximumparticledeformationsformode3atseveralseparationdistances andparticleelasticities,correspondingtouidpressuresobtainedfromrigidspheresimulation,normalizedbysphereradius. Separation, Max:Deformation=ParticleRadius E =200GPa E =20GPa E =2GPa E =0 : 2GPa r= 10 0 0 0 0 r= 100 0 0 0 0 r= 1000 0 0 0 0 r= 10000 0 0 0 0 Table6.3showsthatelasticspheresinmode3motionexperiencenodeformation. Thisisbecausethegaugepressureateverypointonthespheressurfaceisnegligible forthismotion.Sincethereisnodeformation,wecanexpecthydrodynamicforcesthat areequaltothoseofrigidspheres.Themotionforrotationabouttheparticles'lineof centershasnosingularity,andtheonlydeformationwouldbeindirectiontangentto thesurface,whichwouldnotcauseanincreaseinthehydrodynamictorqueforasteady statesolution. 6.4TranslationAlongLineofCenters Spherestranslatingalongtheirlineofcenterscanbesolvedwithatwo-dimensional axisymmetricgeometry.Ageometrysimilartothemode4pressuredependentviscosity simulationsinsection5.4wasused,excepttheproblemisfullycoupledandincludesa domainforthesphere,ratherthanonlyitsboundary.Theotherelasticitysimulations requiredthedecouplingofthecreepingowandsolidmechanicsphysics,buthavingthe two-dimensionalgeometryandasimplemotionallowedforasuccessfulsolutionofthe 59

PAGE 74

fullycoupledproblem.Inthiscase,theresultsfromthecreepingowsimulationareused asinputstothesolidmechanicsmodule,whichpassesitsresultsbacktothecreeping owmoduleforanotheriteration.Thisprocessisrepeateduntilthesolutionreachesa steadystate,atwhichtimetheresultsreectafully-coupledsimulation. Figure6.4:Meshusedforafully-coupledsimulationofmode4withelasticspheres. Aseparatemeshisconstructedforthecreepingowandsolidmechanicsregions, asseeningure6.4.Firsttheuidmeshwasconstructedidenticaltothemeshinthe mode4pressuredependentviscositysimulationsinsection5.4.Next,themeshforthe elasticspherewasconstructed.Theadjacentmeshessharenodesalongtheirboundaries, sothebottomofthespherewhereitmeetstheuidhasadensemeshandtheelements fartherfromthebottomgrowinsize.Thegreatestamountofstressanddeformationis wheretheuidpressuresarehighintheinterstitialregion,soadensemeshforthesphere isnecessaryhere.Therewere66231triangularelementsinthemeshforthesphere,in additiontotheelementsintheuidregion.Simulationsformode4withelasticparticles wereperformedforspheresofamaterialwithaPoisson'sratioof0.4atarangeofparticle elasticitiesandseparationdistances. 60

PAGE 75

Figure6.5:Exaggerateddeformationofmode4elasticparticlesofE=0.2GPa, = rad 10 3 inahighlyviscousuidwith =0 : 5,normalizedbysphereradius. Theareaofthehighestdeformationisshowningure6.5,inthehemisphereclosest totheinterstitialregion.Theareaofhighestpressuresintheuidactsontheadjacent sphereboundary,causingstrainingonthisareaofthesphere.Theblacklineinthegure showstheoriginalcurvatureofthesphere. Table6.4:Maximumnormalizedmode4particledeformationsatseveralseparation distancesandparticleelasticitiesinauidwith =0 : 5 Pa s .Thesimulationis fully-coupledbetweensolidmechanicsandcreepingow. Separation, Max:Deformation=ParticleRadius E =200GPa E =20GPa E =2GPa E =0 : 2GPa r= 10 4.56e-9 4.56e-8 4.56e-7 4.56e-6 r= 100 4.32e-8 4.31e-7 4.31e-6 4.28e-5 r= 1000 6.41e-7 6.33e-6 5.73e-5 3.58e-4 r= 10000 1.11e-5 5.73e-5 N/A N/A Themaximumdeformationsforelasticparticlesinahighlyviscousuidof =0 : 5 Pa s areshownintable6.4.Theabsentvaluesforthelowerelasticityparticlesarea resultoflargedeformationsontheparticlessurface.Thenumericalsolverwasunable toconvergeforthecaseswhendeformationswereverylarge.Thelargestnormalized 61

PAGE 76

deformationof3.58e-4isobservedat = rad 10 3 withaparticleelasticityof E =0 : 2 GPa Atthelargerseparationdistances,thedeformationincreasesbyanorderofmagnitude whentheelasticitydecreasesbyanorderofmagnitude,butthistrendabatesasthe separationdistancedecreases. Figure6.6:Maximumuidpressureintheinterstitialregionformode4motionofelastic spheresinmineraloil,atarangeofparticleseparationdistances.E=20GPa, =0.4. Asthespheresbecomeclosertogether,theuidpressurebetweenthemincreases rapidly,whichisdirectlyrelatedtotheforceontheparticleandtheparticles'deformation.Anorderofmagnitudeseparationdistancedecreasedoesnotresultinanorderof magnitudepressureincreaseatsmallseparationdistances,asweseeingure6.6.The gureshowsthemaximumpressureintheinterstitialregionforparticlesofE=20GPa atarangeofseparationdistances.Therelationshipbetweenpressureandbothparticle deformationandforceisnonlinearatsmallseparationdistances. 62

PAGE 77

Table6.5:Normalizedforcesonmode4elasticbinaryparticlesatarangeofseparation distancesandparticleelasticities. Separation, F Elastic =F Stokes E =200GPa E =20GPa E =2GPa E =0 : 2GPa Rigid r= 10 7.406 7.404 7.403 7.403 7.413 r= 100 53.019 53.016 52.983 52.657 53.422 r= 1000 503.460 499.290 464.414 324.559 504.434 r= 10000 4370.819 2792.957 N/A N/A 5004.90 Wecanseeintable6.5thattheforceonelasticparticlesdecreasessignicantlyfrom rigidparticles.Thisismostevidentatthelowestelasticityvalueof E =0 : 2 GPa ,which wecancomparewiththeforceonrigidspheres. Table6.6:Percentdecreaseinmode4forceonelasticparticles,comparedtorigidparticles. Separation, Percentdecreasefromrigidparticleforce E =200 GPa E =20 GPa E =2 GPa E =0 : 2 GPa r= 10 0.10 0.13 0.13 0.14 r= 100 0.76 0.76 0.82 1.43 r= 1000 0.19 1.02 7.93 35.66 r= 10000 14.51 44.20 N/A N/A Table6.6showsthatthesmallestseparationdistanceshavethegreatestdecreasein dragforce.Thisisbecausethesphereshavethelargestamountofsurfacedeformation attheseclosedistancesasaresultoftheextremeuidpressures. 63

PAGE 78

Figure6.7:Mode4hydrodynamicdragforceswithelasticparticles,wheretheseparation distancedoesnotaccountforsurfacedeformation. Itappearsfromgure6.7thattheparticleswithlowerelasticityhaveasignicantly lowerdragforceatsmallseparationdistances,whenthepressureisthehighest.Thisis becauseasthesurfacesofthespheresdeform,thegapgrowsandthepressurebetween thespheresisrelieved.Thechangeinseparationdistancefromtheparticledeformation isunaccountedforinthisgure. Figure6.8:Mode4hydrodynamicdragforceswithelasticparticles.Separationdistance hasbeenadjustedtoaccountforthesurfacedeformation. Wecanseefromgure6.8thatifthetruesurfaceseparationdistanceisplotted 64

PAGE 79

againstthedragforce,theresultsessentiallycollapsebackontothecurveforrigidspheres. Thisisdonebyadjustingtheseparationdistancetoaccountforthesurfacedeformation ofbothspheres.Thisleadstotheconclusionthatthedragforceisafunctionofthe separationdistanceofthespheresonly,andnotoftheparticleelasticity. 6.5Poisson'sRatio Whenamaterialiscompressedalonganaxis,ittendstoexpandintheothertwo orthogonaldirections.ThePoisson'sratioofanisotropic,linearelasticmaterialdescribes therelationshipbetweentransverseandaxialstrain,byhowmuchthematerialexpands relativetothedeformationintheaxisofcompressionortension.Materialsthatdemonstrateasignicanttransversedisplacementwhencompressedaxially,suchasrubbers, haveaPoisson'sratiocloseto0.5andareclassiedasnearlyincompressible.Materials thatexperienceanegligibletransversedisplacementcomparedtotheaxialdisplacement, suchascork,haveaPoisson'sratioclosetozero.Studyingtheeectsofthemodulusof elasticityisdirectlyrelatedtoPoisson'sratioandwewillexpandupontheinvestigation formode4wheretheparticledeformationswerethelargestofthemodes. Table6.7:Normalizedhydrodynamicforcesat = rad 10 3 andarangeofelasticitiesand Poisson'sratio's. Poisson'sratio, F=F Stokes E =200 GPa E =20 GPa E =2 GPa E =0 : 2 GPa 0 : 25 503.418 498.887 461.539 317.917 0 : 3 503.428 498.993 462.281 319.583 0 : 35 503.439 499.121 463.226 321.758 0 : 4 503.460 499.290 464.414 324.559 0 : 45 503.481 499.502 465.963 328.379 0 : 49 503.513 499.778 468.096 333.992 Table6.7showshowthedierentPoisson'sratiovaluesaectthehydrodynamic 65

PAGE 80

forcesateachmodulusofelasticity.WecanseethatasthePoisson'sratioincreases,the dragforceontheparticlesincreases.Thedierenceindragforcebetween =0 : 25and 0 : 49foranelasticityof0.2GPaisapproximately5%,whichisasmalldierencewhen comparingtotheeectofelasticmodulusonparticledragforcebeforeaccountingfor thesurfacedeformation.Althoughthisisbynomeansadirectcomparison,wecansay thatthePoisson'sratioofthematerialhaslittletonoeectonthedragforce,atleast attheseparationdistances = rad 10 3 .Asmallerseparationdistancewasnotsimulated becausethesolverwouldnotconvergeforthecombinationofseparationsbelow = rad 10 3 andelasticitiesbelow E =20GPa. Figure6.9:ComparisonofPoisson'sratiosvs.hydrodynamicforceforE=0.2GPa, = rad 10 3 ArangeofPoisson'sratiosareplottedagainstthenormalizedforceonparticleswith E =0 : 2GPainmode4motionat = rad 10 3 ingure6.9.Theplotshowsthatthegreater thePoisson'sratio,thegreaterthecorrespondingdragforceis. 66

PAGE 81

Table6.8:Percentdecreaseinmode4forceonelasticparticles,comparedtorigidparticles. Poisson'sratio, Percentdecreasefromrigidparticleforce E =200 GPa E =20 GPa E =2 GPa E =0 : 2 GPa 0 : 25 0.20 1.10 8.50 36.98 0 : 3 0.20 1.08 8.36 36.65 0 : 35 0.20 1.05 8.17 36.21 0 : 4 0.19 1.02 7.93 35.66 0 : 45 0.19 0.98 7.63 34.90 0 : 49 0.18 0.92 7.20 33.79 Thepercentdecreasesindragforcecomparedtoarigidparticleareshownintable 6.8forarangeofparticleelasticityandPoisson'sratiovalues.Foreachvalueof E ,the decreaseinforceislessforthehigherPoisson'sratiomaterials.Thisimpliesthatthe morecompressibleamaterialis,themorethedragforcewillbeaectedbytheparticles' changeinshape. Table6.9:Particledeformationsat = rad 1000 forarangeofelasticitiesandPoisson's ratio's. Poisson'sratio, Maximumsurfacedeformation/radius E =200 GPa E =20 GPa E =2 GPa E =0 : 2 GPa 0 : 25 6.55e-7 6.47e-6 5.85e-5 3.63e-4 0 : 3 6.41e-7 6.33e-6 5.73e-5 3.58e-4 0 : 35 6.23e-7 6.15e-6 5.58e-5 3.51e-4 0 : 4 6.00e-7 5.94e-6 5.40e-5 3.42e-4 0 : 45 5.72e-7 5.66e-6 5.17e-5 3.31e-4 0 : 49 5.37e-7 5.32e-6 4.88e-5 3.17e-4 67

PAGE 82

Table6.9containsthedatawecanusetoseeifthePoisson'sratiooftheparticles hasthesameeectastheelasticity.Themaximumdeformationforeachcombinationof E and areshowninthetable,whereweseethegreatestdeformationforthelowest E andlowest .Thiscombinationofmaterialpropertiescorrespondtothesoftest,most compressiblematerials. Figure6.10:ComparisonofPoisson'sratiosof0.25and0.49forseveralelasticities.Left: deformationunaccountedforinseparation,right:deformationaccountedfor. Theplotsfor =0 : 25and =0 : 49areshownforelasticitiesof20GPa,2GPa,and 0.2GPa,ingure6.10.Ontheleftplot,wecanseethattheforceappearslowerforthe elasticmaterials,whichiswhatwasseeninsection6.4.Italsoappearsthatthelower thePoisson'sratio,thelowertheforce.Theparticledeformationisaccountedforinthe separationdistanceontherightplotandweseethattheforcesforeverycombinationof E and agreewiththeforceonarigidparticle.ParticleswithalowPoisson'sratioof 0.25havealargerdeformationthanthosewithaPoisson'sratioof0.49andthepoint hasalargerrightwardshiftthanthehigherPoisson'sratioof0.49ontherightplot. 68

PAGE 83

Figure6.11:Comparisonof =0.25and =0.49forE=0.2GPa.Left:deformation unaccountedforinseparation,right:deformationaccountedfor. TheforcesonparticleswithPoisson'sratiosof0.25and0.49areshowningure6.11 foraparticleelasticityof0.2GPa.Theleftplotshowsthatastheseparationdistance decreases,thelowerPoisson'sratioparticleshaveaslightlylowerdragforce.Theless compressiblethematerialis,forinstanceamaterialwithaPoisson'sratioof0.49,the higherthedragforceappearstobe.Thisisbecausethedeformationoftheparticles' surfaceisgreatersincethematerialismorecompressible,givingmorespacefortheuid intheinterstitialregiontorelieveitselfofhighpressures.Intherightplot,weadjust theseparationdistancebytheparticles'surfacedeformationandthecurvesessentially collapseontooneanother,whichissimilartotheresultintheelasticitysimulationsin theprevioussections.Thisbehaviorisconsistentwiththeobservationthatthedrag forceisafunctionoftheactualseparationdistancebetweentheparticles'surfaces,and notoftheirPoisson'sratio. 69

PAGE 84

Pressure-dependentruidviscosityandparticleelasticitybothhaveuniqueaectsonthehydrodynamicforcesandtorquesonbinaryspheresinsuspensionsystems.Fluidviscositycanincreasesignifcantlyatenhancedpressures,whichcanaectthemotionofbinaryspheresinviscousruid.Wesawthatbinaryspherestranslatingalongtheirlineofcentershavethelargestincreaseinhydrodynamicforceduetopressure-dependentviscosity,comparedtotheothermodesofmotion.Anapproximate39%increaseinhydrodynamicforcewasobservedforparticlesinthismotionatavelocityof2.15m sinmineraloil.Themajorityoftheforceisfromtheextremeruidpressuresintheinterstitialregionbetweenthespheres.Theotherthreemodesofmotiondidnotsignif-icantlyincreasethehydrodynamicforceortorquewhenintroducingpressure-dependentviscositybecausethepressuresgeneratedareinsucienttoincreasetheruid'sviscosity,signifcantly.Theactualmotionofbinaryparticlesinsuspensionsisacombinationofthefourindividualmotions,somode4wascombinedwiththeotherthreemodestodetermineifthepresenceoflargepressuresaectedthehydrodynamicforceortorqueinmodes1-3.Themode4particlevelocitywas2m sinthesesimulations.Theincreaseintorqueonmode1particlesrotatingat1rad swas5.75%andtheforceincreaseonmode2particlestranslatingat1m swas9.62%.Mode3particlesrotatingat1rad sexperiencednotorqueincreasewhencombinedwithmode4motionbecausethemode3componentsofvelocityarenegligibleintheinterstitialregionwherethehighpressuresarepresent.Thepressure-dependentviscositysimulationsshowthatifparticlevelocitiesaresucient,thehydrodynamicforcesimposedcanbesignifcantlyperturbedbyaruidwithpressure-dependentviscosity,comparedtoconstantviscosity.Particleelasticitywassimulatedforbinarysphericalparticlesinviscousruidswitheachofthefourmodesofmotion.Particlesinmode4motionwerefoundtohavethe70 CHAPTERVIICONCLUSION

PAGE 85

largestdeformationscomparedtotheothermodes.Asignicantdecreaseindragforce wasobservedformode4elasticparticles,comparedtorigidparticles,atcomparable center-to-centerparticledistances.Thelargestforcedecreaseof44%inmode4wasobservedataseparationdistanceof rad 10 4 andaparticleelasticityof20 GPa .Thesimulations forparticleelasticitiesbelow20 GPa didnotconvergeduetolargemeshdeformations. Althougha44%decreasefromtheforceonrigidparticlesappearstobesignicant,a shiftoftheseparationdistancebytheamountofparticledeformationshowsotherwise. Aftermeasuringthesurfacedeformationoftheparticlesandaccountingforthisinthe separationdistance,thegapsizeincreasesandthesepointsshiftontotheoriginalcurve forarigidparticle.Theseresultsdeterminethatparticleelasticitydoesnothaveadirect aectonthehydrodynamicforcebecauseelasticandrigidparticlesincurthesameforces whentheelasticparticleseparationisadjustedtoaccountfordeformation. ThesimulationsusedtoinvestigatetheaectofparticleelasticityassumedaPoisson's ratiovalueof0.4.Becauseofitsrelationshipwiththemodulusofelasticity,thePoisson's ratiooftheparticleswasalsostudied.ArangeofPoisson'sratiosfrom =0.25through =0.49wereconsideredforeachelasticmodulusvalue.Theforcesonparticlesdidnot changesignicantlyatdierentvaluesof foreachvalueof E .Similartotheresults forparticleelasticity,iftheseparationdistanceisadjustedtoaccountforthesurface deformation,theforcesessentiallyfallbackonthecurveforarigidparticle.Wecan concludethatthePoisson'sratioofamaterialdoesnothaveasignicantaectonthe hydrodynamicforcesortorquesonbinaryparticlesinviscousuids. 71

PAGE 86

REFERENCES [1] E.P.Ascoli,D.S.Dandy,L.G.Leal LowReynoldsNumberHydrodynamic InteractionofaSolidParticleWithaPlanarWall J.Numer.Meth.Fluids 9,651688. [2] S.Bair High-PressureRheologyforQuantitativeElastohydrodynamics ,Elsevier, Amsterdam,TheNetherlands,2007. [3] C.Barus Isothermals,IsopiesticsandIsometricsRelativetoViscosity Am.J. Science ,ThirdSeries,Vol.XLV,No.266,1893. [4] A.W.Batchelor,G.W.Stachowiak EngineeringTribology ,secondedition, Butterworth-Heinemann,Oxford,UK,2001. [5] G.K.Batchelor,J.T.Green Thedeterminationofthebulkstressinasuspensionofsphericalparticlestoorder c 2 J.FluidMech. 56 ,401-427. [6] R.B.Bird,E.N.Lightfoot,W.E.Stewart TransportPhenomena ,second edition,JohnWiley&Sons,NewYork,N.Y.,2001. [7] H.Brenner Theslowmotionofaspherethroughaviscousuidtowardsaplane surface ,.Chem.Eng.Sci.16,242-251. [8] H.Brenner,J.Happel LowReynoldsnumberhydrodynamics ,Springer,The Netherlands,1983. [9] S.R.Challa,F.V.Swol Molecularsimulationsoflubricationandsolvation forces ,Phys.Rev.E736,016306. [10] J.S.Chong,E.B.Christiansen,A.D.Baer RheologyofConcentratedSuspensions J.AppliedPolymerSci. 15 ,2997-2021. [11] G.D'Avino,F.Snijkers,R.Pasquino,etal. Migrationofaspheresuspended inviscoelasticliquidsinCouetteow:experimentsandsimulations ,RheologicaActa 51:215. [12] A.Einstein EineneueBestimmungderMolekldimensionen ,.Annalender Physik.19:289. [13] G.E.Fernandez,P.M.Carrica,D.A.Drew BinaryInteractionForce onSpheresI:LowReynoldsNumberFlow Chem.Eng.Comm. [14] E.Guth,R.Simha UntersuchungenberdieViskosittvonSuspensionenundLsungen.3.berdieViskosittvonKugelsuspensionen ,.KolloidZ.74:266. [15] W.Habchi AFull-SystemFiniteElementApproachtoElastohydrodynamicLubricationProblems:ApplicationtoUltra-Low-ViscosityFluids .PhDthesis,INSAde Lyon,Villeurbanne,France;2008. 72

PAGE 87

[16] J.S.Halow,G.B.Wills ExperimentalObservationsofSphereMigrationinCouetteSystems ,.Ind.Eng.Chem.Fundamen.9:603-607 [17] B.J.Hamrock FundamentalsofFluidFilmLubrication ,NasaReferencePublication1255,55-59. [18] M.S.Ingber,S.Feng,A.L.Graham,H.Brenner Theanalysisofselfdiusionandmigrationofroughspheresinnonlinearshearowusingatractioncorrectedboundaryelementmethod ,J.FluidMech. 598 ,267-292. [19] D.J.Jeffrey,Y.Onishi Theforcesandcouplesactingontwonearlytouching spheresinlow-Reynolds-numberow Z.Angew.Math.Phys. 35 ,634-641. [20] I.M.Krieger,T.J.Dougherty Amechanismfornon-Newtonianowinsuspensionsofrigidspheres. Trans.Soc.Rheol. III,137-152. [21] S.R.Majumdar,M.E.O'Neill Asymmetricalslowviscousuidmotionscaused bytherotationortranslationoftwospheres.PartI.Thedeterminationofexact solutionsforanyvaluesoftheratioofradiiandseparationparameters Z.Angew. Math.Phys. 21 ,164-179. [22] O.Reynolds OnthetheoryoflubricationanditsapplicationtoMr.Beauchamp Tower'sexperiments,includinganexperimentaldeterminationoftheviscosityof oliveoil ,Phil.Trans.RoyalSoc., 177 ,157-234. [23] C.J.A.Roelands CorrelationalAspectsoftheViscosity-Temperature-Pressure RelationshipofLubricatingOils ,Ph.D.Thesis,UniversityofTechnology,Delft,1966. [24] H.Sohn,S.Koo PermeabilityofViscousFlowThroughPackedBedofBidisperse HardSpheres ,.KoreanChemicalEngineeringResearch.50:66-71. [25] N.Tetlow,A.L.Graham,M.S.Ingber,S.R.Subia,A.L.Mondy,S.A. Altobelli ,1998. ParticlemigrationinaCouetteapparatus:Experimentandmodeling.J.Rheol. 42,307-327. [26] D.G.Thomas Transportcharacteristicsofsuspension:VIII.AnoteontheviscosityofNewtoniansuspensionsofuniformsphericalparticles ,.J.ColloidSci. 20:267 73

PAGE 88

APPENDIX A.MeshConvergenceStudies Convergencestudiesareimportantinnumericalsimulationstoguaranteethemesh constructionisadequateforthenumericalsimulationtoproducereliableresultsforthe physicsinthesimulation.Aconvergencestudydoesnotguaranteethattheresultsare correct,butratherthatconvergenceisbeingobtainedandtheresultsbetweensuccessive simulationsofvaryingbutsimilarmeshesarenotchangingmorethanapredetermined threshold.Anexampleisforthedeformationofabeam,whereanextremelycoarsemesh isrstdened.Asimulationusingthismeshwillyieldaninitialresultwhichwearenot sureisconvergedornot.Themeshisthenmadeslightlydenserandthenewresultmight bemuchdierentthanthepreviousresult. Thisiterativeprocessiscontinueduntilthemeshisneenoughwherethedierence betweensolutionsisbelowapproximately1%.Atthispoint,wehavedeterminedthe meshisconvergedbelowthethreshold.Aconvergencestudyisprovidedforeachsimulationperformedinthisthesis,atthesmallestseparationdistanceof = rad 10 4 ViscositySimulations Mode1 Rotationaboutaxisperpendiculartolineofcenters TableA.1:Convergencestudyformode1three-dimensionalmeshat = rad 10 4 basedon changeinforcebetweensubsequentmeshes. No.elem. Deg.offreedom Torque[N-m] PercentTchange 264996 2890202 4.104 N/A 298614 3267810 4.249 3.53 315246 3449502 4.326 1.81 354150 3875544 4.357 0.71 74

PAGE 89

Increasingthetotalnumberofelementsinthemeshallowedconvergencetobeobtainedat354150elements,withonlya0.71%changefromthepreviousmesh. Mode2 Translationalongaxisperpendiculartolineofcenters TableA.2:Convergencestudyformode2three-dimensionalmeshat = rad 10 4 basedon changeinforcebetweensubsequentmeshes. No.elem. Deg.offreedom Force[N] PercentFchange 264996 2890202 2.371 N/A 298614 3267810 2.449 3.18 315246 3449502 2.498 1.96 354150 3875544 2.515 0.68 Themode2meshconvergenceissimilartothemode2results.Convergencewas reachedatthesameamountofelementsforbothmodes. Mode3 Rotationaboutlineofcenters TableA.3:Convergencestudyformode3two-dimensionalmeshat = rad 10 4 basedon changeinforcebetweensubsequentmeshes. No.elem. Deg.offreedom Torque[N-m] PercentTchange 4586 32522 1.941 N/A 16551 112643 1.801 6.96 36962 201664 1.793 0.55 44862 210632 1.791 0.11 Mode4 Translationalonglineofcenters Inmode4,themajorityoftheforceisfromthehighuidpressuresintheinterstitial region.Thebulkuidmeshcanberelativelycoarse,withahighdensitymeshinthe interstitialregion. 75

PAGE 90

TableA.4:Convergencestudyformode4two-dimensionalmeshat = rad 10 4 basedon changeinforcebetweensubsequentmeshes. No.elem. Deg.offreedom Force[N] PercentFchange 30195 141010 13956 N/A 58665 274553 13837 0.85 62143 290344 13814 0.16 65284 304581 13814 0.00 Alargedierenceisnotseenwhenincreasingthenumberofelementsinthemode 4meshedbecausetheinterstitialelementsareguaranteedtobeneenoughtotinthe gap.Ifthisconditionwassatised,themeshwassatisfactorysincemostoftheforce cameformtheinterstitialregionpressures. ElasticitySimulations Mode1and2 Thegeometryandmeshwasidenticalformode1and2elasticitysimulations.Once thesimulationwasperformedfortheuidportionoftheproblem,theresultingpressures weremappedontoelasticspheresinanothersolidmechanicssimulation.Themeshforthe uidportionwasidenticaltothemeshforthemode1and2constantviscositynumerical simulations.Thesamesolidmechanicsmeshwasusedinmodes1and2,basedonMode 1convergenceresults. 76

PAGE 91

TableA.5:Convergencestudyformode1three-dimensionalsolidmechanicsmeshat = rad 10 4 andE=0.2GPabasedonchangeindeformationbetweensubsequentmeshes. No.elem. Deg.offreedom Def Norm: Percentdef.change 26442 112537 1.03e-6 N/A 31468 135426 1.09e-6 5.82 76998 331168 1.12e-6 2.75 143676 600459 1.13e-6 0.89 Thedensityofthemeshintheinterstitialregionisimportanttodeterminethecorrect displacementformodes1and2.Thepressuresinthisareaarelarge,somappingthe pressuresontoadensermeshinthisareaisimportant.Arelativelyuniformmeshwas usedforthewholesphere,sotheamountofelementsislargetoguaranteeahighenough densityintheinterstitialregion. Mode3 TableA.6:Convergencestudyformode3two-dimensionalsolidmechanicsmeshat = rad 10 4 andE=0.2GPabasedonchangeindeformationbetweensubsequentmeshes. No.elem. Deg.offreedom Def Norm: Percentdef.change 976 4188 0 N/A 1480 6204 0 0 6522 26628 0 0 Mode3experiencesnodeformationsinceuidpressuresfromthecreepingowsimulationsareverylow,sothedensityofthemeshwasnotcritical. Mode4 Mode4usedafullycoupledapproachwherethemeshofthesolidmechanicsand creepingowmoduleswereadjacentinthesamesimulation. 77

PAGE 92

TableA.7:Convergencestudyformode4two-dimensionalfully-coupledsimulationmesh at = rad 10 4 andE=20GPabasedonchangeinforcebetweensubsequentmeshes. No.solid elem. Totaldeg. offreedom F F Stokes PercentF change 41248 235938 2652.2 N/A 58912 338744 2720.6 2.59 62985 359014 2768.10 1.75 66231 380184 2792.95 0.89 Thedensityoftheuidmeshislimitedbythedensityoftheuidmeshwhichwas constructedrst.Thesolidmeshbeginswiththesizeofelementsfromtheadjacent uidmeshandexpandsthemupwardawayfromtheinterstitialregion.Themeshisnot abletogrowtoocoarsebecauseofthis,sotheforcesforallmeshesarenotfarfromthe convergedmesh. 78