UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Fibre motion in shear flow Wherrett, Geoffrey 1996

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1996-0375.pdf [ 5.93MB ]
JSON: 831-1.0080821.json
JSON-LD: 831-1.0080821-ld.json
RDF/XML (Pretty): 831-1.0080821-rdf.xml
RDF/JSON: 831-1.0080821-rdf.json
Turtle: 831-1.0080821-turtle.txt
N-Triples: 831-1.0080821-rdf-ntriples.txt
Original Record: 831-1.0080821-source.json
Full Text

Full Text

FIBRE MOTION IN SHEAR FLOW By Geoffrey Wherrett B. Eng. (Mechanical) McGill University , 1992  A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F T H E REQUIREMENTS F O R T H E D E G R E E O F M A S T E R OF A P P L I E D SCIENCE  in T H E F A C U L T Y OF G R A D U A T E STUDIES MECHANICAL ENGINEERING  We accept this thesis as conforming to the required standard  T H E UNIVERSITY OF BRITISH C O L U M B I A  May 1996 © Geoffrey Wherrett, 1996  In  presenting  degree  at the  this  thesis  in  University of  partial  fulfilment  of  British Columbia, I agree  freely available for reference and study. I further copying of  this thesis for  department  or  by  his  or  the  representatives.  that the  for  an advanced  Library shall make  agree that permission for  scholariy purposes may be her  requirements  It  is  granted  by the  understood  that  it  extensive  head of copying  my or  publication of this thesis for financial gain shall not be allowed without my written permission.  Department The University of British Columbia Vancouver, Canada  DE-6 (2/88)  Abstract  The motion of a particle at low relative velocity to a shearing fluid is of importance in many areas of engineering. One such application is in the pulp and paper industry where pulp fibres may be separated by means of the fluid forces acting upon them. This process, known as fibre fractionation, divides a pulp fibre stock into fractions with higher percentages of fibres with certain properties. Experimental studies have shown that the motion of rigidfibresin shear may be described with an analytical model. The deformation that occurs in flexible fibres makes it difficult to model their motion in shear. In this thesis, a numerical model of a single fibre was formulated for both rigid and flexible motion. This research is part of a larger project aimed at modelling the hydrodynamic devices which could achieve efficient fractionation. Thefibrewas modelled as a cylindrical rod made up of elements of equal length and diameter. The inertial and structural properties of the model were based on the cylindrical geometry. The hydrodynamic force on each element was calculated for a sphere of the same diameter as the fibre. Equations of motion were formed for each element from the fluid and fibre forces as well as a continuity condition at the boundary of adjacent elements. These equations were simultaneously solved at each time step. The motion of thefibrewas then followed by advancing the time step and recording the location of each element. Three series of trials were conducted with the model. The first determined if the model accurately simulated the motion of a rigid fibre. The second determined the effect of bending stiffness on the fibre's motion. The final series determined the effect of a localized weakness in the fibre. Trials were conducted with solidfibresof aspect ratio 10 and 60 in a Couette flow. The initial position of thefibrewas aligned with the y axis either centred on or extending up from the origin. The first position resulted in purely rotational motion while the second position also included  ii  translation. Initial trials indicated that the model performed in accordance with previous results and within reasonable agreement with theory for both the rotational and translational motion of a rigid fibre. The second series of test was conducted with a 60 aspect ratio fibre of varying stiffness. This aspect ratio was chosen as it falls in the range of typical pulp fibres. The stiffness was measured by the ratio of elastic to shear forces, k i = e  E/TJG.  Here E is the elastic modulus of  the fibre, n is the absolute viscosity of the fluid and G is the rate of shear. The model performed rigidly when k i > 2 X 10 . Some deformation was seen in the fibre above 5 X 10 but the period 6  5  e  of rotation was nearly constant at a value close to that from the rigid case. As the stiffness was reduced below this value the fibre deformed significantly during rotation and its period of rotation decreased. The third series of trials was conducted with the same 60 aspect ratio fibre with k i — 5x 10 . 5  e  A localized weakness, at which the stiffness of the fibre was reduced by 80%, was placed in the fibre at 5 locations along its length. The results showed that the period of rotation did not change from the similar case in the previous series. However the fibre rotated in an asymmetric manner with one end of the fibre bending first. This differed from the previous series where both ends of the fibre behaved symmetrically. A fibre bending number B was derived to determine the type of motion in the fibre model depending on the fluid and fibre properties. This number is made up of three terms. The first is k i, the second is a function of the fibre's hollowness, and the third involves the aspect ratio. e  Comparisons with results showed that there existed a critical value for B, above which the fibre behaved rigidly. This is significant in that a simple rigid fibre model may be used when B is above the critical value. The trials conducted with the present model indicated that the model performed in a reasonable manner with respect to previous experimental and theoretical results. Several improvements could be made to the model. In particular, elements of larger aspect ratio would reduce the computational time.  iii  Table of Contents  Abstract  i i  Table of Contents  iv  List of Figures  vi  List of Tables  viii  List of Symbols  xi  Acknowledgement 1  Introduction  1  1.1  Motivation  1  1.2  Pulp Properties  3  1.3  Pulp Processing  8  1.4  Literature Review  14  1.4.1  Rigid Fibre Motion  14  1.4.2  Flexible Fibre Motion  17  1.4.3  Computer Simulation of Particle Dynamics  21  1.5 2  xii  Scope and Aims of This Thesis  22  Numerical Procedure  23  2.1  24  Model Formulation 2.1.1  Fibre Structural Model  2.1.2  Hydrodynamic Force Model  2.2  Method  2.3  Application  .  24 29 33  . . . ."  37  iv  3  Non-dimensional Groups  37  2.3.2  Computational Algorithm  41  Results and Discussion  44  3.1  Introduction  44  3.2  Results for the 10 Element Model  48  3.2.1  48  3.3  4  2.3.1  Rigid Body Motion - Series 1  Results for 60 Element Model  54  3.3.1  Bending Deformation - Series 3  55  3.3.2  Localised Wall Weakness - Series 4  64  3.4  Discussion  66  3.5  Recomendations for Improvements to the Model  73  Conclusion  75  A Inertial Properties of a Fibre Element  77  B  79  Structural Properties of the Fibre Model  C Dimensionless Parameters and Equations of Motion  81  D Velocity Boundary Condition  85  E  Angular Correction Routine versus Re — Series 2  F  Code listing  p  88 90  Bibliography  96  v  List of Figures  1.1  Wall Thickness  4  1.2  Fibre Structure  7  1.3  Schematic diagram of pressure screen operation  11  1.4  Schematic diagram of a hydrocyclone  13  1.5  Orientation angle of Rotating Particles  15  1.6  Equivalent Ellipsoidal Aspect Ratio, r*  16  1.7  Typical rotational orbits of  2.8  Schematic diagram of a continuous fibre and its representation by a series of  2.9  fibres  19  cylindrical elements  25  Configuration of bonded elements when stretching  26  2.10 Configuration of bonded elements when bending  27  2.11 Shear Flow  30  2.12 Element Geometries: spherical and cylindrical  32  2.13 Hydrodynamic and Tangential Forces  34  2.14 Schematic Diagram showing nonslip condition. 6{ + dj = 2a^-  35  3.15 Initial position of fibres . .  45  3.16 Results - 10 elements, simple rotation  50  3.17 Trial 1(c): Rigid 10 Element Model, Rotation and Translation  53  3.18 Orientation of rigid 60 element model  57  3.19 Fibre Deflection at <j) = 180° for Trials 3(a) to 3 (c)  58  3.20 Fibre Motion for Trials 3(a) to 3 (c)  59  3.21 Fibre Motion for Flexible Trials  61  vi  3.22 Stiffness vs Reynolds Number  62  3.23 Fibre Orbits for Series 3  63  3.24 Fibre Motion for Trials 4(a), 4(b), and 4(c)  65  3.25 a) TG/(TG)  rigid  versus B (Eq. 3.42), b) TG/{TG)  rigid  E.26 Rigid 10 Element Fibre, T G vs Re  versus B (Eq. 3.43)  . . .  69 89  p  vii  List of Tables  1.1  Variation of Dry Relative Density for Douglas Fir trees in British Columbia . . .  5  2.2  Fibre and Flow Parameters  39  2.3  Dimensionless Groups  40  3.4  Constants and Parameters  46  3.5  Constant Values used In Series 1  48  3.6  Results for Trial 1(a)  49  3.7  Results for Trial 1(b)  51  3.8  Constant Values used In Series 3 and 4  54  3.9  Series 3 Elastic Parameters  55  3.10 Results for Trial 3(a)  56  3.11 Constant Values used In Series 4  64  3.12 Results for a flexible 10 element fibre and the two calculated bending ratios . . .  68  3.13 Results from Series 3 and the two calculated bending ratios  68  3.14 Variation of fibre parameters for softwood pulp in British Columbia  71  E.15 Results from Series 2 .  88  viii  LIST OF S Y M B O L S  A.  cross sectional area of fibre wall  Af  cross sectional area of fluid portion of element  B  fibre bending number  E  Young's modulus  ca  C t  rotational constant  CtT  translational constant  ph  hydrodynamic force  ps  spring force  G  shear rate  I> Ics  second moment of area of fibre cross-section  h  element mass moment of inertia  If  mass moment of inertia of fluid portion of element  Ifw  mass moment of inertia of fibre wall portion of element  N  number of elements  Re  Reynolds Number  TO  Re  particle Reynolds Number  s  element specific gravity  Sfw  specific gravity of fibre wall  p  e  rpb  bending torque hydrodynamic torque  u  fluid velocity in x direction  V  fluid velocity in y direction  w  fluid velocity in z direction  ix  a*  wall thickness ratio  ^  radius of inner wall  a  radius of outer wall  di  diameter of inner wall  d  diameter of outer wall  /'  tangential force due to non slip condition  i  element i, acting on element i  ij  acting on element i from element j  ki  elastic constant  k  fibre spring constant  kb  fibre bending constant  I  distance between element centres  l  element length (l = 2a ), equilibrium distance between element centres  Q  a  e  a  0  Q  0  mj  mass of fibre wall portion of element  m/  mass of fluid portion of element  riij  normal pointing from the centre of element i to element j  v  generalized global coordinate  Ti  position vector to centre of element i  r  cylindrical aspect ratio  w  c  r  ellipsoidal aspect ratio  T*  equivalent ellipsoidal aspect ratio  t  time  t*  dimensionless time (tG)  t  torque due to angular non slip condition  u  velocity in x direction  v  velocity in y direction  x  displacement in x direction  y  displacement in y direction  e  9  x  ctij  angle formed by the normal n^-  St  dimensionless time step  e  strain  j]  fluid  u>  viscosity  rotational velocity about z axis  Q  fluid  rotational velocity  cf>  angular position of rigid fibre  p  element density  e  Pf  fibre wall density  Pf  fluid  w  density  a  stress  6  angular displacement about z axis  Oh  bending angle between neighbouring elements  0bo  equilibrium bending angle  (t  translational friction constant  £  rotational friction constant  r  xi  Acknowledgement  I would like to thank my advisors, Dr. Ian Gartshore and Dr. Martha Salcudean. Their constructive input and suggestions were greatly appreciated and made this thesis possible. In addition to my supervisors many individuals in the our research group assisted me in various aspects of the work. Many thanks to Dr. Pingfan He for his help in formulating the numerical code. Mike Savage's sage advice was appreciated throughout. I would like to thank the rest of the gang in the Rusty Hut for their continuing support. Special thanks to Gerry Rohling and Alan Steeves from the Computer Lab in Mechanical Engineering who tirelessly answered my questions and increased my disk quota. James Olson and Dr. Robert Gooding from PAPRICAN clarified many aspects of the fractionation process for me. Finally, I would like to acknowledge Michele who helped me get through this. Thanks.  xii  Chapter 1  Introduction  1.1  Motivation  The challenges of the pulp and paper industry serve as the impetus for research and development in many areas of applied science. Mechanical and chemical engineers, in particular, are required to design and implement the processes which create paper products from virgin timber and recylcled fibres. This role is vital to Canada's economy in which pulp and paper is the largest export product. The global market for pulp and paper products continues to grow. Research has shown that demand increases linearly with population growth. Since most of this population growth is occurring in the developing world, this will be the market for Canadian products. Fortunately, British Columbia has close access to the burgeoning economies of Asia where the largest demand is occurring. The fibre supply needed to satisfy the growing demand for pulp and paper products will come largely from three areas: the boreal forest which stretches across Canada, Scandanavia, and into Siberia; plantations in tropical and sub-tropical regions; and recycled fibres from large urban areas. British Columbia's traditional advantage in this area has been the relatively easy access to high quality softwood trees on the Coast and in the Interior. This advantage is being eroded due to restrictions on logging old growth forests and slow natural regeneration. The high growth rates of timber in tropical climates allows intensive management of plantations where trees are harvested at a rate several times faster than in the boreal forest. Pulp mills built in these regions also have the advantage of access to a pool of low cost labour. This leaves British Columbia with just one of its traditional advantages: low cost power from hydroelectricity.  1  Chapter 1. Introduction  2  W h a t can B r i t i s h C o l u m b i a ' s pulp and paper i n d u s t r y do to improve its position i n the w o r l d wide market? T h e solution that has been put forward by m a n y can be s u m m a r i z e d by the t e r m  value added. T h i s approach suggests replacing p r o d u c t i o n of lower value products  (ie. newsprint) w i t h higher value products (ie. coated paper). T h i s approach generally means m o v i n g i n t o a narrower market where a new volume of product could lower prices a n d reduce profit m a r g i n s . T o succeed i n the new market a product must have a unique appeal w h i c h can be m a r k e t e d or be manufactured w i t h a lower cost technology t h a n competing p r o d u c t s . B r i t i s h C o l u m b i a has an advantage over m a n y areas of the w o r l d . T h e coniferous forests w h i c h cover m u c h of the province produce pulp fibres which are l o n g , slender, a n d flexible. These properties allow the fibres to conform and b o n d well together.  M a n y specialty paper  products require this type of fibre to increase their strength. However, these " g o o d " fibres are m i x e d w i t h less desirable fibre i n the tree. A technology which can separate fibres w o u l d be of great benefit i n creating unique paper products. T h e m a t u r i n g of the recycling industry has also resulted i n changes to the fibre supply. A larger p r o p o r t i o n of fibres have already been through the various m a n u f a c t u r i n g processes at least once. A technology which can separate well processed fibres f r o m v i r g i n fibres could lower refining costs. T h e basis for m o d e r n pulp and paper manufacturing has existed for over 100 years. D u r i n g this t i m e , efforts have focussed on i m p r o v i n g existing processes t h r o u g h analysis and experim e n t a t i o n i n b o t h academia and industry. However there is considerable difficulty i n modelling the complex physical processes involved. In recent years, the advent of high speed computers has allowed m a n y processes to be modelled numerically. T h i s thesis is a first step i n modelling the m o t i o n of a fibre i n a m o v i n g fluid. U s i n g the experience gained f r o m this m o d e l , it m a y be possible to use this t o o l along w i t h fluid m o t i o n models to design devices which separate fibres w i t h different properties. T h i s process is k n o w n as fibre f r a c t i o n a t i o n .  Chapter 1. Introduction  1.2  3  Pulp Properties  The qualities of a paper product are determined by the properties of the pulp fibres from which it is manufactured. The nemesis of consistent paper quality is the variability that occurs in wood fibres due to the natural growing process. This variability between fibres is seen at differing scales. At the largest scale, differences are dependent on the wood structure. At the next level, the type of cells and the cellular structure vary. Finally, the chemical composition of the cell wall becomes important.  Structural Variation in Wood Fibres The primary source for pulpfibresin Canada is softwood or coniferous trees. By far the largest constituent of conifers are the long thin tracheid cells (95% by volume, 99% by weight). These fibres serve as the structural units of the tree. The strength and stiffness required for this role makes them ideal as pulp fibres. Deciduous or hardwood trees are also made up primarily of tracheids (65% by volume, 86% by weight). However, the average length of a hardwood tracheid is 1/4 that of a conifer which inhibits their use in paper products. Other cell types exist in all species of trees. For example, vessels transport nutrients in the form of sap throughout the tree. They have a round structure which is illsuited for papermaking. For this reason, present cleaning technologies are designed to remove vessels from among the tracheids in pulp. The shape of softwood tracheid cells can best be thought of as a hollow cylinder. Average mature fibres are 4 mm long with a diameter of 40 microns [19]. Mature fibres (older than 60 years) are, on average, twice as long as young fibres (less than 20 years). This has two implications. Fibres located lower in the tree (the trunk) are more mature and hence longer than those in the crown. Short rotation tree farming, a practice common in tropic and subtropic climates, will be unable to produce the same quality of fibres as those harvested from naturally regenerating forests. The growth cycle of trees is evident from the annual rings. These rings are visible due to the abrupt change in the thickness of the cell walls. This is shown in Fig. 1.1. Earlywood  Chapter 1. Introduction  4  (springwood) fibres have properties which are quite different than that of latewood (summerwood) fibres. The overall dimensions of earlywood fibres are, on average, larger than latewood fibres. Earlywood fibres also have a much thinner cell wall than latewood fibres. As a result, their dry density may be only one third to one half that of latewood [19]. The properties of the whole fibre depend on the cell wall thickness. For this reason, several measures related to the cell wall are quoted in pulp and paper literature. The pulp density is a good measure of the thickness of the cell wall. This is the ratio of the oven dry mass to the volume of a fibre. The volume may be measured from the green or oven dry value. The unit of measurement used is relative density which is defined as the pulp density divided by the density of water.  Figure 1.1: A cross section of Douglas Fir showing the abrupt change in cell wall thickness at the annual ring and the earlywood-latewood differentiation. Taken from Fig. 2-9 in [19].  Chapter 1. Introduction  5  The cell wall has a density of approximately 1.5 grams per cubic centimetre. This value can be used along with the relative density to determine the proportion of the total volume occupied by the cell wall. Consider the example of dry spruce wood which has a relative density of approximately 0.375. The cell cavity is filled with air when oven dried and assumed to have negligible mass.  cell wall volume _ density of fibre _ 0.375 g/cm _ 1 total volume density of cell wall 1.5 g/cm 4 3  3  Two other measures are frequently used to describe the properties of pulp fibres. The specific surface is a measure of the surface area per unit dry mass of the particle. Similarily, the specific volume is the volume of a particle per unit dry mass of solids. These properties are measured in a statistical manner in order to find the proportions of the fibres which fall into acceptable ranges for a desired product. Variable Between Sites One Site Single Tree Growth Ring  Range 0.44 - 0.47 0.42 - 0.48 0.36 - 0.50 0.28 - 0.70  Table 1.1: Variation of Dry Relative Density for Douglas Fir trees in British Columbia  Table 1.1 shows the variation of dry relative density in Douglas Fir trees in British Columbia. The small variation between averages of several trees from different sites may be attributed to different climatic conditions at different sites. The variation between whole tree averages is due to the relative age of trees in one site. The variation between different parts of the tree shown in the third row is attributable to the differences between mature and juvenile fibres mentioned earlier. The final row shows the variation caused by the earlywood/latewood transition in annual growth ring. The relative density of a fibre gives a good indication of the fibre wall thickness which contributes to the fibre stiffness.  Chapter 1. Introduction  6  C e l l W a l l Chemistry The cell wall derives its strength from the substance cellulose. Cellulose is a carbohydrate made up of repeating glucose units which form extended chains. When these long chain molecules fit together, regions of crystallinity develop. These regions give the cellulose great strength and are difficult to pentrate by solvents. The structure of the cell wall is shown in Fig. 1.2. At a larger scale, aggregates of cellulose molecules form fibrils. The fibrils align together to form the layers of cell walls. The orientation of the fibrils plays a key role in the strength of the fibre. Fibrils become more aligned along the fibre length as a tree becomes more mature. This increases the elastic modulus of the fibre so that more mature fibres have a higher stiffness than less mature fibres. Shorter-chain molecules are also found in the cell wall. These polymers, known as hemicelluloses, are formed from five different repeating sugar units. The matrix which holds contiguous fibres together is a highly polymerized substance known as lignin. The concentration of lignin in the cell walls increases as one moves way from the lumen, the hollow middle section of the fibre. The area between fibres is known as the middle lamella. The objective of chemical pulping is to remove the lignin and release the individual fibres. In practice, a certain amount of cellulose and hemicellulose is also removed. Mechanical pulping attempts to achieve the same result by heating and pulling fibres free from the lignin. This processes leaves more lignin intact in the fibres and removes less cellulose than chemical means. The yield, defined as the percentage of pulp mass after pulping, is less than 50% for some chemical pulping processes. Mechanical pulps may have yields better than 90%. The type of pulping used has many consequences in further processes. In particular, mechanical pulp is much stiffer than chemical pulp. Heartwood refers to the dead cells located in the centre of the tree. Living cells located in the outer region of the tree are known as the sapwood. More organic deposits are found in the heartwood than in the sapwood. As a result, heartwood fibres will be more likely to retain some of their stiffness after pulping because of the insolubility of these molecules. The lignin content in a tree also decreases as the tree matures.  Chapter 1. Introduction  7  Figure 1.2: The structure of the fibre wall. ML refers to the middle lamella, the region where lignin bonds adjacent fibres. P refers to the the primary wall, the thin outer layer of the fibre. Fibrils are randomly oriented as shown. 51, 52, 53 are the three layer of the secondary wall. Each layer has a distinctive fibril alignment. 52 forms the bulk of the fibre wall. W refers to the lumen, the cental canal of the fibre which is empty. Taken from Fig. 2-5 in [19].  Chapter 1. Introduction  1.3  8  Pulp Processing  The raw material for the pulp and paper industry comes predominantly from forest product operations in the form of wood chips. Automated mills process the chips into pulp which may be used for a variety of paper products. Pulp, the fibrous material from the tree, is released through mechanical or chemical means or by a combination of these methods. The processing of pulp occurs at high capacities. As a result, mills require large capital investments and consume massive quantities of energy. No single process may be allowed to bottleneck the operation of the entire mill. Pulp processing consists of several major steps. The raw material, wood, must first be converted into pulp. The pulp is then cleaned to remove any undesirable materials. Finally, the pulp is refined and blended with other pulp stocks before being sent to the papermaking machines. Manufacturers attempt to reduce processing costs by selecting pulp supplies according to the desired physical properties of the paper product. However, there is a wide variance in the natural properties of wood fibres across both growth season and species. Different pulping techniques also result in changes to the pulp.  Refining Several pulp properties are necessary to produce acceptable paper. A minimumfibrelength is required for bonding between fibres. This property is especially critical for tear strength (ie. localized stress in paper). The key to pulp quality is fibre flexibility. Earlywoodfibreswhich have thinner walls and lower density have an advantage in this area. They collapse easily into ribbons during sheet formation. This increases the area available for bonding between adjacent fibres. However, manyfibresdo not,have these properties and must be refined to change their physical properties to those required for the desired end use. This involves a mechanical treatment called beating which presses fibres between a stator and a rotor. This flattens the fibres and creates nicks on the surface. These effects increase the ability of the fibres to bond together to  Chapter 1. Introduction  9  form a strong paper product. Mechanical treatment affects pulp physical properties in different ways. The external surface area is increased as fibrils (cellulose surrounded by hemicellulose) are released from the surface. These serve to increase the bonding surface between neighbouring fibres. The inner wall delaminates, increasing the surface area and reducing the bending strength of the fibre. Many fibres are broken resulting in a change in the length distribution and a lower mean fibre length in treated fibres. Thefibresare also deformed. Changes in geometry occur such as flattening, curling, and kinking. Flattening reduces a fibre's area moment of inertia. This allows it to conform more easily when bonding with neighbouring fibres. Bonding is also aided by the flat surfaces. Curl is a permanent curve set in a fibre. Kinks are local areas of wall damage and usually result in the formation of a hinge where the fibre will be much weaker in one direction. On a smaller scale, crimps, dislocations and microcompressions may be induced or removed from the walls. An excellent collection of photomicrographs which show these changes in structure may be found in [18]. The damage to fibres (bends, kinks) are caused by the refining rotor pressing its edge into a fibre when it is lying on a stator.  Chapter 1. Introduction  10  Fibre Fractionation Obviously, there is a strong desire for bringing the properties of the pulp supply closer to those desired for the finished product. This is presently being achieved in industry using fibre fractionation.  Fibre fractionation is the mechanical separation of fibres from a mixture to  produce at least two fractions which have higher percentages of fibres with certain properties. This process occurs between the pulping and refining stages in the mill. Each fraction may be refined separately before being recombined. This has two advantages. First, less energy is expended in refining the fraction which is closest to the desired properties. This in turn reduces the amount of structural damage to fibres. Fibres which are not damaged during refining are able to retain more of their natural strength. This added strength can then be used to advantage in many products. Reduced refining energy also decreases the amount of choppedfibres(known as fines). Reducing the amount offinesincreases the freeness, or ability of pulp to drain. This allows a higher throughput in the papermaking machine. This consideration is of added importance for recycled fibres. Having been through the manufacturing process at least once, these fibres are more susceptible to breakage. Pulp mills have a wide variety of devices to remove unwanted particles from papermaking fibres. Many of them can also be used to fractionate fibres. Two common types of devices for fractionation are pressure screens and centrifugal cleaners (hydrocyclones). They serve to remove impurities such as dirt or plastic as well as shives (bunches of unpulped fibres) from the slurry. By further understanding the mechanisms which control their operation, refinements may be found which will allowfibresto be separated on the basis of small differences in physical properties. Pressure screens utilize a cylindrical, perforated screen located inside a vessel. The pulp slurry is introduced into the vessel so that it must either pass through the screen (accepts) or remain outside and exit through a discharge (rejects). Pressure screens may operate on either of two mechanisms. In barrier screening, the perforations are smaller than the width of the unwanted particles. The screen serves as the barrier. The  Chapter 1. Introduction  11  INLET SIDE  ACCEPTSSIDE Figure 1.3: Schematic diagram of pressure screen operation. The rotor moves past the screen creating a pressure pulse. Taken from Fig. 9-25 in [19] disadvantage of this method is the low throughput of slurry caused by the small perforations required to remove all the undesirable particles. Probability screening is more widely used in industry. In this method the perforations are shaped as slots in the screen. Their dimensions are larger than the smallest diameter of the desired particles. In the case of long particles, such as fibres, it is the length and orientation of the particle as it approaches the slot which determines whether it will pass through. The ability of a particle to conform to a streamline as it approaches the slot will also determine its chance of passing through. In this way, probability screening fractionates on the basis of particle length and flexibility. Screening is complicated by the fibre mats which form over the screen perforations. These are removed by a rotating foil (Fig. 1.3), which causes a pressure pulse. The reversed flow through the slots purges the openings of fibres.  Chapter 1. Introduction  12  The influence of the fluid flow on the orientation and deformation of the fibre near the the slot opening is clear. An understanding of this area of screening is vital for further improvements in fractionation technology. Hydrocyclones are another type of equipment which could be adapted to fibre fractionation. Fig. 1.4 is a schematic diagram showing the general principles of operation of a hydrocyclone. A slurry is introduced tangentially to the hydrocyclone wall at the top of the inverted conical vessel. As the slurry travels in a spiral down the vessel, the centrifugal forces cause the denser particles to migrate to the walls of the vessel. The denser particles then become entrained in a downward flow which takes them to the reject outlet at the bottom of the vessel. The less dense pulp fibres remain in the swirling flow, are carried inward and upward where they are removed through an outlet at the centre of the top of the hydrocyclone. Hydrocyclones use a combination of centrifugal force and fluid shear to separate particles on the basis of particle density and size. Centrifugal forces separate fibres of differing densities as described above. Shear forces are required to pull apart the fibre networks which are known as floes. This thesis is a first step in modelling the motion of a flexible or rigid fibre in a shear flow, and is aimed at exploring the effects of fibre properties on fibre motion in complex fluid fields such as those found in a hydrocyclone.  Chapter 1. Introduction  Figure 1.4: Schematic diagram of a hydrocyclone. Taken from Fig. 9-30 in [19]  13  Chapter 1. Introduction  1.4  14  Literature Review  The dynamics of particles suspended in a fluid are of fundamental importance to many areas of industry and technology. As such, a large body of literature exists on this topic. This section will review the previous work conducted on the motion of rigid and flexible individualfibresin a fluid velocity field. Particular emphasis will be placed on those aspects which relate to pulp fibres in dilute suspensions.  1.4.1  Rigid Fibre Motion  Jeffery derived the equations of motion for an isolated rigid neutrally buoyant ellipsoidal particle in Stokes flow in a Newtonian fluid in 1922 [12]. According to his theory, the centroid of the particle assumes the translational velocity of the fluid it displaces. A Couette flow (ie. two dimensional flow field with a linear profile) causes the particle to rotate about the centroid. The formulation for the period of rotation is referred to as Jeffery's equation in the literature.  TG = 2Tr(r + r: )  (1.1)  1  e  where T is the period of rotation, G is the shear rate, and r is the ratio of the length of the e  ellipsoid to its maximum diameter. For large values of the ellipsoidal axis ratio, r , Jeffery's e  equation simplifies to a linear fuction,  TG  « 2irr  e  (1.2)  The rate of rotation was formulated as,  where <p is the orientation angle of the ellipsoid as it rotates in a clockwise direction from the y axis as shown in Fig. 1.5.  Chapter 1. Introduction  15  Figure 1.5: Orientation angle of Rotating Particles When integrated with time, an equation for the angular displacement is formed,  <t> = i a n [ r tan(Gt _1  e  T  e  * +1  e  )]  (1.4)  The equations of motion of the particle are approximated by neglecting the terms involving the squares of the fluid velocity. This assumption of Stokes flow means that the fluid viscous forces predominate over inertia! forces. The condition for the validity of this assumption is that the Reynolds number is less than one,  Re= -^<l P  (1.5)  that is, the product of the relative velocity of the ellipsoid by its length, ul must be smaller than the kinematic viscosity of the fluid, v = r\\p.  Chapter 1. Introduction  16  70  60  50 * r  e  40  30  20  10  , , i , , 20  i  i  i • , , , 1 , , ,' , 1 . . . . 1 . . 40  60  80  100  120  Cylindrical Aspect Ratio, r  c  Figure 1.6: Equivalent Ellipsoidal Aspect Ratio, r* versus Cylindrical Aspect Ratio, r . This figure was adapted from data in [22], [9], and [14]. c  Bretherton [3] subsequently extended JefFery's work to describe the motion of any axisymmetric particle by using the concept of an equivalent ellipsoidal axis ratio. In 1950, Mason [13] stated that pulp fibres may be regarded as cylindrical particles. Trevelyan and Mason [21] showed that JefFery's equation could describe the motion of a cylindrical particle by substituting an equivalent ellipsoidal axis ratio, r* for r in JefFery's equation. The value of e  r* was calculated by comparing measured periods of rotation from experiments with cylindrical rods of aspect r with the periods obtained from JefFery's equation. Fig. 1.6 exhibits the c  variation of r with r*. The difference in shape between the rods used in the experiments and c  the cigar shaped ellipsoids upon which the theory is based results in differing rates of rotation.  Chapter 1. Introduction  17  The assumption of Stokes flow for the motion of a suspended pulp fibre allows Jeffery's equation to be used to describe its motion. This assumption is supported by the small size of a typical softwood fibre (length = 4 mm). It can also be assumed that the fibre will move at a low velocity relative to the surrounding fluid due to the similarity in densities. The fibre wall has a specific gravity of approximately 1.5. Depending on the thickness of the fibre wall the fibre will have an overall specific gravity « 1.1 when immersed in water. Jeffery's equation assumes that, apart from the disturbance produced in the immediate vicinity of the particle, the fluid motion is steady and varies in space on a scale which is large compared with the dimensions of the particle. The assumption of an isolated fibre in a Newtonian fluid also serves as the basis for the computational model developed in this thesis. For this reason, Jeffery's equation, using r*, will serve as a benchmark for the rigidfibremodel.  1.4.2  Flexible Fibre M o t i o n  The modelling offlexiblefibremotion in a suspension is much more difficult than that of rigid fibre motion. Even a theoretical model of an isolated fibre will be limited by the fact that the fibre shape, and hence forces acting on it, is not constant.  Much of the literature has  focussed on experiments conducted with pulp and other types offibresin order to describe the characteristics offlexiblefibremotion. A large body of work on the motion of particles in suspension was created by Mason and others after 1950 at the Pulp and Paper Research Institute of Canada (PAPRICAN). Much of it was published as a series of articles in the Journal of Colloid Science [21, 1, 6, 7]. This work has particular application to the present research due to its focus on the theory and classification of pulp fibre motion.  Chapter 1. Introduction  18  Mason [14] showed that the motion of a single wood fibre is more complex than the model in Jeffery's equation due to the fibre flexibility. Forgacs, Robertson, and Mason [8] divided the range of motion of fibres into the four "orbit classes" shown in Fig. 1.7. The fibre orbits are described below: • Type I - Rigid Rotation The fibre remains rigid throughout the rotation and its period and rate can be described by Jeffery's equations. • Type II - Springy Rotation The fibre begins unstressed in a horizontal position. As it rotates, a compressive force acts on the fibre due to the fluid. This force reaches a maximum at 45° and the fibre responds by bending. The deformed rotation continues until the fibre approaches 135°. At this point the fibre straightens and continues its rotation. • Type III - Flexible Rotation As the flexibilty of thefibreincreases, one or both ends of the fibre bend inward soon after it has rotated from the horizontal position. In a loop turn both ends bend in to create an S-shape. This shape persists in the fibre as it rotates before straightening along the x axis near 180°. A snake turn occurs when the forward end of the fibre bends upward and undulates down the fibre until it is facing backward. This turn is more common than the loop turn in actual shear flows due to asymmetries in the fibre. • Type IV - Complex Orbit When the fibre is sufficiently flexible it will continue to rotate in a snake turn motion without straightening at any point in the rotation.  Chapter 1. Introduction  t/ 1  —  19  I —  \  /  r -  x  II  -  v (  —  s  —  III  —  IV  < _  O  C L _  s  J  c  c  -  —  -  —  c  Figure 1.7: Typical rotational orbits of fibres in a one dimensional shear field as shown on the top left of the figure. The fibre flexibility increases from Classes I to IV. Adapted from Fig. 14 in [8]. The authors [8] also identified factors in actual pulp fibres which account for their behaviour. In particular, the nonuniformity offibrestiffness was noted. There are sections along the length of the fibre which are more flexible than others. There often exist localized areas of weakness which are caused by fractures in the cell wall.  Chapter 1. Introduction  20  Onset of Buckling The ability of a fibre to resist deforming in a shear flow depends on its stiffness and the flow parameters: viscosity and shear rate. Forgacs and Mason [6] put forward a theory for the deformation of cylindrical particles rotating in a velocity field. Mason had stated in an earlier paper [13] that wood fibres could be regarded as cylindrical in shape. Equations were developed which calculated the critical value of the ratio E/Grj above which the shear-induced axial compression causes the particle to buckle. The theory was based on an approximate theory of the forces acting on a thin rod in laminar shear developed by Burgers [4]. For a rigid particle rotating in a one dimensional shear, Burgers showed that the total axial force, F acting on the central cross section of the fibre is approximately,  F  =  , ,  f  TG  ln(2r ) - 1.75  M  (1.6) V  ;  c  (1.7) where I is half the fibre length, r is the cylindrical aspect ratio (length to diameter ratio), c  M = sine)) cos<j>, and </» is the angle of the rod measured clockwise from the y axis. The sign of F depends upon whether M is negative or positive, and determines whether the fibre is under compression or tension. The maximum compressive force is reached at d> = —45°. By considering Euler's classical equation for the shape of a solid rod suffering small deformations under compressive forces, the critical value for the first mode of buckling was found,  E W " *  2rt " l<2r ) - 1.75  (  L  8  j  c  where E is the Young's modulus of the rod. Experimental results showed reasonable agreement with the theory. This equation identifies two dimensionless parameters which determine the bending of a solid cylindrical particle: r and E/Grj, the ratio of the fibre stiffness to the fluid shear. c  is a function of the fibre geometry (T" ) and is independent of the size of the particle. c  (E/Gr])^  Chapter 1. Introduction  21  Change in T G due to Flexibility and Permanent Bend In 1959, Forgacs and Mason [7] showed that for a particle undergoing shear-induced deformation, the flow pattern around it will depend at any instant on the deforming forces. In particular, they found that the product T G decreased, for a given fibre geometry and flow condition, as the flexibilty was increased. They also found that the period showed a high sensitivity to fibre curl. They reported that a permanent 20° curl caused a 9 fold decrease in the period due to the increased velocity gradient over the fibre.  1.4.3  Computer Simulation of Particle Dynamics  Forgacs and Mason [7] noted that there was no adequate theory to describe the motion of flexible fibres in suspensions. In recent years, computer simulations have been carried out which allow this limitation to be overcome. Brady and Bossis [2] proposed a Stokesian dynamics method for the motion of spheres that accurately dealt with the hydrodynamic interaction between particles. Yamamoto and Matsuoka [22, 23, 24] have proposed a method of simulating the motion of arbitrarily shaped particles accompanied by deformation by modeling the particle as a string of spheres. These authors have applied their numerical method to the rheology of injection molding offibrereinforced resins. Initial results from their model utilize a free draining approxiation that ignores the hydrodynamic interaction between the flow fields of neighbouring spheres. Their results for rigid cylindrical rods follow Jeffery's equation closely within 8% for a 10 sphere model. Results for flexible rods [22, 23] follow that of Forgacs and Mason [7] qualitatively. The approach taken in the present research is based on the work of Yamamoto and Matsuoka. Their approach has the advantage of allowing localized changes infibrestucture to be modelled so that the motion of the fibre may be determined with a variety of flexibilities, both from one fibre to another and within a single fibre. The development of this technique is found in the next chapter.  Chapter 1. Introduction  1.5  22  Scope and A i m s of This Thesis  It is clear that there is a recurring desire to further understand the motion of fibres in shear flows. This thesis will develop a numerical model of fibre motion in a shear flow. A series of numerical tests will be conducted which identify how physical properties change the fibre movement. The characteristics of a fibre which will be studied in particular are: • Aspect Ratio (length to diameter ratio) • Bending Stiffness • Kinks (areas of localized wall weakness) At this point, the model will be restricted to a Couette flow. This flow will be fully decribed by two parameters; G the shear rate and 77 the viscosity. In the future, the model could be extended to three dimensional flows. The results from the trials will be compared by developing a dimensionless number. This number will relate the behaviour of fibres of different aspect ratios and stiffnesses. In this way, a general parameter describing the model behaviour will be available.  Chapter 2  Numerical Procedure  This chapter outlines the development and implementation of a numerical model offibremotion in a shear flow. The defining equations of the model are developed in the Model Formulation section. First, a model defining the inertial and elastic properties of the fibre is described. This model is based on a fibre constructed of cylindrical elements. Next, the flow field is defined and the interaction between the fibre and the flow field is considered. The hydrodynamic force exerted on each element is calculated for a sphere of the same diameter as the cylindrical elements. The apparent contradiction between these two conceptions of the element geometry is simply a matter of convenience in the present calculations; other models of the elastic, inertial, and hydrodynamic properties could be substituted. The equations of motion of the fibre are developed in the Method section. The conditions required to maintain continuity between elements of thefibreare presented. The finite difference technique and algorithms used to follow the fibre through time are then discussed in light of the continuity condition. The Application section describes how these algorithms are applied in a Fortran computer code. The dimensionless equations used in the code are introduced and the parameters which control the model's performance are discussed.  23  Chapter 2.  2.1  24  Numerical Procedure  Model Formulation  A m o d e l has been developed to describe a fibre and its d y n a m i c behaviour i n a t w o dimensional shearing flow field. In order to provide reasonable values for its mechanical properties, a single fibre is considered to be a cylindrical r o d . T h e fibre is discretized i n t o a series of elements of equal length a n d diameter. E a c h element can stretch a n d bend w i t h respect t o its neighbours by changing its b o n d distance a n d angle. T h e continuity of a real fibre is preserved i n the m o d e l by m a i n t a i n i n g a contact point at the interface of elements. T h e m o t i o n of the fibre as a whole is determined by solving the translational and rotational equations of i n d i v i d u a l elements under the h y d r o d y n a m i c force a n d torque exerted by the relative m o t i o n of the fluid. T h e fibre is tracked for a specified time period by a L a g r a n g i a n scheme, the m o t i o n over a longer time being composed of a series of steps for consecutive time increments.  2.1.1  Fibre Structural Model  T o represent its structural properties, the fibre is made up of N cylindrical elements of outer radius a  0  a n d length 2a .  diameter 2a  D  Q  In this manner, the model forms a cylindrical r o d of length 2a N, Q  a n d aspect ratio N. T h e m o d e l represents the hollowness of pulp fibres t h r o u g h  the use of a n inner radius, a{. T h e thickness of the fibre walls are set by the ratio of ai a n d a . a  T h e walls are assumed to have a constant density. T h e interior space of the fibre is assumed to be filled w i t h water. T h i s p o r t i o n of the element adds mass to the m o d e l based o n the density of water b u t has no effect o n the structural properties.  Chapter 2.  25  Numerical Procedure  F i g u r e 2.8: Schematic d i a g r a m of a continuous fibre and its representation by a series of c y l i n d r i c a l elements T h e elastic properties of the fibre m o d e l are represented i n the m o d e l t h r o u g h the use of spring a n d bending constants. T h e product of these constants and the b o n d distance and angle determine the m a g n i t u d e of the spring force and bending torque respectively. T h i s is better understood i f one considers a pair of connected elements i and j.  T h e spring force exerted  between elements depends on the b o n d distance I between the centres of the elements as shown i n F i g . 2.9. If the elements are stretched from their e q u i l i b r i u m distance (taken as l  Q  an equal a n d opposite force F  s  = 2a ), D  is exerted between the elements d r a w i n g t h e m back to the  e q u i l i b r i u m distance:  F = s  where k is the spring force constant. s  -k  s  (I  - y  (2.9)  26  Chapter 2. Numerical Procedure  I  J  F i g u r e 2.9: C o n f i g u r a t i o n of bonded elements when stretched. T h e spring force F function of I, the distance between the centres of the elements.  s  T h i s force acts along the n o r m a l  is a linear  , a unit vector p o i n t i n g from the centre of element i to  the centre of element j ,  (2.10) T h e magnitude of the spring constant is based on the linear deformation of a c y l i n d r i c a l r o d . It is evaluated by using the Y o u n g ' s modulus E and the cross sectional area of the m o d e l . T h e full development of the spring force equation and constant is found i n A p p e n d i x B , where it is shown t h a t  k. = *-Ea (1 - g )  (2.11)  0  where  is the radius of the hollow inner p o r t i o n of the fibre and a is the outer radius. a  Chapter 2. Numerical Procedure  27  F i g u r e 2.10: C o n f i g u r a t i o n of bonded elements when bent. T h e bending torque T  b  is a linear  f u n c t i o n of Ob, the b o n d angle formed between the elements. A s i m i l a r approach is taken i n c a l c u l a t i n g the bending between elements.  T h e bending  torque is a f u n c t i o n of the b o n d angle Ob between the elements. A b o n d angle is f o r m e d when neighbouring elements rotate relative to each other away f r o m the e q u i l i b r i u m angle Obo (taken as zero i n the case of straight fibres). A n equal and opposite torque T  b  is then exerted on each  element,  T  b  = -h(0  b  - Obo)  (2.12)  where fc& is the bending torque constant. T h i s torque acts about a vector which is perpendicular to  . T h e m a g n i t u d e of the bending  constant is based on the bending deformation of a c y l i n d r i c a l r o d . It is evaluated by using the  Chapter 2. Numerical Procedure  28  Y o u n g ' s m o d u l u s E and the second moment of area of the m o d e l about about the n e u t r a l axis. T h e f u l l development of the bending torque equation and constant is f o u n d i n A p p e n d i x B , where it is shown t h a t :  k = \Eal 8 b  (1 a%  (2.13)  A m o d e l constructed i n this manner allows the flexibility of the fibre to be changed t h r o u g h the b o n d parameters k  s  and k\,. T h e overall flexibility m a y be adjusted b y v a r y i n g the stiffness,  E, or the inner a n d outer r a d i i , ai and a . a  L o c a l i z e d weaknesses such as kinks m a y be modelled  by adjusting the b o n d parameters at i n d i v i d u a l locations along the fibre length.  Chapter 2. Numerical Procedure  2.1.2  29  Hydrodynamic Force Model  The formulation of a model for the drag and torque exerted on each element by the fluid flow is considered next. The fibre model is assumed to be immersed in a Newtonian fluid with an absolute viscosity n. Theflowfieldis modelled as a simple two dimensional shear flow as shown in Fig. 2.11. The fluid has a velocity in the x direction which changes linearly in the y direction,  U(y) = Gy  (2.14)  where G is the shear rate and y is the distance above the origin. The flow field is steady with time. Theflowfieldhas no velocity components in either the y or z direction and does not change in the x direction,  U(x)  = constant  (2.15)  V(r)  =  0  (2.16)  = 0  (2.17)  W(r) where r is the global coordinate.  The macroscopic rotational velocity fi of the fluid is based on the shear-induced vorticity (defined as positive in the counter-clockwise direction),  (2.18) (2.19)  Chapter 2. Numerical Procedure  30  r i  '  U = Gy  i  t=  t>0  0 •  / F i g u r e 2.11:  T w o dimensional shear flow. T h e fluid velocity i n the x direction, U is a linear  f u n c t i o n of the shear rate G and the y position. W h e n a fibre is suspended i n a shear flow, a h y d r o d y n a m i c force a n d torque are exerted on it by the fluid. T h i s phenomena is modelled through the choice of suitable elements. T h e h y d r o d y n a m i c m o d e l involves several i m p o r t a n t assumptions and simplifications. A t the largest scale, the fibre is considered to exist under the condition of an infinitely dilute system. In other words, a single fibre is being modelled i n a flow field unaffected by walls, or other  fibres.  A  "free d r a i n i n g " a p p r o x i m a t i o n is applied to the h y d r o d y n a m i c force acting on each i n d i v i d u a l element. T h i s means that the force acting on an element is considered t o be solely a f u n c t i o n of the bulk fluid m o t i o n at the location of the centre of the element. Changes i n the flow field induced by the m o t i o n of neighbouring elements are neglected.  Chapter 2. Numerical Procedure  31  A low Reynolds number flow regime over the element is assumed. The rationale for this assumption with respect to entire fibres was developed in the previous chapter. The forces and torques exerted by the fluid on the elements are based on Stokes flow. The inertial forces exerted by the fluid are assumed to be negligible; in other words, forces exerted on the elements are due to shear forces alone. The hydrodynamic force F  h  and torque T  h  are calculated according the  magnitude and direction of the relative flow at the centre of the element as shown in Fig. 2.13. The magnitude of the forces is based on Stokes formulations for spheres of radius a . The Q  translational friction force is then  F = -( [u-U(y)} h  t  where ( = t  67r77a0  (2.20)  is the translational friction constant. Similarly, the angular friction torque is  T  h  = -C [u, - n(r)] r  (2.21)  where £ = 8irrjag is the rotational friction constant r  There may appear to be a contradiction in using a cylinder as the basis for choosing the mass and elastic constants of the model and a sphere for choosing the hydrodynamic constants. This has been done to make the calculation of the fluid forces as simple as possible. In this model, the fluid force always acts in the x direction (ie. only a drag force is assumed). If a cylinder had been used for the calculation of the hydrodynamic force, the force acting along the axis of the element would have to have been considered separately from the force acting normal to the element. The relative magnitude of these forces would be related to the orientation of the element to the direction of the flow. This would also be true for the torque exerted on the element by the fluid. The use of a sphere to model the fluid drag and torque was adopted for simplicity.  32  Chapter 2. Numerical Procedure  Hydrodynamic Model  Structural Model  Figure 2.12: Element geometries: spherical and cylindrical. The outer radius a defines the geometry in both the hydrodynamic and structural models. The inner radius defines the fibre wall thickness in the structural model. 0  A model forfibremotion in a two dimensional shear flow has been developed. The structural properties of the fibre are based on those of a hollow or solid cylindrical rod. Individual elements have a diameter and length 2a . The hydrodynamic modelling is based on free draining spheres a  of radius a . These geometries are shown in Fig. 2.12. The hydrodynamic forces are linearly Q  related to the relative velocity between the fluid and element. This allows them to be easily broken into directional components.  At this stage of the model development, only physical  (spring or bending) and hydrodynamic forces are considered in the model. Body forces such as gravity have not been included.  Chapter 2. Numerical Procedure  2.2  33  Method  The model simulates fibre motion by solving the governing equations of motion for each element over a series of time steps. The terms of the governing equations are derived from the two models proposed in the previous section. The results are produced as data sets which list the location of the centre of each element in the model. A plot of these data points serves as a "snapshot" of the fibre's centreline at a particular time step. However, at this stage of development, the model does not reflect the continuity of an actual fibre. The structural force, F ensures continuity in the axial direction. However, no analagous 8  force exists for the tangential direction. This is resolved by adding a tangential force /£• between elements which acts perpendicular to THJ. This force is exerted on element i by element j at their contact point as shown in Fig. 2.13. It acts in an equal and opposite manner on element j. Adding this term, the general equations of motion for translation and rotation are then formed for each element,  (2.22)  (2.23) where m and I are the mass and mass moment of inertia of each element respectively. The e  e  formulae defining these inertial terms are developed in Appendix A. The summation terms in the equations imply consideration of neighbouring elements on both sides of element i but do not extend to more than the immediate neighbours in any case. The magnitude of /* is found by considering a non-slip condition at the point of contact between each element. This condition requires that the translational velocities of neighbouring elements be the same at their contact point. In terms of the generalized cordinates T{ and 6{, this condition is written as,  34  Chapter 2. Numerical Procedure  F i g u r e 2.13: H y d r o d y n a m i c forces exerted on the elements cause t h e m to move apart as seen i n the two outer elements. Tangential forces act at the interface to prevent this.  dri dt  dOi dt  J  drj dt  d6j  T h e contact point is assumed to occur at a distance a  0  (2.24) f r o m the centres of b o t h elements  along n-ij. T h e non-slip condition is applied so that the velocities perpendicular to n{j are equal at the contact point. T h e magnitude of / * is not directly calculated i n this procedure. Instead, a set of simultaneous equations for u*, the velocity due to jf* must be solved to ensure an equal velocity at each contact point. These equations are formed by first solving E q s . 2.22 and 2.23 to  find  the accelerations while neglecting the f* terms. New velocities at t + 6t are found by adding the contributions of the accelerations using a finite difference technique. A new set of velocity equations including the u n k n o w n velocity u  l  due to jf* are formed. It is these new velocities  which are substituted into E q . 2.24. T h e terms involving v} are then isolated so that they may be solved. F o r a fibre made up of N elements, there w i l l be (N — 1) interfaces and (N — 1) simultaneous equations which must be solved to determine the velocity contribution of jf*.  Chapter  2.  Numerical  e + Qj = 2a i  35  Procedure  ij)  tf = 0  Q + Qj<2a j, i  iJ!  h <o g  F i g u r e 2.14: Schematic D i a g r a m showing nonslip condition. Oi + Oj = 2ctij  T h e use of finite difference solutions to differential equations introduces n u m e r i c a l errors w h i c h cause the m o d e l to violate precise continuity. In particular, an angular non-slip condtion may be broken when the m o d e l bends. T h e geometric form of the condition is:  Bi + Oj = 2a  i3  where Q ^ J is the angle of  (2.25)  from the fixed axes.  T h i s condition requires that nieghbouring elements " r o l l " w i t h respect to each other so that their contact point does not slip. A d i a g r a m of this condition can be seen i n F i g . 2.14. T h e use of spheres i n the figure implies that the contact point remains at a constant distance a  D  from the centres of the elements. If longer elements were used i n the m o d e l the condition would change to reflect the new geometry. A n o t h e r possibility i n the future development of the m o d e l w o u l d be to use a coordinate system which guarantees continuity of the fibre automatically.  Chapter 2. Numerical Procedure  36  In the present model the correction is applied by adding the contribution of a torque,  h[(0j + Oi) - 2aij]/C  rot  (2.26)  where k is a constant that is adjusted according to the degree of correction desired. In practice g  the value of the bending constant kf, was also used for k . The rotational acceleration induced g  by t results in an additional contribution to the angular velocity of each element over the time g  step 6t. Average velocities (u, v, and u) over the time step St are calculated from the mean of the velocities at the start and end of the time step. New translational and rotational displacements are found by adding the contribution of the velocities over the time step. These values are recorded for each element and the time step is advanced. New equations of motion are then formed by calculating the hydrodynamic and fibre structural forces. The motion of the whole fibre may be followed by continuing this process over a series of time steps. The computational algorithm followed in the application of the procedure is detailed in the next section.  Chapter 2. Numerical Procedure  2.3 2.3.1  37  Application Non-dimensional Groups  In order to implement the numerical procedure outlined in this chapter, a computer code was written in the Fortran programming language.  First the equations were expressed in  dimensionless form. It is thefluidforces which are responsible for the motion of the fibre in this model. As such, a r\G was chosen as the unit of force in the simulation. The units of length and 2  0  time were selected as a and 1/G respectively. Other dimensional units used in the simulations a  are identified in Appendix C as well as the dimensionless equations of motion. These equations are as follows,  %  =  ^  ^;(^  =  +E ^  ^(E^  +£ n  (2-28)  + E n  % = ^( ' * + E * + E / r l  (2-27)  r6  7  t  *X^ +  E *) r9  (2-29)  where * refers to the dimensionless value. C . and C t are the translational and rotational t7  TO  constants which relate the inertial and hydrodynamic forces in the model. Their formulation is shown below.  C t T  =  Crot  27r^^  =  ^  (2.30)  (2.31)  Chapter 2. Numerical Procedure  C  tr  38  is proportional to Re , the particle Reynolds number defined as Re = p  p  a  °y . G  This  dimensionless group is unlike the traditional Reynolds number since it is the inertia of the element that is considered in its formulation; not that of the fluid surrounding it. Fluid inertia! forces are neglected because of the assumption of Stokes flow. The element density, p includes e  the contribution of the fibre wall and the fluid immobilized inside the walls. The fluid density thus has an effect on the inertia of the element. The term particle Reynolds number will be used in this thesis to emphasize the use of the particle's inertia as opposed to that of the fluid. C t is related to Re in a more complicated manner than Ct because of the form of the ro  p  r  element's mass moment of inertia I which is, e  U =  (3a + lal) + ^ 2  t  where rrif is the mass of the fibre wall and m w  w  (3a + 4a*) 2  t  (2.32)  is the the mass of the fluid immobilized inside  the hollow fibre. As can be seen from the equations for J , the inner radius a{ of an element has an effect e  on its inertial properties. It does so by changing the wall thickness of the fibre as well as the amount of immobilized fluid. The wall thickness is represented nondimensionally in the model through the ratio of inner to outer radius in the fibre,  a* — ai/a  a  (2.33)  The elastic properties of the fibre are defined by the product EI, where E and I refer to the Young's Modulus and second moment of area respectively. Changing the wall thickness ratio, a*, affects J as shown in Appendix B. A list of the parameters and the forces which they influence is given in Table 2.2.  Chapter 2. Numerical Procedure  39  Force Inertia!  Model Fibre  Hydrodynamic  Fluid Fibre Fluid  Elastic  Fibre  Parameter Outer Radius Inner Radius Fibre Wall Density Fluid Density Outer Radius Shear Rate Viscosity Outer Radius Inner Radius Elastic Modulus  Table 2.2: The fibre and fluid parameters which define the model. Note that the fluid density which effects the fibre inertia refers to the fluid immobilized inside the fibre walls. Since the fibres are immersed in water it is assumed that they fill with water.  A dimensionless group involving E appears in the process of forming the equations. The elastic coefficient,  k  el  = E/Gri  (2.34)  relates the elastic force to the hydrodynamic force. This group is found in the dimensionless forms of the bending and stretching constants, A£ and A:*. These constants determine the rigidity of the fibre. They are functions of the Young's modulus, and the cross-sectional properties (ie. area, second moment of area). The forms of these constants are:  *:  =  EA /a r,G  (2.35)  =  f(k i,a* )  (2.36)  k* =  EI /atvG  (2.37)  =  /(W )  (2.38)  b  2  C3  0  2  e  ca  4  Chapter 2. Numerical Procedure  40  Group  Meaning  Re  Inertial / Hydrodynamic  a*  Elastic / Hydrodynamic Aspect Ratio Wall Thickness Ratio  p  Definition V  E/r/G r =N a* — a.i/a c  0  Table 2.3: Dimensionless groups and the parameters involved in their formulation  The dynamic behaviour of a fibre is influenced strongly by its shape. Higher aspect ratio fibres rotate at a slower rate in a shear flow than lower aspect ratio fibres. This geometric factor is represented by r , the cylindrical aspect ratio which is equal to the length of the fibre c  divided by its width. Since the elements used in the model have equal length and dimameter the aspect ratio will simply be equal to N, the number of elements.  Chapter 2. Numerical Procedure  2.3.2  41  Computational Algorithm  This section outlines the algorithm that was followed to perform simulations. It closely follows the code listing given in Appendix F. 1. PreUminary Declarations Before beginning calculations, all variables other than counters are declared as double precision. The number of elements, N is set as a parameter. Array declarations are then based on this parameter. The number of iterations, along with the dimensionless time step, are chosen to simulate the desired total time for the fibre motion. The largest time step allowable is chosen in order to reduce computational time. The issues of calculation time and time step limitations will be discussed in the Results & Discussion chapter. 2. Parameters The parameters which define the fibre (wall density, inner and outer radius, and Young's modulus) and the fluid model (water density, shear rate and viscosity) models are set according to the trial being performed. These values are then used to calculate the crosssectional area, second moment of area, mass, mass moment of inertia, density, bending and spring constant, and the translational and rotational constants. 3. Initial Conditions Before entering the computational loop, the initial conditions of the fibre model are set. In all cases, the fibre is positioned initially with its centreline aligned with the y axis. The importance of the initial position will be discussed in the next chapter. The vertical placement of the fibre is defined by the location of the centre of the first element. Remaining elements are connected above the first. Internal stresses are set to zero as the fibre is straight and therefore undeformed. The velocity and rotation of the elements are initialized to zero.  Chapter 2. Numerical Procedure  42  4. Data Files Before calculations begin, the data files are created and opened. The primary file stores the position of each element at selected time intervals. Other files may be opened to store data on the forces, torques, et cetera acting on the fibre. If the model is continuing the results from a previous trial, a file will be read which contains the data from the last time step of the previous trial. 5. Computational Loop (a) At t = 0, the shear flow begins to act on the motionless fibre and the computational loop begins. A series of calculations are performed which depend on values from the previous time step (or initial conditions, at t = 0). The normal directions (Eq. 2.10), spring forces (Eq. 2.9), and bending torques (Eq. 2.12) are determined from the position of the elements relative to each other. The hydrodynamic force and torque are calculated by Eqs. 2.20 and 2.21 from the relative velocity between the elements and the fluid. (b) Equations of motion (Eqs. 2.22 and 2.23) are then formed. The accelerations are due to the forces and torques calculated above. The velocities from the previous time step are saved. New velocities are found by adding the contribution of the accelerations over the time step,  u(i) = u (i) + u(i)St 0  (2.39)  where u (i) is the velocity of the element from the previous time step. 0  (c) The velocity v} due to the tangential force f is found by solving (N — 1) simultaneous t  equations. These equations are constructed by finding ti* from Eqs. 2.22 and 2.23 and substituting it into Eq. 2.24. A matrix is then formed from the coefficients of u*. The development of this matrix is found in Appendix D. The tridiagonal structure of the matrix decreases the computations required in the inversion subroutine tridiagonal.  Chapter 2. Numerical Procedure  43  The contributions of /* to the translational and angular velocity of each element are then added. A further adjustments are made to the angular velocity so that the model does not break the geometric condition in Eq. 2.25. A torque, t is calculated by applying 9  Eq. 2.26 to the angular values from the beginning of the time step. An additional rotational velocity, u  9  due to the torque is found. The angular velocities are then  updated to include this contribution. When the velocities have been found at the new time step, the displacement of each element over the time step may be calculated. This is done by considering the average velocity over the time step. For the x direction this calculated as,  x(i) = x(i) + 0.5 [u (i) + u(i)] St Q  (2.40)  where u (i) is the velocity at the beginning of the time step and u(i) is the velocity 0  at the end of the time step. The new position of each element may then be recorded in a data file. 6. Continue Loop or End At this point, the program returns to the beginning of the loop or finishes. The time step is advanced, and a new series of calculations and matrix manipulations are begun. Once the loop is finished, the values of all variables used in the program are stored in a final data file. This file may be used if the simulation is to be continued. By repeating the computational loop over an adequate number of time steps, the motion of the elements and hence the fibre may be followed.  Chapter 3  Results and Discussion  3.1  Introduction  A series of trials have been conducted to simulate rotational and translational motion of a single fibre using the numerical model developed in the previous chapter. A preliminary series of tests was conducted with a rigid 10 element model to compare the present calculations with previously published results. These trials were also run to determine the behaviour of the model over a range of values for the inertial constants Ct and C t- Subsequent trials were performed T  ro  with a 60 element model. The second series of trials was conducted over a range of values for the elastic constant k i. Results were determined for fibres with and without bending deformation. e  A final series was completed where the fibre model had localized weaknesses (kinks) in the bending stiffness along its length. The results of these trials as well particular aspects of the model's behaviour will be discussed in this chapter. Motion is induced in afibreby placing it in a Couette flow where it begins at rest with zero velocity as described earlier. After time zero, the shear flow exerts a force on the fibre which responds by rotating in a clockwise direction. The subsequent motion of the fibre is determined by its initial position. A preliminary test was conducted with the fibre centred on the origin along the y-axis. However the remaining tests were run with the fibre aligned with the y-axis extending upward from the origin. This resulted in both rotational and translational motion. A schematic diagram of the two initial positions is found in Figure 3.15. The recurring feature in all trials is the presence of a shearing flow which causes the fibre model to rotate. It follows that the fundamental determinant of fibre motion is G, the shear rate of the flow field and that the dimensionless time unit is t* = tG.  44  Chapter 3. Results and Discussion  45  F i g u r e 3.15: T h e i n i t i a l position of fibres used i n the trials. T h e fibre i n a) w i l l rotate about the origin when a shear is exerted on i t . T h e fibre i n b ) extending up from the origin, w i l l translate and rotate. T h e response of the fibre to the shear can be measured by T , the period of r o t a t i o n . F o r this reason, m a n y of the results discussed i n this chapter will be related to TG, the dimensionless time required to perform one period of r o t a t i o n . T h e c o m p u t a t i o n a l time for the trials was determined by the number of time steps required to calculate the p e r i o d .  M o s t trials were  executed for the time required to complete 360° of r o t a t i o n from the starting position.  Constants and  Parameters  T h e behaviour of the fibre m o d e l is determined by the forces and torques exerted on i t . In the previous chapter, it was shown that three forces are involved: i n e r t i a l , h y d r o d y n a m i c , and elastic. T h e influence of the these forces is determined by the value of the the dimensionless groups shown i n Table 2.3 i n the previous chapter. These groups i n t u r n take their values from the parameters listed i n Table 2.2.  Chapter 3. Results and Discussion  46  Parameter Fluid Density, pf Viscosity, 77 Shear Rate, G Outer Radius, a Fibre Wall Density, pf Young's Modulus, E Aspect Ratio, N  Values Used in trials  1000 kg/m 1Q~ kg/ms 10, 1000 s' 3  3  1  10,100 /xm  Q  1 0 - 10 kg/m 100 - 2 X 10 Pa 4  w  3  6  10, 60  Table 3.4: The constants and values chosen for the parameters. The fibre wall density, pf kept constant at IbWkg/m in the trials with the 60 element model.  w  was  3  The values chosen for the parameters depended on the fibre and flow properties being modelled. Several parameters were kept constant for most trials. The fluid properties were those of water at 20° Celsius so that the fluid density, pf and viscosity, 77 were set at 1000kg/m  3  and 10~ kg/ms respectively for all trials. Again, it should be emphasized that the fluid density 3  plays no role in the hydrodynamic force exerted on the fibre. It solely affects the fibre's inertia due to the fluid immobilized inside the fibre walls. The fluid velocity field was then determined by one parameter G, the shear rate. A fibre wall density, pf  of 1500%/m was used in most trials. The fibre therefore had a 3  w  specific gravity, S of between 1.0 and 1.5, depending on the wall thicknes. The overall dimension of the fibre was set by the outer radius, a . Cylindrical elements had a length and diameter of Q  2a matching those of spherical elements. The fibre length is then equal to 2a iV where N is a  0  the number of elements in the fibre model. Table 3.4 summarizes the constants and parameter ranges used in the trials. A similar table will be shown at the beginning of each section giving the parameters used in that particular trial.  Chapter 3. Results and Discussion  47  Format of Results Data was collected from the model in a Lagrangian fashion. The main data file recorded the position of the fibre. Since the fibre's position is represented by the midpoint of each element, the results displayed here appear half an element shorter than the actual fibre model at each end. Other datafilesrecorded the forces acting in the fibre during the trial. At the end of each trial, a final data file was created which stored all variable values so that the model could be restarted if needed. Each data file was divided into separate zones for each specified time step. In theory, the number of zones is limited by the size of the dimensionless time step, St and the length of the trial. In practice, a zone was written to the data file every 0.1 t* in order to keep the file at a manageable size. This resolution was found to be high enough to observe the behaviour of the model in sufficient detail. The data representations in this chapter were created using Tecplot. This program processes ASCII data files to produce binary plot files allowing the data to be graphed in the manner desired by the user. This software was also used to create animated results of the fibre motion. Tecplot limits each plot file to a maximum of one hundred zones. This limitation required that plot files to be formed by taking zones from a portion of, or at regular intervals from the parent data file for each trial. Results from each trial were analyzed to determine the model's dynamic behaviour. For rigid fibres, the angle <j> formed by the fibre with the y axis was plotted against time and the period TG was determined. The period offlexiblefibreswas determined by inspecting the fibre plots. The degree of deformation seen inflexiblefibreswas found by plotting the orbits formed by the fibre tips as they rotate.  Chapter 3. Results and Discussion  3.2  48  Results for the 10 Element Model  3.2.1  Rigid Body Motion - Series 1  A preliminary series of simulations was conducted to study the rigid body motion of a 10 element fibre model. The purpose of these trials was to: a) reproduce the results of Yamamoto and Matsuoka [22, 23], and b) determine the characteristics of the rotational and translational motion for a rigid model.  Rotational Motion: Trials 1(a) and 1(b) Trials 1(a) and 1(b) were executed using the model developed by Yamamoto and Matsuoka [22]. This model differs from the model developed in the previous chapter in that a solid sphere is used as the basis for the inertial properties of the elements as opposed to a solid cylinder. The hydrodynamic and structural properties of both models are identical. The trials were run at a particle Reynolds number, Re  p  Re  p  — 0.1 This is the same value for  used by Yamamoto and Matsuoka. This number was formulated by choosing the values  for the sphere radius, shear rate, and fibre wall density shown in Table 3.5. The fluid viscosity was set at the value listed in Table 3.4. The elastic coefficient k i was set at 10 . This results 4  e  in a bending constant kb of 3930. This value was well within the range of rigid fibre motion cited in the previous work [22].  100 nm 10 s'  Element Radius, a Shear Rate, G Fibre Wall Density, pf Elastic Modulus, E 0  1  1000 kg/m 100 Pa  3  w  Table 3.5: These values were kept constant in all the trials run in Series 1  The trials were executed for m = 5 X 10 steps using a dimensionless time step of St = 10 . 6  -5  Thus, the motion of the fibre was followed over a dimensionless time, t* of 50 and 500 data zones were collected.  Chapter 3. Results and Discussion  49  The number of time steps, m executed by the code for a particular trial is inversely related to the time step St. The maximum time step allowable was chosen in each trial to reduce m and therefore the computational time required. This was especially true as computational times increased with the higher aspect ratio (60 element) fibres. Trial 1(a) required 15 minutes of processing on a Silicon Graphics "Indy" workstation. The time step was limited by an instability in the model which manifested itself when St went above a critical value. The maximum value for St in each trial was determined by trial and error. The causes of the instability and its relationship to the time step will be discussed at the end of this chapter. The fibre's initial position was aligned with the y axis centred at the origin as shown in Fig. 3.15(a). This resulted in the fibre motion being limited to rotation about the origin. Results for Trial 1(a) are listed in Table 3.6 at four locations from the fibre's initial vertical position: a) 90°, b) 180°, c) 270°, and d) 360°. The values obtained from JefFery's equation with r* = 7.76 are also listed. The value used for the equivalent axis ratio was obtained from Ref.[23] based on the observations of Trevelyan and Mason [21]. Angle 90° 180° 270° 360°  Trial 1(a)  JefFery's Equation  t*  At*  t*  12.5 23.8 36.1 47.4  12.5 11.3 12.3 11.3  12.39 24.78 37.18 49.57  Table 3.6: Results for Trial 1(a). t* refers to the time to complete the specified rotation. At* refers to the time interval over one quadrant of rotation.  The model shows reasonably good agreement with JefFery's equation. It was observed that the model required more time to rotate from the vertical to horizontal position than vice versa (ie. At^o > A i J ) This disagrees with JefFery's equation which predicts equal time over any 80  quadrant of rotation. A delay of t* = 0.2 is apparent in Atg . This may be attributed to the 0  time required by the fibre to come up to speed from its initial angular velocity of zero. The  Chapter 3. Results and Discussion  360 315  50  r  '-  i  i / / / i i 1 i 1  J  270 225 180  (deg) i  i f / / i i i i J i i 1 i J  Trial 1(a)  J  Jeffery's equation  135  90 45 ,  0  .  i  10  .  .  .  ,  i  20  .  .  .  .  i  30  .  .  .  .  i  40  .  .  .  .  i  50  Figure 3.16: One rotation of the ten element model compared with the results obtained from Jeffery's equation using r* = 7.76. "start-up" represents about 1% of the time to complete 180° of rotation and was not considered to be important. Fig. 3.16 shows the plot of these results. Several features of the fibre motion are evident. The pattern of motion repeats every half turn or 180°. This means that, neglecting the "startup" time, the fibre motion may be fully described by following it for 180° of rotation from the vertical. The fibre spends most of the time aligned with the flow direction (ie. 90° or 270°). The model spent 57.2% of its rotation between 80° and 100° or 260° and 280° during one full rotation. For Jeffery's equation, this percentage was slightly higher at 59.9%. The period for Trial 1(a), 47.4 was higher than the value, 45.53 computed by Yamamoto and Matsuoka [23] for a rigid 10 element fibre. A possible cause of this discrepancy was the angular correction procedure outlined in Eqs. 2.25 and 2.26 of Chapter 2. In Yamamoto and Matsuoka's paper [22] the angular position of individual elements was adjusted after the new  Chapter 3. Results and Discussion  51  position had been calculated at each time step. In the method employed here, the correction is applied as a torque before the new position is found for each time step. Trial 1(b) was conducted under conditions identical to Trial 1(a) while neglecting the angular correction routine. This reduced the processing time by 20% to 12 minutes. The results for Trial 1(b) are shown in Table 3.7.  Angle 90° 180° 270° 360°  Trials 1(a) Ai* t* 12.5 12.5 23.8 11.3 36.1 12.3 47.4 11.3  Trial i* 11.2 22.8 34.1 45.6  1(b) Ai* 11.2 11.6 11.3 11.5  JefFery's Equation i* 12.39 24.78 37.18 49.57  Table 3.7: Results for Trial 1(b). i* refers to the time to complete the specified rotation. A i * refers to the time interval over one quadrant of rotation.  The period for this trial, 45.6, was nearly identical to the value published in [23], 45.53. This demonstrated that the angular correction routine employed in Trial 1(a) had a much more significant effect on the results than that employed by Yamamoto and Matsuoka. This observation is reinforced by the values of Ai* which are almost constant through the four quadrants of rotation. The period calculated from JefFery's equation using r* = 7.76 was 8.0% longer than the period of Trial 1(b). It was concluded from this trial that neglecting the angular correction routine had a minimal effect on the motion of the fibre model for the particular conditions tested here. This routine could be neglected to save processing time. This issue was studied further in Trial 1(c).  Chapter 3. Results and Discussion  52  Rotational and Translational Motion: Trial 1(c) Trial 1(c) was conducted under the same conditions as Trial 1(b). However, in this trial the fibre was initially aligned along the y axis extending upward from the origin as shown in Fig. 3.15(b) in order to produce both rotational and translational motion. This initial position was used in all subsequent trials. In order to more accurately represent a real fibre, solid cylindrical elements were used as the basis for the inertial properties of the model as described in the Method chapter. The cylinders had the same length and diameter so that there aspect ratio was 1. The formulation of Re = p a G/rj is the same for both spherical and cylindrical geometries. p  e  0  The change to the cylindrical geometry is reflected in the value of the constants Ct and C tT  TO  Ctr is proportional to the element's mass which increases by 50% for the cylindrical geometry. C t is proportional to the mass moment of inertia of an element and increases by 218% for the ro  cylindrical geometry. The assumption of Stokes flow in this model implies that thefibrecentroid once in motion will translate at the bulk fluid velocity. For the initial position chosen in this trial, the fibre centroid should advance forward at a rate of one fibre length every At* = 2. The centroid should remain at a constant position half afibrelength above the x axis. Fig. 3.17 plots of the results of Trial 1(c). As expected the fibre model rotated at the same rate as in Trial 1(b). Thefibrealso translated at the expected rate of one fibre length every At* = 2 once it attained its forward speed. The centroid of the fibre was found to be approximately 1/100 of afibrelengthbehind its expected position due to the time required by the fibre to come up to speed. The period for Trial 1(c) increased to 50.8 when the correction routine was used. This increase over the period for Trial 1(a), 47.4 was due to the change in C tTO  Chapter 3. Results and Discussion  53  1  \  o  2  o  4  6  8  10  12  Fibrelengths  Figure 3.17: Trial 1(c): A rigid 10 element model translating and rotating through approximately 220°. The fibre centrelines are plotted at intervals of At* = 2 ending with t* = 24. In order to study the effect of the angular correction routine, Trial 1(c) was repeated over a range of Re . p  The constants Ct and C t are proportional to Re for solid elements. The T  ro  p  results of these trials (Series 2) may be found in Appendix E . The conclusions drawn from the results are discussed here. Trials were run in Series 2 for Re ranging from 0.001 to 1.0. The results showed that the p  period of rotation increased dramatically above Re = 0.1 when the correction routine was in p  use. The period for Re  p  = 1.0 was nearly 150. The trials run without the correction routine  remained nearly constant over Re with the period for Re = 1.0 dropping slightly to 44.6. p  p  It was concluded from Series 2 that subsequent trials would be performed without the correction routine due to the problems seen at Re  p  above 0.1. Further development in the  model will be required to create a properly functioning routine. However the results from the preliminary trials indicate that the model performed adequately without the routine. The translational behaviour of the model was unaffected by the angular correction routine in Series 2. It was found that the lag from zero velocity was proportional to the value of CtrFor Re = 0.1, the lag was found to be 0.0165 of a fibre length. For Re = 1, the lag increased p  p  to 0.165 of a fibre length. In summary, the two preliminary series of trials showed that the model performed within reasonable agreement of previous results and the theory for rotational and translational motion.  Chapter 3. Results and Discussion  3.3  54  Results for 60 Element Model  Two series of trials (3 and 4) were conducted using a model constructed of 60 elements. The model was formulated in the same manner as in Trial 1(c); solid cylinders as the basis for the elastic and inertial properties and spheres as the basis for the hydrodynamic properties. The purpose of these trials was to determine the range of motion induced in the fibre model by changing its elastic properties. An aspect ratio of 60 was chosen for these trials as it is representative of typical softwood pulp fibres [19]. The fibre was initially aligned along the y axis extending upward from the origin as shown in Fig. 3.15(b). The Reynolds number, Re in Series 3 and 4 was kept constant at 0.15. The parameters p  used in the formulation of the number are listed in Table 3.8.  10 microns  Outer Radius, a Shear Rate, G Fibre Wall Density, pf 0  1000 s-  1  1500 kg/m  6  w  Table 3.8: These values were kept constant in all the trials run in Series 3 and 4  The trials in Series 3 and 4 were conducted without the angular correction routine. Processing time was reduced by 18% by neglecting the routine. The amount of processing time required by the code became an important consideration in these trials. The 60 element model required hours as opposed to minutes to obtain results. The 10 element model required 2:26 to process 10 time steps while the 60 element model 6  required 15:31. Higher aspect fibres also have a longer period of rotation which increases the dimensionless time, t* for which the motion must be followed.  Chapter 3. Results and Discussion  3.3.1  55  Bending Deformation — Series 3  Series 3 consisted of five trials where the elastic constant, k i was varied from 2 X 10 to 10 6  4  e  to determine the range of motion induced in the fibre model.  Comments will be made in  this section on the qualitative similarities between the model's behaviour and the behaviour of flexible fibres seen in previous work [6].  Trial 3(a) 3(b) 3(c) 3(d) 3(e)  k = E/rjG 2 X 10 10 5 x 10 10 10  6  6  5  5  4  b  K  el  7.85 3.93 1.96 3.93 3.93  x x x x x  10 10 10 10 10  5  5  5  4  3  3.14 1.57 7.85 1.57 1.57  *: x x x x x  St 10 10 10 10 10  6  6  5  5  4  2 x lO" 5 x 1(T 1 x 1(T 5 x 1(T  6 6  5  5  2 x itr  4  Table 3.9: Series 3 Elastic Parameters  It was found in the trials that the dimensionless time step, St had to be reduced as the stiffness of the fibre was increased. Table 3.9 lists the bending and spring constants for each trial, k* and k* along with the corresponding St. Trial 3(a) had to be processed a hundred times longer than Trial 3(e) in order to obtain data for the same period of actual time.  Chapter 3. Results and Discussion  56  Rigid Body Motion - Trial 3(a) The purpose of Trial 3(a) was to observe rigid body motion in the 60 element model. The elastic constant, k i was increased to 2 X 10 . The deflection of the fibre from its initial staight 6  e  position was limited to 0.001 fibrelengths as shown in Fig. 3.19. Angle 90° 180° 270° 360°  Trial 3(a) At* t* 65.6 65.6 133.4 67.8 198.8 65.4 266.6 67.8  JefFery's Equation  t* 62.56 125.11 187.68 250.23  Table 3.10: Results for Trial 3(a). t* refers to the time to complete the specified rotation. At* refers to the time interval over one quadrant of rotation.  The results for Trial 3(a) as listed in Table 3.10. The model required slightly more time to rotate from the horizontal to the vertical than vice versa. This behaviour was consistent with results seen in earlier trials conducted without the angular correction routine. The period resulting from the trial, 266.6, was 6.5 % higher than the value obtained from JefFery's equation, 250.23. This period was calculated using an equivalent ellipsoidal axis ratio, r* = 39.8 obtained from the experimental measurements in Trevelyan and Mason [21]. This suggests that the present hydrodynamic model underestimates the force exerted on the fibre at higher aspect ratios. This could be a result of not considering the interaction of flow fields between neighbouring elements (free draining approximation).  Chapter 3. Results and Discussion  450  57  r  405 360 315  4>  270  r  225  -  (deg) 180 135  90  Trial 3(a) r  Jeffery's equation  F-  45 0  J  50  100  i_  J 150  I  I  1  L_  200  •  i  i  i  i  i  i  i  250  i  i  300  *  t  Figure 3.18: Orientation of rigid 60 element model over 450° compared with results obtained from Jeffery's equation using r* = 39.8. The orientation of the model is plotted in Fig. 3.18. The tendency of the 60 element model to align with the flow was more pronounced than for the 10 element case seen in Fig. 3.16. The motion of the fibre can be described as long periods of translation in a horizontal position interspersed with short "tumbles". The model spent 91.5% of the time within 10° of the flow direction and 83.2% within 5° of the flow direction. The corresponding values from Jeffery's equation were 91.0% within 10° and 82.2% within 5°.  Chapter 3. Results and Discussion  58  Slightly Flexible Motion - Trials 3(b) to 3(c) Trials 3(b) and 3(c) were conducted with elastic constants of k i = 10 and 5 X 10 respectively. 6  5  e  The results showed that the reduced stiffness in Trials 3(b) and 3(c) caused a slight reduction in the period of rotation from Trial 3(a). The periods in Trials 3(b) and 3(c) are 266.2 and 264.2 respectively.  a)t*=133.4  66.6745  b)t*= 133.2  66.6750 66.6755  66.4 66.5 66.6 66.7 66.8  c ) t / = 132.2  65.9 66.0 66.1 66.2  Figure 3.19: Fibre Deflection at cp = 180° for the three stiffest trials: (a) k = 2 X 10 , (b) k i — 10 , and (c) k i = 5 X 10 . The axes are measured infibrelengths. The axes in a) are not to scale in order to magnify the deflection. 6  ei  6  e  5  e  The deflections of thefibresin Trials 3(a) to 3(c) at <p = 180° are exhibited in Fig. 3.19. The time taken to reach <p = 180° is listed above each plot. The distinguishing feature of Trials 3(b) and 3(c) was the bending that occurred in the model. The amount of deflection in Trial 3(b) amounted to 0.184 of afibrelengthand increased to 0.242 of afibrelengthin Trial 3(c). These deflections are due to the increase in the component of fluid forces exerted along the length of the fibre as it rotates from the horizontal to the vertical.  Chapter 3.  59  Results and Discussion  (a)  t =138  t =128 1.0 0.5 0.0  64  65  66  67  68  69  64  65  66  67  68  69  (b)  (c)  Figure 3.20: The motion of the fibre for the three stiffest trials: (a) k and (c) k = 5 X 10 . The axes are measured in fibre lengths.  = 2 X 10 , (b) k 6  ei  = 10 , 6  ei  5  el  Fig. 3.20 plots the fibre centrelines for Trials 3(a), 3(b), and 3(c) as they pass through the "tumbling" region (100° < <p < 260°) between t* = 128 and 140. The fibres are recorded at intervals of A i * = 1 to show the detail of the motion. The motion of Trial 3(b) could be classified as "springy" in the terminology of Forgacs and Mason [6]. A region of bent motion occurs between t* — 132 and 134. The majority of its motion was identical to Trial 3(a). The rotational behaviour of Trial 3(c) could be classified as a "loop turn". Bending increased between t* = 128 and 133. The fibre rotated as a rigid body between t* = 133 and 134 before recovering its straight shape at t* = 13. Due to the flexibility in Trial 3(c), a smaller fibre profile is exposed to the fluid flow in the x direction. This may have significance for the hydrodynamic drag of aflexiblefibre.  Chapter 3. Results and Discussion  60  Very Flexible Motion - Trials 3(d),  3(e)  The elastic constant was reduced to k i = 10 and 10 for Trials 3(d) and 3(e) respectively. 5  4  e  Fig. 3.21 plots the fibre centrelines at ten evenly spaced locations over A i * = 18 during the tumbling region. The dimensionless time at the beginning and end of the region are shown above the fibres. The results from Trials 3(d) and 3(e) differ considerably from the trials in Series 3. Most notably, deformation was more prominent in the fibres and was not confined to the tumbling region but occurred over a wide range of the motion. This made it dimcult to identify the exact angle of rotation of the fibre. In Fig. 3.21 the fibre estimated to represent the (p = 180° position is identified by the time displayed above it. The periods for Trials 3(d) and 3(e) were thus estimated as T G = 210 and 148 respectively compared to 266.6 for Trial 3(a). The motion of the fibre in Trial 3(d) is similar to that seen in Trial 3(c) and could be described as a more compact loop turn. The reduction in the period between the trials is likely due to the increased deformation. In Fig. 3.21(a), it can be seen that the fibre ends were pushed further apart into the flow field by the buckling action around (p — 90°. This induced the tumbling motion in Trial 3(d) that is not seen in Fig. 3.20(c) until the fibre has rotated past <p = 100°.  Chapter 3. Results and Discussion  61  (a) 1.0  t =110  t =98  t =116  0  0.5 0.0 49  50  51  52  53  54  55  56  57  58  (b) 1.0  t=62  t' = 74  t' = 80  0.5 0.0 31  32  33  34  35  36  37  38  39  40  Figure 3.21: The motion of the fibre for the two most flexible trials: (a) k \ = 10 and (b) k i = 10 . The axes are measured in fibre lengths. 5  e  4  e  The motion of Trial 3(e) seen in Fig. 3.21(b) was dominated by fibre deformation. Almost no rigid body motion was visible. The fibre ends bent inwards and migrated toward the middle of the fibre. After passing the midpoint of the fibre at t* = 70 the fibre ends moved apart until the fibre assumed a straight geometry at approximately t* = 88. At this point the pattern seen between t* = 0 and 62 is repeated. The high flexibility of the model in this trial would likely result in a more complex asymmetric orbit for an actual wood fibre with the same properties. However, the symmetry inherent in the formulation of the model causes it to perform in this manner. This symmetry is difficult to create even in highly idealized experimental situations. This represents one of the simplifications of the model which was addressed in the next series of trials. Another weakness of the model is the free draining approximation. The close proximity of the various parts of the fibre in Trial 3(e) would result in some manner of hydrodynamic interaction. The model in its present form does not take this into account.  Chapter 3. Results and Discussion  62  Conclusions for Series 3 Five trials were conducted in Series 3 to determine the range of motion for a 60 element fibre model with stiffness varying from k i = 2 X 10 to 10 . 6  4  e  The rigid 60 element fibre model, Trial 3(a), showed an increase in the tendency to align with the flow direction over that seen with the 10 element model. This behaviour follows that seen in experimental and theoretical work. The period of rotation for Trial 3(a) was 266.6. A slightly decrease was seen in Trials 3(b and 3(c) where T G = 266.2 and 264.2 respectively. Rigid body motion predominated in these trials. When the stiffness was reduced in Trials 3(d) and 3(e) the period dropped to T G = 210 and 148 respectively.  This resulted from the deformation of the fibre which caused the  tumbling motion to occur at an earlier point than in the stiffer trials. The period calculated from JefFery's equation was 250.23 using r* — 39.8. The variation of the period with flexibility is plotted in Fig. 3.22.  -o  250 200  •o  r  150 TG  •  100 h  Results form Series 3  - • Jeffery's equation  50 h 10  s  Elastic Constant (E / n G )  Figure 3.22: Stiffness vs Reynolds Number for the five trials in Series 3.  Chapter 3. Results and Discussion  63  i  •j/ =  3(a)  - - 3(b) 3(c)  i  i  i  i  Li  3(d)  t  3(e)  Figure 3.23: Fibre Orbits for Series 3. The axes are marked in tenths of a fibrelength. It was observed in the trials that the fibre profile exposed to the flow direction during the tumbling region was reduced for the flexible fibres. This phenomena is shown here by plotting the orbits traced by the fibre tips. The orbits were found by having the point of view translate at the expected velocity for the fibre's centroid (ie. the average fluid velocity). The orbits shown in Fig. 3.23 are shifted slightly to the left due to the translational lag of the fibre from zero velocity. The maximum profile of thefibreexposed to the flow direction as a percentage of fibrelength decreases from 100% for the rigidfibrein Trial 3(a) to 90.2, 57.7, 32.4, and 17.1% for Trials 3(b) to 3(e) respectively. This would suggest that the most flexible fibres have the lowest profile. However, if the profile is averaged over a whole rotation it is likely that the mostflexiblefibres will have a higher profile due to the deformation that exists throughout the rotation. Stiffer fibres would have a much smaller profile due to the large proportion of time they spend aligned with the flow. Another consequence of the orbits is the smaller effective area occupied byflexiblefibres while rotating. This could reduce the possibility of interaction between fibres in a suspension of flexible fibres in a simple shear flow compared with rigid fibres at the same concentration.  Chapter 3. Results and Discussion  3.3.2  64  Localised Wall Weakness - Series 4  Series 4 consisted of five trials where a localized weakness (kink) in the bending stiffness was placed at five evenly spaced locations along the fibre. These trials were conducted to address the deficiencies of the symmetric formulation of the model. The conditions for these trials were the same as Trial 3(c). The parameters used in this series are listed in Table 3.11. Thefibrekinks were modelled by reducing the bending stiffness at one of 5 locations along the fibre for each trial. The bending stiffness, k* at the kink was reduced to 20% of the global value. 10 \xm 1000 s~  Element Radius, a Shear Rate, G Wall Density, pj  1500 kg/m  E/ G  5 x 10  0  w  V  l d  5  Table 3.11: These values were kept constant in the trials run in Series 4  The results for the trials are shown in Fig. 3.24 for the tumbling region between t* = 128 and 140. It can be seen from the figure that Trials 4(c) and 3(c) yielded similar results. The period of all the trials in Series 4 was unchanged from Trial 3(c). Figs. 3.24 shows that a kink caused the model to reproduce a "snake turn" as described in the literature as opposed to the symmetric "loop turn" seen in Trial 3(c). The motions in Trials 4(a) and 4(e) were identical but 180° of phase. The same was true for Trials 4(b) and 4(d). The kink in Trial 4(a) allowed one end of the fibre to buckle first but the bend did not occur at the kink as was the case in Trial 4(b). The importance of these trials lies in demonstrating that the model can simulate the actual motion seen in wood fibres by applying a local modification to the bending stiffness. Further tests could be run at the lower bending stiffnesses used in Trials 3(d) and 3(e) to determine if the periods of the rotation found in these trials are effected by a kink.  Chapter 3. Results and Discussion  65  Figure 3.24: Fibre motion over the tumbling region for the kink located at: a) 1/6, b) 1/3, c) 1/2 d) 2/3, and e) 5/6 from the left end of the fibre. The arrow points to the location of the kink. The axes are measured in fibrelengths.  Chapter 3. Results and Discussion  3.4  66  Discussion  To this point, results have been presented separately for fibres of different aspect ratios with rigid orflexiblemotion. A single dimensionless number is required which can identify the range of motion in a fibre for any fluid or fibre conditions. A simple dimensional analysis is used in the development of this number below.  Fibre Bending Number, B First, consider the viscous forces exerted on the fibre due to the shear flow,  where n is the viscosity, a is the outer radius of the fibre,and G is the shear rate. The viscous c  force per unit length on each element is then,  ^  oc  Ga  V  0  Now, consider the deflection associated with the buckling or bending of the fibre where it is considered a beam or column,  0  oc  wL* EI  ——  where 8 is the deflection, w is the load per unit length, L is the length of the fibre, and EI is the stiffness. For a critical deflection, (^) the load per unit length is,  EI W  K  V  Chapter 3. Results and Discussion  67  A dimensionless bending ratio, B can then be formed,  B  =  ^  (3.4!)  nGa  0  = &&T?  '  (3 42)  Hence, a fibre is stiff (will not deform) if B exceeds a critical value. A more exact ratio for B may be formed by considering the theory of Forgacs and Mason [13]. This theory, developed in the Introduction, considers the case of a fibre where the viscous forces due to the fluid vary along its length. Eq. 1.8 is repeated below.  E  2r  %  >G-n  ,crit  4  fn(2r )-1.75 c  where r is the aspect ratio of the cylindrical fibre and a-a refers to the value above which the c  fibre will remain rigid. This equation can modified to include the effect of a hollow fibre on the stiffness EI by adding an additional term involving J . After rearranging terms, a new form of the bending ratio, B is found,  As in the previous ratio, a fibre is stiff if B exceeds a critical value. Both of these ratios are made up of three terms. The first is a ratio of the fibre stiffness to the fluid viscous forces. The second measures the effect of hollowness on the second moment of area, / . The third term is a measure of the aspect ratio effect. It was determined from the results of the 60 element model that the period of rotation decreased with stiffness for flexible fibres. The period is a useful manner of measuring the flexibilty of a fibre. If the flexible period is normalized with respect to the period for a rigid fibre, TG/(TG) igid, r  results may compared over different aspect ratios.  Chapter 3. Results and Discussion  68  An additional series of trials was conducted with the 10 element model over a range of stiffness k i = E/r/G ranging from 50 to 10000. This allowed the flexible motion to be compared e  for both aspect ratios N=10 and 60. The data for the N=60 fibre comes from Series 3. The results of the trials are tabulated in Table 3.12 and 3.13.  kl e  50 100 200 500 1000 2000 4000 10000  B (Eq. 3.42) 6.25 X IO" 1.25 X IO" 2.50 x IO" 6.25 x 10~ 0.125 0.250 0.500 1.25  3  2  2  2  B (Eq. 3.43) 3.11 X IO" 6.23 X IO" 1.25 x IO" 3.11 x I O 6.23 x 10~ 0.125 0.250 0.623  3  3  2  - 2  2  TG 39.7 40.4 41.0 43.3 45.0 45.5 45.6 45.6  % of TG igid r  87.1 88.6 89.9 95.0 98.7 99.8 100 100  Table 3.12: Results for a flexible 10 element fibre and the two calculated bending ratios  kl e  10 10 5 x 10 1 x 10 2 x 10 4  5  5  6  6  B (Eq. 3.42) 5.79 X IO" 5.79 X IO" 0.289 0.579 1.16  3  2  B (Eq. 3.43) 1.17 X IO" 1.17 x IO" 5.86 X IO" 0.117 0.234  3  2  2  TG 148 210 264.2 266.2 266.6  % Of TGrigid 55.5 78.8 99.1 99.8 100  Table 3.13: Results from Series 3 and the two calculated bending ratios  TG/(TG) igid T  number B.  versus B is plotted in Fig. 3.25 for the two different forms of the bending  Figure 3.25:  TG/(TG) g ri  id  versus B a) Eq. 3.42  and b) Eq.  3.43  Chapter 3. Results and Discussion  70  It appears from the plots that the form of the bending ratio in Fig. 3.25 (b) collapses the data from the two aspect ratios more effectively. In this plot the two vertical lines divide the range of motion observed into three regions. In region I the fibre is rigid and tyhe period remains constant. In region 77, the fibre deforms but the period is reduced by only a small amount. This could be described as "springy" motion. Region III indicates the range in which the period is reduced by a significant amount. The fibre performs a "loop turn" in this region. The most important observation from these plots is that the location of the regions appears to be independent of the aspect ratio as the theory suggests. Another feature of the plots is the slow change in the period for the the first three data points of the N=10 fibre. This was likely a result of the axial stretching which was observed in the fibre. The development of a fibre bending number and the determination of critical values for this number is a new manner in which to to analyze the problem of fibre motion in a shear flow. If, for given fibre and fluid conditions, the fibre motion lies in the rigid range, a simple model based on JefFery's equation can be used to describe its motion. Otherwise a more complex model, such as the one developed in this thesis, would have to be used.  Chapter 3. Results and Discussion  71  Practical Application The results from the model have been presented in a dimensionless form. Some comments will now be made on the relation of the results to the properties of fibres and flow fields seen in practice. The properties of softwood pulp fibres were discussed in the Introduction. The range of dimensions for pulpwoods found in British Columbia are outlined here according to [19].  Fibre Parameters Diameter (microns) Wall Thickness (microns) Moment of Area, I (m ) Length (mm) Aspect Ratio Stiffness, EI (Nm ) 4  cs  4  C3  Lower Limit 30 2 1.73 x I O < 2 40 1 X IO" (Kraft) - 2 0  12  Upper Limit 50 10 2.67 x 10~ 4 100 100 x I O (TMP) 19  - 1 2  Table 3.14: Variation of fibre parameters for softwood pulp in British Columbia. Kraft and T M P (thermo-mechanical pulp) refer to the pulping processes used.  Reynolds Number The flow field in which the fibre model is immersed is defined by two parameters: shear rate and viscosity. The shear rate of flow fields seen in actual hydrocyclones is in the order of 100 5  _ 1  . The region around a slot in a screen may see shear rates in the order of 1000  The  viscosity used in the trials was that of water at 20 °C, 0.001 kg/ms. In actual fractionation processes, the concentration of fibres may be high enough to increase the effective viscosity of the fluid due to the interaction of the flow fields of neighbouring fibres. This model considers only isolated fibres. This model is based on the assumption of low relative velocity between the fibre and the fluid. This assumption was tested by recording the relative velocity at each element for the two rigid cases where N = 10 and 60. For the majority of the motion, where the fibre was closely aligned with the fluid vectors,  Chapter 3. Results and Discussion  72  the relative velocities were small such that actual the Reynolds number was much less than 1. The largest relative velocities occurred at the tips when the fibre was aligned at 45° to the flow during the tumbling motion. For the N = 10 case the maximum dimensionless relative velocity was 3.84. For the 60 element case, the value was approximately 6 times larger at 23.6. This suggests that the dimensionless relative velocity is proportional to the aspect ratio. Consider the case for the maximum Reynolds number (a = 25 microns, G = 1000 s , -1  c  relative velocity = 23.6). The actual relative velocity is,  Ui re  =  23.6(25 x 10 m)(1000 s )  (3.44)  =  0.59m/s  (3.45)  _6  _1  The Reynolds number is then calculated using standard values for the density and viscosity of water,  Re  =  (3.46)  _  (0.59)(25xl0- )(1000)  =  14.75  6  (3.48)  This value is well outside the range of Stokes flow. However, it should be emphasized that the Reynolds number only reaches this maximum value for a very small portion of the fibre's motion.  Fibre Bending Number Consider the case for the minimum bending number where a = 25 microns, G = 1000 s , _1  D  and EI = 1 X 10 JVm . For an aspect ratio of 100, as calculated by Eq. 3.43, B is equal to _12  4  0.058. This lies near the transition between region II and III on Fig. 3.25 indicating that an isolated fibre would begin to deform significantly under these conditions.  Chapter 3. Results and Discussion  3.5  73  Recomendations for Improvements to the Model  Wall Thickness T h e m o d e l as presently formulated uses a cylinder as the basis for the elastic a n d i n e r t i a l properties of t h e fibre. T h e varying w a l l thickness of real w o o d fibres is represented i n the m o d e l t h r o u g h t h e ratio of inner t o outer radius i n the fibre, a * = T h e a x i a l stiffness,  EA  CS  and the bending stiffness,  EI  CS  ai/a . 0  of the fibre m o d e l are a f u n c t i o n  of the Y o u n g ' s modulus and the cross section of the fibre w a l l calculated f r o m a*. T h e present f o r m u l a t i o n assumes that the fibre w a l l strains linearly. In reality, pulp fibres behave differently whether stretched or bent due t o their complex geometry of various layers and fibril angles. It was also f o u n d i n the trials that varying the a x i a l stiffness of a fibre h a d no significant effect on its d y n a m i c behaviour for the range of stiffnesses tested. T h e available d a t a for pulp fibres usually specifies the bending stiffness, EI , the outer CS  diameter, a , and the density, p . a  p  These considerations suggest some simplifications that could be made i n the m o d e l . F i r s t , the calculation of a x i a l forces i n the fibre could be neglected. T h e y could be replaced b y a nonslip force, similar t o that which acted tangentially between the elements i n the present m o d e l . T h e new non-slip force would have components i n b o t h the a x i a l and tangential directions t o ensure t h a t neighbouring elements do not break continuity. T h e w a l l thickness as presently defined makes the formulation of several parameters more complicated t h a n necessary. It is recommended that the properties of the m o d e l be b r o u g h into line w i t h the available d a t a . T h i s would remove the ratio a* f r o m the f o r m u l a t i o n of t h e fibre model's properities.  Chapter 3. Results and Discussion  74  Stability of the Model It was discovered in the trials that certain conditions caused the numerical code which defined the model's behaviour to become unstable.  This instability was manifested initially when  adjacent elements moved apart so that the fibre centreline became crooked. The phenomena then developed to the point where adjacent elements would move apart to distances that were orders of magnitude larger than the fibre itself. The instability occurred under two conditions. For a given elastic constant, k i, there existed e  a critical time step, St above which the instability occurred. The relation between k l and St e  varied with the aspect ratio of the model. Instability also developed when the fibre was not initially aligned vertically. The factor which caused the limitation in both cases was the axial force term, F" in the equations of motion. This was the first term to grow without bound. This problem would likely be reduced by neglecting the axial forces as suggested in the previous section. This would allow the time step to be increased making the model less computationally intensive.  Reducing Processing Time Several steps could be taken to improve the efficency of the numderical code. An adaptive time step could be implemented to allow larger steps to be taken while the fibre is nearly aligned with the flow. This could reduce computational times significantly. Using elements with a higher aspect ratio would also achieve the the same goal. This would require reformulation of the inertia! and hydrodynamic forces in the model.  Chapter 4  Conclusion  A numerical model of fibre motion in shear flow has been developed. A series of trials have been conducted to determine the behaviour of the model over a range of parameters. In particular, the effect of model's stiffness on its motion was studied. The following conclusions can be made from these trials: • Rigid Rotation The model was deemed to be rigid if its deflection was less than 1% from linear. The period of rotation for the N=10 fibre was 45.6. This was 8.0% less than the value obtained from JefFery's equation using r* =7.76. This value for r* was used by Yamamoto and Matsuoka [23]. The period for the N=60 model was 266.6 which was 6.5% more than the value obtained from JefFery's equation using r* = 39.8. This value for r* was obtained from experimental measurements in Trevelyan and Mason [21].  The period increased  linearly with the aspect ratio of the model. • Translation The fibre translated at the average velocity of the fluid over its initial position. The lag resulting from the initial zero velocity was proportional to C . For the values of C t r  t r  used  in the trials the lag was insignificant. • Flexible Rotation and Kinks For a solid 60 element fibre, the period was constant at the rigid value for k i > 2 X 10 . 6  e  For 1 X 10 > k i > 5 X 10 , the period was reduced by less than 1%. For k i < 5 X 10 , the 6  5  5  e  e  period reduced at a rate proportional to log(C i). The fibre model deformed in a manner e  75  Chapter 4. Conclusion  76  similar to that seen in previous experimental results for pulp fibres. Kinks or localized weaknesses in the fibre walls had negligible effect on the period of rotation. However, The motion of kinked fibres while rotating was asymmetric. • Fibre Bending Number B A dimensionless ratio B was derived to determine the type of motion in the fibre model depending on the fluid and fibre properties. Comparisons with results showed that there existed a critical value for B, independent of aspect ratio, above which the fibre behaved rigidly. The behaviour of the model in the trials exhibited several areas which could be improved. These were discussed in detail in the previous chapter. Several recommendations for the future development of the model are made: • Use Longer Elements The elements used in the model had the same length and diameter. Longer elements could be used to reduce the computational time requirements. The hydrodynamic and inertial models of the fibre would need to be reformulated to reflect this change. • Remove Spring Force, F" This force should be replaced by an axial component of the non-slip force as described in the previous chapter. This should reduce the stability problems which were found in the model. Neglecting stretching in the fibre would have a neglible affect on the model's dynamic behaviour for the range of values tested in these trials. • Simplify Parameters EI  C3  and p  e  The present formulation of these parameters makes them a function of a*, the wall thickness ratio. This formulation overly complicates these parameters. The model would be simplified by specifying these parameters in the same manner as data published for fibre properties.  Appendix A  Inertial Properties of a Fibre Element  Mass, m  e  m-element  m  e  " V i b r e wall  =  = m,f  = A J pf =  0  7r ( a  "Vluid  + mf  w  c  +  + Afl pf  w  2  -  0  a?)  (2a )  p  0  7ra  +  fw  2  (2a ) 0  2  =  27ra (1 - -a | - ) ^ / u , + 2 7 r a a /?/  =  27ra [ p  3  2  0  l  2  -  3  fw  -± (p  - p) }  fw  f  o  a  Density, p and Specific Gravity, S e  e  Pe  =  Pfw ~ \  (Pfw  ~  Pf)  Pf C  Jfw  Pf  _  w  -  Pf  a, ,pfw _ 2  Pfw  CX  1  x  2  -S/u; -  (Sfw -  a?.  «2' 77  1)  a?  Appendix A. Inertial Properties of a Fibre Element  Mass Moment of Inertia, I  e  lelement Ie.  Ifw  » outer  -  —  Izz  =  I fibre, wall  =  Ifw  +  ^  lo Pfw) ( 3 d + All)  2  2  2  {\<  (2ira p ) 0  1_  ~  ^  {\dj  (3a + 4a ) -  ^  (27ra, a p )  2  3  12  If  rrii ^ ( 3 a * + 4Z )  + 4Z ) -  rrif  Ifluid  h r  ^ ( 3 ^  12  +  fw  2  [27ra ( a - a?)  l  p )  a  fw  2  0  (3a + 3a + 4a )  2  2  0  2  2  (7a + 3a ) 2  2  =  Ie = ^  ^(3a  2 t  (3a + 7 a ) + ^ 2  t  2  + 4a ) 2  (3a + 4a ) 2  2  fw  (3d + 2  412)  (3aj + 4 a ) 2  Appendix B  Structural Properties of the Fibre Model  Spring Force Equation  (Ee) A  ca  EA — <° 2a  Al  EA  0  =  k Al a  For a hollow cylindrical rod, the cross sectional area is:  = *(«2-«?) = ™o(i-4) The spring constant, k is then: s  *.  =  7T „  2  a?  -Ea (l-^) . al v  79  0  Appendix B. Structural Properties of the Fibre Model  80  Bending Torque Equation  T  b  =  EI  =  EL  r  dd dx AO lo  EI , A6 2a C  0  =  k AO b  For a hollow cylindrical rod, the second moment of area, I  cs  The bending constant, kb is then:  k = ^Ea\ b  (1  -  about the neutral axis is:  Appendix C  Dimensionless Parameters and Equations of Motion  Dimensional Parameters Unit Time Length Velocity Acceleration (translational) Acceleration (rotational) Force Torque  81  Parameter  Value  1/G a aG aG G alvG a nG  s m m/s m/s  a  D  2  a  2  2  0  kg m/s kg m /s  2  2  2  Appendix C. Dimensionless Parameters and Equations of Motion  Dimensionless Equations of Motion Translation x Direction  du m — dt e  Am*  du"  Similarily for the T/ direction  Rotation  r _ <ft e  (0) (T * + E * + E /*• x « « + E h  dt*  TB  Appendix C. Dimensionless Parameters and Equations of Motion  Hydrodynamic Terms  F  -  -6irr)a (u - V)  F*  =  -67T77a (tt - Gy)  =  -6TT (u-Gy)  =  - 6 T T (ti* -  h  h  0  (^jg)  0  (—p;)  y*)  Sirr/a , (UJ — fi) 3  -8*r,al(u +  rph*  G/2)(—)  -8* ( « + f ) ( ^ ) =  -8?r («* + 1/2)  Appendix C. Dimensionless Parameters and Equations of Motion  Structural Terms F'  = -k  F*  =  0  s  . a. 2^(1-^) 2  7r  2 ~  2a )  = -k* {I* - 2)  3  k,  (I -  s  a r/(j  2 77G  T  a  0  =  6  T*  a2  U  0  j  -ft AO t  = -k* AO  h  H = §*:ci-$ v  -  *  E  n  "<i  84  Appendix D  Velocity Boundary Condition  Generalized Coordinates  dr{  _  + 0 o  d6{  _  X n  dOj  d-Ti  ,. _ =  +  a o  _  X 7 l  .  i  Generalized Non-Dimensionalized Coordinates  dr*  dOi  dt*  dt*  3  dr*  d0j  dt*  dt*  3  Component of velocity perpendicular to 71^+1 (* is assumed)  Ui X nyi — Vi X nxi — oji = v.i i X nyi — Vi+\ X nxi + tv^+i +  Additional Velocities due to unknown force f = F X St t  X nyi-  t  / * _ ! X 7l1/i_i)/C'trans  Ui  =  Ui + {fi  Vi  =  Vi - (/,• X nXi - / / _ ! X  "i  =  Ui-(fi+fU)/C t ro  85  nXi-x)IC rans t  Appendix D. Velocity Boundary Condition  86  Application of the Condition at Interface i —> i + 1  [«i + (/' X.njfc - /i-i X nW-iJ/Ctron,] ~  [i  ~ (fi  v  X  n  x  i  ~ fi-1  X  n x  Xny  {  XTIX;  i-l)/ tran ] C  3  — cy? -»=  [u i + (/ '  -  [v  i+  i+1  t  +1  x ny  - ft x nyj/Ctrans]  i+1  - (fl X nx +1  +  + ft x n x i ) / C  i+1  xnyi  4rans  ] Xnxi  [ui i-{fi+fU)ic \ +  TOt  Collect Coefficients of /*  fi-i /' fi+l  X nyi - nxi-i x nxi)/C'trans  [(-nyi-i  [{ Vi + 2 n x ? ) / C - nXi X  n  = (itii - Ui) X  /i-i  trans  nX )/Ctran i+1  X nyi - nxi_i X n x i ) / C [2/C  /• [{-nyi  X ny  i+1  trans  - nxi  = {u i - Ui) X nyi - (v i+  + 1/Crot]  3  a  n a  l/C  +  r o t  ]  +1  t  + l/C  - Vi) X nxi +  fi-i X bbi + ft x ddi +  + +  rot  nxi )/C rans  +  + a;;)  + 2/C ]  X  i+1  t T  +  +  rot  - Vi) X nii +  -  +  r o  + 2/C ]  2n  [{- yi X  l/C t]  +  r o t  ]  x aai = b  + Ui)  {  +  Appendix D. Velocity Boundary Condition  87  A Tridiagonal Array is formed,  dd aa  0  •  bb dd aa  0  0  bb dd aa 0  0  •  0  0  bb dd aa  0  0  bb dd aa  •  0  bb dd  Special cases occur at each end of the fibre model: element 1 and element m  Appendix E  Angular Correction Routine versus Re - Series 2 p  A series of trials was conducted to determine the effect of the angular correction routine on the motion of the fibre model over a range of particle Reynolds numbers, Re . All parameters p  except pf  were kept constant from Trial 1(c). Re was varied by adjusting from 10 kg/m  3  w  p  to  10000 kg/m . 3  As in Series 1, the behaviour of the model was determined by the time required to complete 90°, 180°, 270° , and 360° of rotation. Thefibrewas initially aligned along the y axis extending upward from the origin as shown in Fig. 3.15(b). Six trials were conducted with the angular correction routine in place for values of Re  p  ranging from 0.001 to 1. Three trials were completed without any angular correction for Re  p  = 0.001, 0.1, and 1. Trial  2(c) IO"  2(d)  2(e)  2(f)  Trial  2(h)  2(i)  0.30  1.0  Re  2(g)  0.15  10"  io-  1.0  11.6  14.5  17.1  25.3  67.0  '90  11.4  11.3  10.0  11.4  11.4  16.8 83.8  Ai*  11.4  4*  22.8  11.5 22.8  12.5  23.0  12.6 37.9  57.8  Ai*  11.4  11.2  141.6  '270  34.2  34.0  9.6 32.1  16.8 158.4  Ai*  11.4  360  45.6  11.6 45.6  12.5 44.6  Re  2(a) IO  10~  '90  11.5 11.4 22.9  -3  p  Ai*  2(b) 2  1  25.9  11.5 28.6  11.5  13.6  15.0  +*  11.5 34.4  34.5  39.5  43.6  19.5 57.4  Ai*  11.4  11.4  11.3  ••360  45.8  45.9  50.8  11.5 55.1  12.6 70.0  l  180  Ai* '270  p  l  l  180  3  1  22.5  Table E.15: Results from Series 2. Trials 2(a) to 2(f) used the angular correction routine. Trials 2(g) to 2(i) were conducted without it.  The results for all trials are listed in Table E.15. The baseline test, Re  p  = 0.001 showed a  consistent rate of rotation over each quadrant for both trials 2(a) and 2(f). Both trials yielded  88  Appendix E. Angular Correction Routine versus Re - Series 2  89  p  160  r  Figure E.26: Plot of dimensionless period, T G versus Re for trials in Series 2 p  nearly identical results. As Re was increased above 0.1, several trends were visible from these results. Trials 2(c) p  to 2(f), with the correction routine intact, exhibited a rapid increase in the time to complete the initial quadrant of rotation, i*, , and in rotating from the vertical to horizontal, At*, . The 0  70  trials conducted without the routine, 2(g) to 2(i) showed a decrease in the time to complete these quadrants of rotation. In all trials the time required to rotate from the horizontal to the vertical increased moderately as Re increased. p  These factors combined to cause the period to increase rapidly for trials 2(a) to 2 (f) as Re  p  increased above 0.1. Trials 2(g) to 2(i) showed a slight decrease in the period. The variation of period versus the Reynolds number is plotted in Fig. E.26. The period was found by calculating TG = 2(£*. - i j ) . This ensured that the rotational delay due to the inital conditions did not 60  80  have an effect on the value.  Appendix F  Code listing  ************ MAIN PROGRAM ************* c********************************************************************** c implicit double precision (a-h,o-z) parameter (m=10) c common aa(m-l),bb(m-l),dd(m-l),b(m-l) common + uO(m),vO(m),rotvelO(m),accelx(m),accely(m),rotacc(m), + fs(0:m),torbz(0:m),alpha(0:m),spr(m),fh(m),ben(m),th(m) common / l l / npi(m) common /xyvel/ + x(m),y(m),u(m),v(m),angle(m),rotvel(m),tg(0:m),ft(0:m) + ,cnx(0:m),cny(0:m), cks(0:m), ckb(0:m) c  c c c c c  c  c  Assign i n i t i a l constants mml=m-l ao radrat g rhop rhof E eta shrate pi  = 100.0d-6 = 0.0d+0 = 9.81d+0 = 1.0d+3 = 1.0d+3 = 100.0d+0 = 1.0d-3 = 10.0d+0 = 3.14159d+0  area arin amf amp amc amin  = = = = = =  pi*(ao**2)* (1.0d+0-radrat**2) pi/4.0d+0*(ao**4)* (1.0d+0-radrat**4) 2.0d+0*pi*(ao**3)*(radrat**2)*rhof 2.0d+0*pi*(ao**3)* (1.0d+0-radrat**2)*rhop amp + amf amp*(ao**2)/12.0d+0* (7.0d+0 + 3.0d+0*radrat**2) + + amf*(ao**2)/12.0d+0* (4.0d+0 + 3.0d+0*radrat**2) rhoc = rhop - radrat**2*(rhop-rhof) re = (ao**2)*rhoc*shrate / eta cksl = pi/2.0d+0 * E / (eta*shrate) * (1 - radrat**2) ckbi = pi/8.0d+0 * E / (eta*shrate) * (1 - radrat**4) ctrans = 2.0d+0*pi*re crot = amin*shrate / (ao**3*eta) mm i s the number of iterations  90  Appendix F. Code listing  c c c  c  c  t = O.Od+0 delta = 1.0d-05 mm = 5.Oe+6 I n i t i a l Conditions f o r Elements 1 - m x(l)= 0.0d+0 y(l)= 1.0d+0 dx=2.Od+0 dy=2.Od+0 do i=2, m x(i)=x(i-l) + 0.0d+0 * dx y(i)=y(i-i) + 1.0d+0 * dy enddo do i=l,m u(i)= 0.0d+0 v(i)=0.0d+0 angle(i)= 0.0d+0 rotvel(i)= 0.0d+0 npi(i)= 0 enddo  c c c c  c c c  c c c c c c c c c c c c  I n i t i a l i z e variables to 0 do i=0,m cnx(i)=0.0d+0 cny(i)=0.0d+0 alpha(i)=0.Od+0 fs(i)=0.0d+0 torbz(i)=0.0d+0 ft(i)=0.0d+0 tg(i)=0.0d+0 enddo I n i t i a l i z e spring and bending constants do i=0,m cks(i) = cksl ckb(i) = ckbl enddo Specify kink location ckb(il) = 0.2d+0*ckbl Open Data F i l e s open (unit=31,file='main.dat',status='unknown') open (unit=41,file='fs.dat',status='unknown') open (unit=5i,file='torb.dat',status='unknown') open (unit=61,file='ft.dat',status='unknown') open (unit=7i,file='torg.dat',status='unknown')  91  Appendix F. Code listing  c c  c c c c c c c c  Input values from previous t r i a l ( i f needed) open(unit=8,file='flag.in',status='old',form='formatted') read(8,*) i f l a g close(8) if(iflag.eq.1) then open(unit=9,f ile='final.dat',status='old',form='formatted') read(9,*) x,y,angle,alpha,npi,u,v,rotvel,t close(9) end i f Record I n i t i a l Position of Elements write(31,*) 'zone' do i=l,m write(31,ll) ao*x(i),ao*y(i),t/shrate enddo format(dl6.8,2x,dl6.8,2x,dl0.5)  11  c c c c  Begin time do loop do 500 ii=l,mm  c c c c  t= t + delta Calculate normals, spring & bending forces  +  c c c  Calculate accelerations  + + c c c  do i=l,mml r r = dsqrt((x(i+l)-x(i))*(x(i+l)-x(i)) +(y(i+l)-y(i))*(y(i+l)-y(i))) cnx(i) = ( x ( i + l ) - x ( i ) ) / r r cny(i) = ( y ( i + l ) - y ( i ) ) / r r f s ( i ) = - cks(i)*(rr-2.0d+0) torbz(i) = ckb(i) * (angle(i+l)-angle(i)) enddo  +  and update velocities  do i=l,m accelx(i) = - 6.0d+0*pi* (u(i)- y(i))/ctrans - (cnx(i)*fs(i)-cnx(i-l)*fs(i-l))/ctrans accely(i) = - 6.0d+0*pi*v(i)/ctrans - (cny(i)*fs(i)-cny(i-l)*fs(i-i))/ctrans rotacc(i) = (-8.0d+0*pi* (rotvel(i)+0.5d+0) + torbz(i)-torbz(i-l))/crot  store velocities from start of time step and calculate new values u0(i)=u(i) v0(i)=v(i) rotvelO(i)=rotvel(i) u(i)=u(i)+accelx(i)*delta v(i)=v(i)+accely(i)*delta rotvel(i)=rotvel(i)+rotacc(i)*delta  Appendix F. Code listing  c c c c c c c c c c c c c c c c c c c c  enddo output bending and spring values if(mod(ii,mm/1000).eq.O) then write(41,*) 'zone ' write(51,*) 'zone ' do i=l,m f h ( i ) = -6.0d+0*pi* (u(i)- y ( i ) ) spr(i) = -cnx(i)*fs(i)+cnx(i-l)*fsXi-1) t h ( i ) = -8.0d+0*pi* (rotvel(i)+0.5d+0) ben(i) = t o r b z ( i ) - t o r b z ( i - l ) write(41,ll) s p r ( i ) , f h ( i ) , t write(51,ll) ben(i), t h ( i ) , t enddo endif f i n d Non-Slip force define tridiagonal terms  + c 1 c  93  1  1 1 1 c 1 1  do i=l,mml b(i)= -((u(i+l)-u(i))*cny(i)-cnx(i)*(v(i+l)-v(i)) + (rotvel(i+l)+rotvel(i))) enddo dd(l)=-2.0d+0*cnx(l)*cnx(l)/ctrans -2.Od+0*cny(1)*cny(1)/ctrans-2.Od+O/crot aa(l)= cnx(l)*cnx(2)/ctrans + cny(2)*cny(l)/ctrans - 1.Od+O/crot do i=2,mml-l dd(i) =-2.0d+0*cny(i)*cny(i)/ctrans - 2.Od+0*cnx(i)*cnx(i)/ctrans-2.Od+O/crot bb(i) = cny(i-l)*cny(i)/ctrans + cnx(i)*cnx(i-l)/ctrans - 1.Od+O/crot aa(i) = cny(i+l)*cny(i)/ctrans + cnx(i)*cnx(i+l)/ctrans - 1.Od+O/crot enddo i=mml dd(i) =-2.0d+0*cny(i)*cny(i)/ctrans - 2.0d+0*cnx(i)*cnx(i)/ctrans-2.Od+O/crot bb(i) = cny(i-l)*cny(i)/ctrans + cnx(i)*cnx(i-l)/ctrans - 1.Od+O/crot  c c c  c a l l subroutine tridiagonal  c c c  define f t from solved non-slip  c c c  update v e l o c i t i e s from f t  call  tridiagonal(1,mml,bb,dd,aa,b)  do i=l,mml ft(i)=b(i) enddo  condition  Appendix F. Code listing  c c c c c c c C  c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c  c c c  c  do i=l,m u(i)=u(i) + cny(i)*ft(i)/ctrans - c n y ( i - l ) * f t ( i - l ) / c t r a n s v(i)=v(i) - cnx(i)*ft(i)/ctrans + c n x ( i - l ) * f t ( i - l ) / c t r a n s rotvel(i)=rotvel(i) - f t ( i ) / c r o t - f t ( i - l ) / c r o t enddo f i n d another tg from geometric conditions do i=l,mml alold=alpha(i) if(abs(x(i+l)-x(i)).lt.l.0d-10) then alpha(i)=-dble(npi(i))*pi ©Xs 6 alpha(i)=-(pi/2.0d+0-datan((y(i+l)-y(i))/(x(i+l)-x(i)))) + -dble(npi(i))*pi endif if(abs(alpha(i)-alold).gt.pi*2.0d+0/3.0d+0) then alpha(i)=alpha(i)-pi npi(i)=npi(i)+l endif tg(i)=ckbl/crot*(angle(i+l)+angle(i)-2.0d+0*alpha(i)) enddo do i=l,m rotvel(i)=rotvel(i)-tg(i)*delta-tg(i-l)*delta enddo output f t & tg  values  if(mod(ii,mm/1000).eq.O) then write(61,*) 'zone ' write(71,*) 'zone ' do i=l,m spr(i) = c n y ( i ) * f t ( i ) - c n y ( i - l ) * f t ( i - l ) ben(i) = - t g ( i ) * d e l t a - t g ( i - l ) * d e l t a write(61,llj s p r ( i ) , t write(71,ll) ben(i), t enddo endif Calculate new displacements and rotation do i=l,m x(i)=x(i)+0.5d+0 *(uO(i)+u(i))*delta y(i)=y(i)+0.5d+0 *(vO(i)+v(i))*delta angle(i)=angle(i)+0.5d+0 *(rotvelO(i)+rotvel(i))*delta enddo Record displacements and time if(mod(ii,mm/500).eq.O) then write(31,*) 'zone ' do i=l,m write(31,ll) x(i)/20.0d+0, y(i)/20.0d+0, t enddo endif  94  Appendix F. Code listing  c c c  c c c c c c c c c  c c c  c  500  cont inue  Output f i n a l values open(unit=9,f ile='f inal.dat',form='formatted') write(9,*) x,y,angle,alpha,npi,u,v,rotvel,t close(9) Close Data Files close(31) close(41) close(51) close(61) close(71) stop end subroutine t r i d i a g o n a l ( i l , i u , bb, dd, aa, cc) double precision aa, bb, cc, dd dimension aa(iu), bb(iu), dd(iu), cc(iu) lp=il+l do i = l p , i u r=bb(i)/dd(i-l) cc(i) = cc(i) - r * c c ( i - l ) enddo Back substitution cc(iu)=cc(iu)/dd(iu) do i = l p , i u j = iu - i + i l cc(j)=(cc(j)-aa(j)*cc(j+l))/dd(j) enddo return end  95  Bibliography  [1] W. Bartok and S. G. Mason. "Particle Motions in Sheared Suspensions V. Rigid Rods and Collision Doublets of Spheres". Journal of Colloid Science, 12:243-262, 1957. [2] J . F . Brady and G. Bossis. "The Rheology of Concentrated suspensions of spheres in simple shear flow by numerical simulation". Journal of Fluid Mechanics, 155:105-129, 1985. [3] F.P. Bretherton. "The Motion of Rigid Particles in a Shear Flow at Low Reynolds Number". Journal of Fluid Mechanics, 14:284-304, 1962. [4] J. M . Burgers. "Second Report on Viscosity and Plasticity". North Holland Publishing Company, Amsterdam, 1938. [5] J . Feng, H. H. Hu, and D. D. Joseph. "Direct Simulation of Initial Value Problems for the Motion of Solid Bodies in a Newtonian Fluid Part 1. Sedimentation". Journal of Fluid Mechanics, 261:95-134, 1994. [6] O. L. Forgacs and S. G. Mason. "Particle Motions in Sheared Suspensions IX. Spin and Deformation of Threadlike Particles". Journal of Colloid Science, 14:457-472, 1959. [7] 0. L. Forgacs and S. G. Mason. "Particle Motions in Sheared Suspensions X. Orbits of Threadlike Particles". Journal of Colloid Science, 14:473-491, 1959. [8] 0. L. Forgacs, A. A. Robertson, and S. G. Mason. "The Hydrodynamic Behaviour of Paper-Making Fibers". Pulp and Paper Magazine of Canada, 59(5):117-128, 1958. [9] H. L. Goldsmith and S. G. Mason. "Rheology, Vol. 4, edited by F. R. Eirich". Academic, New York, 1964. [10] R. W. Gooding. "The Motion of Fibres Near a Screen Slot". Master's thesis, University of British Columbia, 1990. [11] R. W. Gooding and R. J. Kerekes. "The Motion of Fibres Near a Screen Slot". Journal of Pulp and Paper Science, 15(2):J59-J62, March 1989. [12] G. B. Jeffery. "The Motion of Ellipsoidal Particles Immersed in a Viscous Fluid". Proceedings of the Royal Society (London), A102:161-179, 1922. [13] S. G . Mason. "The Motion of Fibers in Flowing Liquids". Pulp and Paper Magazine of Canada, 51(5):93-100, 1950. [14] S. G. Mason. "Fiber Motions and Flocculation". TAPPI Journal, 37(11):494-501, 1954.  96  Bibliography  97  [15] S. G. Mason and R. St J. Manley. "Particle Motions in Sheared Suspensions: Orientations and Interactions of Rigid Rods". Proceeding of the Royal Society (London), A238:117-131, 1956. [16] J. A. Olson. "Fibre Fractionation in Flow Through Narrow Apertures". Research proposal, University of British Columbia, March 1994. [17] T. Rehmat and R. Branion. "Fibre Fractionation in Hydrocyclones". Article To Be Submitted, August 1993. [18] R. S. Seth, D. W. Francis, and C. P. J. Bennington. "The Effect of Mechanical Treatment During Medium Concentration Fluidization on Pulp Properties". APPITA Journal, 46(l):54-60, 1992. [19] G. A. Smook. "Handbook for Pulp and Paper Technologists". TAPPI, CPPA, 1992. [20] C. A. Stover and C. Cohen. "The Motion of Rodlike Particles in the Pressure-Driven Flow Between Two Flat Plates". Rheologica Acta, 29:192-203, 1990. [21] B. J. Trevelyan and S. G. Mason. "Particle Motions in Sheared Suspensions I. Rotations". Journal of Colloid Science, 6:354-367, 1954. [22] S. Yamamoto and T. Matsuoka. "A Method for Dynamic Simulation of Rigid and Flexible Fibers in a Flow Field". Journal of Chemical Physics, 98(l):644-650, 1993. [23] S. Yamamoto and T. Matsuoka. "Viscosity of Dilute Suspensions of Rodlike Particles: A Numerical Simulation Method". Journal of Chemical Physics, 100(4):3317-3324, 1994. [24] S. Yamamoto and T. Matsuoka. "Dynamic Simulation of Fiber Suspensions in Shear Flow". Journal of Chemical Physics, 102(5):2254-2260, 1995.  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items