@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix dc: . @prefix skos: . vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Mechanical Engineering, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Wherrett, Geoffrey"@en ; dcterms:issued "2009-02-12T19:25:43Z"@en, "1996"@en ; vivo:relatedDegree "Master of Applied Science - MASc"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """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 rigid fibres in 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. The fibre was 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 the fibre was 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 solid fibres of aspect ratio 10 and 60 in a Couette flow. The initial position of the fibre was 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 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, kei = E/ɳG. Here E is the elastic modulus of the fibre, ɳ is the absolute viscosity of the fluid and G is the rate of shear. The model performed rigidly when kei > 2 X 10⁶. Some deformation was seen in the fibre above 5 X 10⁵ but the period 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 kei — 5 x 10⁵. 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 kei, the second is a function of the fibre's hollowness, and the third involves the aspect ratio. 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."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/4545?expand=metadata"@en ; dcterms:extent "6216216 bytes"@en ; dc:format "application/pdf"@en ; skos:note "F I B R E M O T I O N I N S H E A R F L O W 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 OF THE REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF A P P L I E D SCIENCE in THE FACULTY OF GRADUATE STUDIES MECHANICAL E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA May 1996 © Geoffrey Wherrett, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholariy purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia Vancouver, Canada Department 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 rigid fibres in 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. The fibre was 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 the fibre was 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 solid fibres of aspect ratio 10 and 60 in a Couette flow. The initial position of the fibre was 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, kei = 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 kei > 2 X 106. Some deformation was seen in the fibre above 5 X 105 but the period 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 kei — 5x 105. 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 kei, the second is a function of the fibre's hollowness, and the third involves the aspect ratio. 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 rea-sonable manner with respect to previous experimental and theoretical results. Several improve-ments 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 xii 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 Scope and Aims of This Thesis 22 2 Numerical Procedure 23 2.1 Model Formulation 24 2.1.1 Fibre Structural Model . 24 2.1.2 Hydrodynamic Force Model 29 2.2 Method 33 2.3 Application . . . .\" 37 iv 2.3.1 Non-dimensional Groups 37 2.3.2 Computational Algorithm 41 3 Results and Discussion 44 3.1 Introduction 44 3.2 Results for the 10 Element Model 48 3.2.1 Rigid Body Motion - Series 1 48 3.3 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 4 Conclusion 75 A Inertial Properties of a Fibre Element 77 B Structural Properties of the Fibre Model 79 C Dimensionless Parameters and Equations of Motion 81 D Velocity Boundary Condition 85 E Angular Correction Routine versus Rep — Series 2 88 F Code listing 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 fibres 19 2.8 Schematic diagram of a continuous fibre and its representation by a series of cylindrical elements 25 2.9 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 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 Rep particle Reynolds Number se element specific gravity Sfw specific gravity of fibre wall 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 aQ radius of outer wall di diameter of inner wall da 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 kei elastic constant ka fibre spring constant kb fibre bending constant I distance between element centres l0 element length (lQ = 2a0), equilibrium distance between element centres mjw 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 rc cylindrical aspect ratio re ellipsoidal aspect ratio T* equivalent ellipsoidal aspect ratio t time t* dimensionless time (tG) t9 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 x ctij angle formed by the normal n^ -St dimensionless time step e strain j] fluid viscosity u> rotational velocity about z axis Q fluid rotational velocity cf> angular position of rigid fibre pe element density Pfw fibre wall density Pf fluid density a stress 6 angular displacement about z axis Oh bending angle between neighbouring elements 0bo equilibrium bending angle (t translational friction constant £ r rotational friction constant 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 Br i t i sh Columbia 's pulp and paper industry do to improve its posi t ion i n the wor ld wide market? The solution that has been put forward by many can be summarized by the te rm value added. Th is approach suggests replacing product ion of lower value products (ie. newsprint) w i th higher value products (ie. coated paper). Th is approach generally means moving into a narrower market where a new volume of product could lower prices and reduce profit margins. To succeed in the new market a product must have a unique appeal which can be marketed or be manufactured wi th a lower cost technology than compet ing products. B r i t i sh Co lumb ia has an advantage over many areas of the wor ld. The coniferous forests which cover much of the province produce pulp fibres which are long, slender, and flexible. These properties allow the fibres to conform and bond 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 mixed w i th less desirable fibre in the tree. A technology which can separate fibres would be of great benefit in creating unique paper products. The matur ing of the recycl ing industry has also resulted in changes to the fibre supply. A larger propor t ion of fibres have already been through the various manufactur ing processes at least once. A technology which can separate well processed fibres f rom virg in fibres could lower refining costs. The basis for modern pulp and paper manufactur ing has existed for over 100 years. Dur ing this t ime, efforts have focussed on improving exist ing processes through analysis and experi-mentat ion in bo th academia and industry. However there is considerable diff iculty i n model l ing the complex physical processes involved. In recent years, the advent of h igh speed computers has allowed many processes to be modelled numerically. Th is thesis is a first step in modell ing the mot ion of a fibre in a moving f lu id. Us ing the experience gained f rom this model , it may be possible to use this too l along w i th f luid mot ion models to design devices which separate fibres wi th different properties. Th is process is known as fibre f ract ionat ion. Chapter 1. Introduction 3 1.2 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 pulp fibres in 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 (summer-wood) 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/cm3 _ 1 total volume density of cell wall 1.5 g/cm3 4 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 Range Between Sites 0.44 - 0.47 One Site 0.42 - 0.48 Single Tree 0.36 - 0.50 Growth Ring 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 Cel l Wal 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 8 1.3 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 ma-chines. 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 minimum fibre length 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. Earlywood fibres which 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, many fibres do not,have these properties and must be refined to change their phys-ical 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 neighbour-ing 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. The fibres are 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 chopped fibres (known as fines). Reducing the amount of fines increases 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 allow fibres to 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 perfora-tions 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 13 Figure 1.4: Schematic diagram of a hydrocyclone. Taken from Fig. 9-30 in [19] Chapter 1. Introduction 14 1.4 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 individual fibres in 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(re + r:1) (1.1) where T is the period of rotation, G is the shear rate, and re is the ratio of the length of the ellipsoid to its maximum diameter. For large values of the ellipsoidal axis ratio, re, Jeffery's equation simplifies to a linear fuction, TG « 2irre (1.2) The rate of rotation was formulated as, where

= ian _ 1 [r e tan(Gt *e )] (1.4) T e + 1 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=P-^, and = —45°. By considering Euler's classical equation for the shape of a solid rod suffering small defor-mations under compressive forces, the critical value for the first mode of buckling was found, E 2rt W \" * \" l<2rc) - 1.75 ( L 8 j 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: rc and E/Grj, the ratio of the fibre stiffness to the fluid shear. (E/Gr])^ is a function of the fibre geometry (T\"c) and is independent of the size of the particle. 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 of fibre reinforced 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 in fibre stucture 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 22 1.5 Scope and Aims 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 of fibre motion 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 the fibre are 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. Numerical Procedure 24 2.1 Model Formulation A model has been developed to describe a fibre and its dynamic behaviour in a two dimensional shearing flow f ield. In order to provide reasonable values for its mechanical propert ies, a single fibre is considered to be a cyl indr ical rod . The fibre is discretized into a series of elements of equal length and diameter. Each element can stretch and bend w i th respect to its neighbours by changing its bond distance and angle. The continuity of a real fibre is preserved i n the model by mainta in ing a contact point at the interface of elements. The mot ion of the fibre as a whole is determined by solving the translat ional and rotat ional equations of ind iv idua l elements under the hydrodynamic force and torque exerted by the relative mot ion of the f lu id. The fibre is t racked for a specified t ime period by a Lagrangian scheme, the mot ion over a longer t ime being composed of a series of steps for consecutive t ime increments. 2.1.1 Fibre Structural Model To represent i ts st ructural properties, the fibre is made up of N cy l indr ical elements of outer radius a0 and length 2aD. In this manner, the model forms a cyl indr ical rod of length 2aQN, diameter 2aQ and aspect rat io N. The model represents the hollowness of pulp fibres through the use of an inner radius, a{. The thickness of the fibre walls are set by the rat io of ai and aa. The walls are assumed to have a constant density. The interior space of the fibre is assumed to be fi l led w i th water. Th is por t ion of the element adds mass to the model based on the density of water but has no effect on the structural properties. Chapter 2. Numerical Procedure 25 Figure 2.8: Schematic d iagram of a continuous fibre and its representation by a series of cy l in -dr ica l elements The elastic properties of the fibre model are represented in the model through the use of spring and bending constants. The product of these constants and the bond distance and angle determine the magnitude of the spring force and bending torque respectively. This is better understood i f one considers a pair of connected elements i and j. The spring force exerted between elements depends on the bond distance I between the centres of the elements as shown i n F i g . 2.9. If the elements are stretched from their equi l ibr ium distance (taken as lQ = 2a D ) , an equal and opposite force Fs is exerted between the elements drawing them back to the equi l ib r ium distance: Fs = -ks (I - y (2.9) where ks is the spring force constant. Chapter 2. Numerical Procedure 26 I J Figure 2.9: Configurat ion of bonded elements when stretched. The spring force Fs is a linear function of I, the distance between the centres of the elements. T h i s force acts along the normal , a unit vector point ing from the centre of element i to the centre of element j , (2.10) The magnitude of the spring constant is based on the linear deformation of a cy l indr ica l rod . It is evaluated by using the Young's modulus E and the cross sectional area of the model . The full development of the spring force equation and constant is found in Append ix B , where i t is shown that k. = *-Ea0 (1 - g ) (2.11) where is the radius of the hollow inner por t ion of the fibre and aa is the outer radius. Chapter 2. Numerical Procedure 2 7 Figure 2.10: Conf igurat ion of bonded elements when bent. The bending torque Tb is a l inear funct ion of Ob, the bond angle formed between the elements. A s imi lar approach is taken in calculat ing the bending between elements. The bending torque is a funct ion of the bond angle Ob between the elements. A bond angle is formed when neighbouring elements rotate relative to each other away f rom the equi l ib r ium angle Obo (taken as zero in the case of straight fibres). A n equal and opposite torque Tb is then exerted on each element, Tb = -h(0b - Obo) (2.12) where fc& is the bending torque constant. Th i s torque acts about a vector which is perpendicular to . The magni tude of the bending constant is based on the bending deformat ion of a cy l indr ica l rod . It is evaluated by using the Chapter 2. Numerical Procedure 28 Young 's modulus E and the second moment of area of the model about about the neutra l axis. The fu l l development of the bending torque equation and constant is found in Append i x B , where it is shown that : kb = \\Eal (1 - (2.13) 8 a% A model constructed in this manner allows the f lexibi l i ty of the fibre to be changed through the bond parameters ks and k\\,. The overall f lexibi l i ty may be adjusted by varying the stiffness, E, or the inner and outer rad i i , ai and aa. Local ized weaknesses such as kinks may be model led by adjust ing the bond parameters at ind iv idual locations along the fibre length. Chapter 2. Numerical Procedure 29 2.1.2 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. The flow field is 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. The flow field has 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) W(r) = 0 (2.17) 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 ' U = Gy i i t = 0 t>0 • / Figure 2.11: T w o dimensional shear flow. The fluid velocity in the x d irect ion, U is a l inear funct ion of the shear rate G and the y posi t ion. W h e n a fibre is suspended in a shear flow, a hydrodynamic force and torque are exerted on it by the fluid. Th is phenomena is modelled through the choice of suitable elements. The hydrodynamic model involves several important assumptions and simpli f icat ions. A t the largest scale, the fibre is considered to exist under the condit ion of an infinitely di lute system. In other words, a single fibre is being modelled in a flow field unaffected by walls, or other fibres. A \"free dra in ing\" approximat ion is applied to the hydrodynamic force act ing on each ind iv idua l element. Th is means that the force act ing on an element is considered to be solely a funct ion of the bulk fluid mot ion at the locat ion of the centre of the element. Changes in the flow field induced by the mot ion 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 Fh and torque Th 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 aQ. The translational friction force is then Fh = -(t[u-U(y)} (2.20) where (t = 67r77a 0 is the translational friction constant. Similarly, the angular friction torque is Th = -Cr [u, - n(r)] (2.21) where £ r = 8irrjag is the rotational friction constant 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. Chapter 2. Numerical Procedure 32 Hydrodynamic Model Structural Model Figure 2.12: Element geometries: spherical and cylindrical. The outer radius a0 defines the geometry in both the hydrodynamic and structural models. The inner radius defines the fibre wall thickness in the structural model. A model for fibre motion 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 2aa. The hydrodynamic modelling is based on free draining spheres of radius aQ. These geometries are shown in Fig. 2.12. The hydrodynamic forces are linearly 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 33 2.2 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, F8 ensures continuity in the axial direction. However, no analagous 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, where me and Ie are the mass and mass moment of inertia of each element respectively. The 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, (2.22) (2.23) Chapter 2. Numerical Procedure 34 Figure 2.13: Hydrodynamic forces exerted on the elements cause them to move apart as seen in the two outer elements. Tangent ial forces act at the interface to prevent this. dri dOi drj dt dt J dt d6j (2.24) The contact point is assumed to occur at a distance a0 f rom the centres of both elements along n-ij. The non-slip condit ion is applied so that the velocities perpendicular to n{j are equal at the contact point. The magnitude of / * is not directly calculated in this procedure. Instead, a set of simul-taneous 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 Eqs. 2.22 and 2.23 to find the accelerations while neglecting the f* terms. New velocities at t + 6t are found by adding the contr ibut ions of the accelerations using a finite difference technique. A new set of velocity equations including the unknown velocity ul due to jf* are formed. It is these new velocities which are subst i tuted into E q . 2.24. The terms involving v} are then isolated so that they may be solved. For a fibre made up of N elements, there wi l l be (N — 1) interfaces and (N — 1) simultaneous equations which must be solved to determine the velocity contr ibut ion of jf*. Chapter 2. Numerical Procedure 35 ei + Qj = 2aij) tf = 0 Qi + Qj<2aiJ! j, hg formed by the fibre with the y axis was plotted against time and the period TG was determined. The period of flexible fibres was determined by inspecting the fibre plots. The degree of deformation seen in flexible fibres was found by plotting the orbits formed by the fibre tips as they rotate. Chapter 3. Results and Discussion 48 3.2 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, Rep — 0.1 This is the same value for Rep 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 kei was set at 104. This results in a bending constant kb of 3930. This value was well within the range of rigid fibre motion cited in the previous work [22]. Element Radius, a0 100 nm Shear Rate, G 10 s'1 Fibre Wall Density, pfw 1000 kg/m3 Elastic Modulus, E 100 Pa Table 3.5: These values were kept constant in all the trials run in Series 1 The trials were executed for m = 5 X 106 steps using a dimensionless time step of St = 10 - 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 Trial 1(a) JefFery's Equation t* At* t* 90° 12.5 12.5 12.39 180° 23.8 11.3 24.78 270° 36.1 12.3 37.18 360° 47.4 11.3 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 > AiJ 8 0 ) This disagrees with JefFery's equation which predicts equal time over any quadrant of rotation. A delay of t* = 0.2 is apparent in Atg 0. This may be attributed to the time required by the fibre to come up to speed from its initial angular velocity of zero. The Chapter 3. Results and Discussion 50 360 r i / / / 315 '-i i 1 i 1 i J 270 225 180 (deg) 135 90 / / i i i i J i i 1 i J i J i f Trial 1(a) Jeffery's equation 45 , . i . . . , i . . . . i . . . . i . . . . i 0 10 20 30 40 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 \"start-up\" 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 Trials 1(a) Trial 1(b) JefFery's Equation t* Ai* i* Ai* i* 90° 12.5 12.5 11.2 11.2 12.39 180° 23.8 11.3 22.8 11.6 24.78 270° 36.1 12.3 34.1 11.3 37.18 360° 47.4 11.3 45.6 11.5 49.57 Table 3.7: Results for Trial 1(b). i* refers to the time to complete the specified rotation. Ai* 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 Rep = pea0G/rj is the same for both spherical and cylindrical geometries. The change to the cylindrical geometry is reflected in the value of the constants CtT and CTOt-Ctr is proportional to the element's mass which increases by 50% for the cylindrical geometry. Crot is proportional to the mass moment of inertia of an element and increases by 218% for the cylindrical geometry. The assumption of Stokes flow in this model implies that the fibre centroid 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 a fibre length 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). The fibre also 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 a fibrelength behind 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 CTOt-Chapter 3. Results and Discussion 53 1 \\ o o 2 4 6 8 10 12 Fibrelengths Figure 3.17: Trial 1(c): A rigid 10 element model translating and rotating through approxi-mately 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 Rep. The constants CtT and Crot are proportional to Rep for solid elements. The 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 Rep ranging from 0.001 to 1.0. The results showed that the period of rotation increased dramatically above Rep = 0.1 when the correction routine was in use. The period for Rep = 1.0 was nearly 150. The trials run without the correction routine remained nearly constant over Rep with the period for Rep = 1.0 dropping slightly to 44.6. It was concluded from Series 2 that subsequent trials would be performed without the correction routine due to the problems seen at Rep 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 Ctr-For Rep = 0.1, the lag was found to be 0.0165 of a fibre length. For Rep = 1, the lag increased 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 54 3.3 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, Rep in Series 3 and 4 was kept constant at 0.15. The parameters used in the formulation of the number are listed in Table 3.8. Outer Radius, a0 10 microns Shear Rate, G 1000 s-1 Fibre Wall Density, pfw 1500 kg/m6 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. Pro-cessing 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 106 time steps while the 60 element model 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 55 3.3.1 Bending Deformation — Series 3 Series 3 consisted of five trials where the elastic constant, kei was varied from 2 X 106 to 104 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 kel = E/rjG Kb *: St 3(a) 2 X 106 7.85 x 105 3.14 x 106 2 x lO\"6 3(b) 106 3.93 x 105 1.57 x 106 5 x 1(T6 3(c) 5 x 105 1.96 x 105 7.85 x 105 1 x 1(T5 3(d) 105 3.93 x 104 1.57 x 105 5 x 1(T5 3(e) 104 3.93 x 103 1.57 x 104 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, kei was increased to 2 X 106. The deflection of the fibre from its initial staight position was limited to 0.001 fibrelengths as shown in Fig. 3.19. Angle Trial 3(a) JefFery's Equation t* At* t* 90° 65.6 65.6 62.56 180° 133.4 67.8 125.11 270° 198.8 65.4 187.68 360° 266.6 67.8 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 pe-riod 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 57 4> (deg) 450 r 405 360 315 270 r 225 -180 135 r 90 F-45 0 J i_ J I I 1 L _ Trial 3(a) Jeffery's equation • i i i i i i i i i 50 100 150 200 250 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 kei = 106 and 5 X 105 respectively. 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 b)t*= 133.2 c ) t /= 132.2 66.6745 66.6750 66.6755 66.4 66.5 66.6 66.7 66.8 65.9 66.0 66.1 66.2 Figure 3.19: Fibre Deflection at cp = 180° for the three stiffest trials: (a) kei = 2 X 106, (b) kei — 106, and (c) kei = 5 X 105. The axes are measured in fibre lengths. The axes in a) are not to scale in order to magnify the deflection. The deflections of the fibres in Trials 3(a) to 3(c) at

G-n,crit fn(2rc)-1.75 where rc is the aspect ratio of the cylindrical fibre and a-a refers to the value above which the 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)rigid, 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 kei = E/r/G ranging from 50 to 10000. This allowed the flexible motion to be compared 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. kel B (Eq. 3.42) B (Eq. 3.43) T G % of TGrigid 50 6.25 X IO\"3 3.11 X IO\"3 39.7 87.1 100 1.25 X IO\"2 6.23 X IO\"3 40.4 88.6 200 2.50 x IO\"2 1.25 x IO\"2 41.0 89.9 500 6.25 x 10~2 3.11 x I O - 2 43.3 95.0 1000 0.125 6.23 x 10~2 45.0 98.7 2000 0.250 0.125 45.5 99.8 4000 0.500 0.250 45.6 100 10000 1.25 0.623 45.6 100 Table 3.12: Results for a flexible 10 element fibre and the two calculated bending ratios kel B (Eq. 3.42) B (Eq. 3.43) T G % Of TGrigid 104 5.79 X IO\"3 1.17 X IO\"3 148 55.5 105 5.79 X IO\"2 1.17 x IO\"2 210 78.8 5 x 105 0.289 5.86 X IO\"2 264.2 99.1 1 x 106 0.579 0.117 266.2 99.8 2 x 106 1.16 0.234 266.6 100 Table 3.13: Results from Series 3 and the two calculated bending ratios TG/(TG)Tigid versus B is plotted in Fig. 3.25 for the two different forms of the bending number B. Figure 3.25: TG/(TG)rigid 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 Lower Limit Upper Limit Diameter (microns) 30 50 Wall Thickness (microns) 2 10 Moment of Area, Ics (m4) 1.73 x I O - 2 0 2.67 x 10~19 Length (mm) < 2 4 Aspect Ratio 40 100 Stiffness, EIC3 (Nm4) 1 X IO\" 1 2 (Kraft) 100 x I O - 1 2 (TMP) Table 3.14: Variation of fibre parameters for softwood pulp in British Columbia. Kraft and TMP (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 c = 25 microns, G = 1000 s - 1 , relative velocity = 23.6). The actual relative velocity is, Urei = 23.6(25 x 10_6m)(1000 s _ 1 ) (3.44) = 0.59m/s (3.45) The Reynolds number is then calculated using standard values for the density and viscosity of water, Re = (3.46) _ (0.59)(25xl0-6)(1000) = 14.75 (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 aD = 25 microns, G = 1000 s _ 1 , and EI = 1 X 10 _ 1 2JVm 4 . For an aspect ratio of 100, as calculated by Eq. 3.43, B is equal to 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 73 3.5 Recomendations for Improvements to the Model Wall Thickness The model as presently formulated uses a cylinder as the basis for the elastic and inert ia l properties of the fibre. The varying wal l thickness of real wood fibres is represented in the model through the rat io of inner to outer radius in the fibre, a* = ai/a0. The ax ia l stiffness, EACS and the bending stiffness, EICS of the fibre model are a funct ion of the Young 's modulus and the cross section of the fibre wal l calculated f rom a*. The present formulat ion assumes that the fibre wal l strains linearly. In reality, pulp fibres behave differently whether stretched or bent due to their complex geometry of various layers and fibril angles. It was also found in the trials that varying the axia l stiffness of a fibre had no significant effect on its dynamic behaviour for the range of stiffnesses tested. The available data for pulp fibres usually specifies the bending stiffness, EICS, the outer diameter, aa, and the density, pp. These considerations suggest some simplif ications that could be made i n the model . F i r s t , the calculat ion of ax ia l forces in the fibre could be neglected. They could be replaced by a non-slip force, similar to that which acted tangential ly between the elements in the present model . The new non-sl ip force would have components in both the axia l and tangent ial directions to ensure that neighbouring elements do not break continuity. The wal l thickness as presently defined makes the formulat ion of several parameters more compl icated than necessary. It is recommended that the properties of the model be brough into l ine w i th the available data. Th is would remove the rat io a* f rom the formulat ion of the fibre model 's properit ies. 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, kei, there existed a critical time step, St above which the instability occurred. The relation between kel and St 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 t r . For the values of C 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 kei > 2 X 106. For 1 X 106 > kei > 5 X 105, the period was reduced by less than 1%. For kei < 5 X 105, the period reduced at a rate proportional to log(Cei). The fibre model deformed in a manner 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 EIC3 and pe The present formulation of these parameters makes them a function of a*, the wall thick-ness 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, me m-element = \"Vibre wall + \" V l u i d me = m,fw + mf = AcJ0pfw + Afl0pf = 7r ( a 2 - a?) (2a0) pfw + 7ra 2 (2a 0 ) 2 = 27ra 3 (1 - - | - ) ^ / u , + 27ra 0 a 2 /?/ al 2 = 27ra 3 [ pfw - -± (pfw - pf) } ao Density, pe and Specific Gravity, Se Pe = Pfw ~ \\ (Pfw ~ Pf) Pf C _ Pfw Jfw -Pf Pfw a,2 ,pfw _ 1 x CX2 -S/u; - (Sfw - 1) a ? . a? «2' 77 Appendix A. Inertial Properties of a Fibre Element Mass Moment of Inertia, Ie lelement — Izz Ie. = I fibre, wall + Ifluid = Ifw + If Ifw » outer ^ ( 3 ^ + 4Z 2) - ^ ( 3 a * + 4Z 2) ^ { \\ < lo Pfw) (3d 2 + All) ~ ^ {\\dj la pfw) (3d 2 + 412) (2ira30 pfw) (3a 2 + 4a 2 ) - ^ (27ra, 2a 0 pfw) (3aj + 4a 2 ) 1_ 12 rrif - h r rrii 12 [27ra0 ( a 2 - a?) ( 3a 2 + 3 a 2 + 4 a 2 ) (7a 2 + 3a 2 ) = ^ ( 3 a t 2 + 4a 2 ) Ie = ^ (3a t 2 + 7a 2 ) + ^ (3a 2 + 4a 2 ) Appendix B Structural Properties of the Fibre Model Spring Force Equation (Ee) Aca EA — EA<° Al 2a0 = kaAl For a hollow cylindrical rod, the cross sectional area is: = *(«2-«?) = ™ o ( i - 4 ) The spring constant, ks is then: 7T „ a? *. = -Ea0(l-^) 2 . v al 79 Appendix B. Structural Properties of the Fibre Model 80 Bending Torque Equation Tb = EIr = EL EIC, 2a0 dd dx AO lo A6 = kb AO For a hollow cylindrical rod, the second moment of area, Ics about the neutral axis is: The bending constant, kb is then: kb = ^Ea\\ (1 -Appendix C Dimensionless Parameters and Equations of Motion Dimensional Parameters Unit Parameter Value Time 1/G s Length aa m Velocity aDG m/s Acceleration (translational) aaG2 m/s2 Acceleration (rotational) G2 Force alvG kg m/s2 Torque a0nG kg m2/s2 81 Appendix C. Dimensionless Parameters and Equations of Motion Dimensionless Equations of Motion Translation x Direction du me — dt Am* du\" Similarily for the T/ direction Rotation (0) (Th* + E T B* + E /*• x « « + E r _ e i + 1 [«i + (/' X.njfc - /i-i X nW-iJ/Ctron,] Xny{ ~ [vi ~ (fi X n x i ~ fi-1 X n x i - l ) / C t r a n 3 ] XTIX; — cy? -»-= [ui+i + (/ t ' + 1 x n y i + 1 - ft x nyj/Ctrans] xnyi - [vi+1 - (fl+1 X nxi+1 + ft x n x i ) / C 4 r a n s ] Xnxi + [ui+i-{fi+fU)icTOt\\ Collect Coefficients of /* fi-i [(-nyi-i X nyi - nxi-i x nxi)/C'trans + l / C r o t ] + /' [{2nVi + 2nx?)/C t r a n s + 2 /C r o t ] + fi+l [{-nyi X - nXi X nXi+1)/Ctran3 + 1/Crot] + = (iti+i - Ui) X - - Vi) X nii + + a;;) / i - i X nyi - nxi_i X n x i ) / C t T a n a + l / C r o t ] + /• [2 / C t r a n s + 2 /C r o t ] + [{-nyi X nyi+1 - nxi X nxi+1)/Ctrans + l / C r o t ] + = {ui+i - Ui) X nyi - (vi+1 - Vi) X nxi + + Ui) fi-i X bbi + ft x ddi + x aai = b{ Appendix D. Velocity Boundary Condition 87 A Tridiagonal Array is formed, dd aa bb dd 0 bb 0 0 0 • • 0 aa 0 dd aa 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 Rep - Series 2 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, Rep. All parameters except pfw were kept constant from Trial 1(c). Rep was varied by adjusting from 10 kg/m3 to 10000 kg/m3. As in Series 1, the behaviour of the model was determined by the time required to complete 90°, 180°, 270° , and 360° of rotation. The fibre was 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 Rep ranging from 0.001 to 1. Three trials were completed without any angular correction for Rep = 0.001, 0.1, and 1. Trial 2(a) 2(b) 2(c) 2(d) 2(e) 2(f) Trial 2(g) 2(h) 2(i) Rep I O - 3 10~2 IO\"1 0.15 0.30 1.0 Rep 10\" 3 i o - 1 1.0 '90 11.5 11.6 14.5 17.1 25.3 67.0 '90 11.4 11.3 10.0 Ai* 11.4 11.4 11.4 11.5 12.6 16.8 Ai* 11.4 11.5 12.5 l180 22.9 23.0 25.9 28.6 37.9 83.8 4* l180 22.8 22.8 22.5 Ai* 11.5 11.5 13.6 15.0 19.5 57.8 Ai* 11.4 11.2 9.6 +* '270 34.4 34.5 39.5 43.6 57.4 141.6 '270 34.2 34.0 32.1 Ai* 11.4 11.4 11.3 11.5 12.6 16.8 Ai* 11.4 11.6 12.5 ••360 45.8 45.9 50.8 55.1 70.0 158.4 l360 45.6 45.6 44.6 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, Rep = 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 Rep - Series 2 89 160 r Figure E.26: Plot of dimensionless period, T G versus Rep for trials in Series 2 nearly identical results. As Rep was increased above 0.1, several trends were visible from these results. Trials 2(c) to 2(f), with the correction routine intact, exhibited a rapid increase in the time to complete the initial quadrant of rotation, i*,0, and in rotating from the vertical to horizontal, At*,70. The 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 Rep increased. These factors combined to cause the period to increase rapidly for trials 2(a) to 2 (f) as Rep 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(£*.60 - i j 8 0 ) . This ensured that the rotational delay due to the inital conditions did not have an effect on the value. Appendix F Code listing c************ 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 Assign i n i t i a l constants c mml=m-l c ao = 100.0d-6 radrat = 0.0d+0 g = 9.81d+0 rhop = 1.0d+3 rhof = 1.0d+3 E = 100.0d+0 eta = 1.0d-3 shrate = 10.0d+0 pi = 3.14159d+0 c area = pi*(ao**2)* (1.0d+0-radrat**2) arin = pi/4.0d+0*(ao**4)* (1.0d+0-radrat**4) amf = 2.0d+0*pi*(ao**3)*(radrat**2)*rhof amp = 2.0d+0*pi*(ao**3)* (1.0d+0-radrat**2)*rhop amc = amp + amf amin = 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) c mm is the number of iterations 90 Appendix F. Code listing 91 t = O.Od+0 delta = 1.0d-05 mm = 5.Oe+6 c c I n i t i a l Conditions for Elements 1 - m c x(l)= 0.0d+0 y(l)= 1.0d+0 dx=2.Od+0 dy=2.Od+0 c do i=2, m x(i)=x(i-l) + 0.0d+0 * dx y(i)=y(i-i) + 1.0d+0 * dy enddo c 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 Initialize variables to 0 c 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 c c Initialize spring and bending constants c do i=0,m cks(i) = cksl ckb(i) = ckbl enddo c c Specify kink location c c ckb(il) = 0.2d+0*ckbl c c Open Data Files c open (unit=31,file='main.dat',status='unknown') c open (unit=41,file='fs.dat',status='unknown') c open (unit=5i,file='torb.dat',status='unknown') c open (unit=61,file='ft.dat',status='unknown') c open (unit=7i,file='torg.dat',status='unknown') c Appendix F. Code listing c Input values from previous t r i a l (if needed) c open(unit=8,file='flag.in',status='old',form='formatted') read(8,*) iflag 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 c c Record I n i t i a l Position of Elements c c write(31,*) 'zone' c c do i=l,m c write(31,ll) ao*x(i),ao*y(i),t/shrate c enddo 11 format(dl6.8,2x,dl6.8,2x,dl0.5) c c c Begin time do loop c do 500 ii=l,mm c t= t + delta c c Calculate normals, spring & bending forces c do i=l,mml rr = 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))/ rr cny(i) = (y(i+l)-y(i))/ rr fs(i) = - cks(i)*(rr-2.0d+0) torbz(i) = ckb(i) * (angle(i+l)-angle(i)) enddo c c Calculate accelerations and update velocities c 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 c c store velocities from start of time step and calculate new values c 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 93 enddo c c output bending and spring values c c if(mod(ii,mm/1000).eq.O) then c write(41,*) 'zone ' c write(51,*) 'zone ' c do i=l,m c fh(i) = -6.0d+0*pi* (u(i)- y(i)) c spr(i) = -cnx(i)*fs(i)+cnx(i-l)*fsXi-1) c th(i) = -8.0d+0*pi* (rotvel(i)+0.5d+0) c ben(i) = torbz(i)-torbz(i-l) c write(41,ll) spr(i), fh ( i ) , t c write(51,ll) ben(i), th(i), t c enddo c endif c c find Non-Slip force c c define tridiagonal terms c 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 c dd(l)=-2.0d+0*cnx(l)*cnx(l)/ctrans 1 -2.Od+0*cny(1)*cny(1)/ctrans-2.Od+O/crot aa(l)= cnx(l)*cnx(2)/ctrans 1 + cny(2)*cny(l)/ctrans - 1.Od+O/crot c do i=2,mml-l dd(i) =-2.0d+0*cny(i)*cny(i)/ctrans 1 - 2.Od+0*cnx(i)*cnx(i)/ctrans-2.Od+O/crot bb(i) = cny(i-l)*cny(i)/ctrans 1 + cnx(i)*cnx(i-l)/ctrans - 1.Od+O/crot aa(i) = cny(i+l)*cny(i)/ctrans 1 + cnx(i)*cnx(i+l)/ctrans - 1.Od+O/crot enddo c i=mml dd(i) =-2.0d+0*cny(i)*cny(i)/ctrans 1 - 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 - 1.Od+O/crot c c c a l l subroutine tridiagonal c ca l l tridiagonal(1,mml,bb,dd,aa,b) c c define f t from solved non-slip condition c do i=l,mml ft(i)=b(i) enddo c c update velocities from f t c Appendix F. Code listing 94 do i=l,m u(i)=u(i) + cny(i)*ft(i)/ctrans - cny(i-l)*ft(i-l)/ctrans v(i)=v(i) - cnx(i)*ft(i)/ctrans + cnx(i-l)*ft(i-l)/ctrans rotvel(i)=rotvel(i) - ft(i)/crot - f t ( i - l ) / c r o t enddo c c find another tg from geometric conditions c c do i=l,mml c alold=alpha(i) c if(abs(x(i+l)-x(i)).lt.l.0d-10) then c alpha(i)=-dble(npi(i))*pi C ©Xs 6 c alpha(i)=-(pi/2.0d+0-datan((y(i+l)-y(i))/(x(i+l)-x(i)))) c + -dble(npi(i))*pi c endif c c if(abs(alpha(i)-alold).gt.pi*2.0d+0/3.0d+0) then c alpha(i)=alpha(i)-pi c npi(i)=npi(i)+l c endif c tg(i)=ckbl/crot*(angle(i+l)+angle(i)-2.0d+0*alpha(i)) c enddo c c do i=l,m c rotvel(i)=rotvel(i)-tg(i)*delta-tg(i-l)*delta c enddo c c output f t & tg values c c if(mod(ii,mm/1000).eq.O) then c write(61,*) 'zone ' c write(71,*) 'zone ' c do i=l,m c spr(i) = cny(i)*ft(i) - cn y ( i - l ) * f t ( i - l ) c ben(i) = -tg(i)*delta-tg(i-l)*delta c write(61,llj spr(i), t c write(71,ll) ben(i), t c enddo c endif c c Calculate new displacements and rotation c 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 c c Record displacements and time c 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 c Appendix F. Code listing 95 500 cont inue c c Output fin a l values c open(unit=9,f ile='f inal.dat',form='formatted') write(9,*) x,y,angle,alpha,npi,u,v,rotvel,t close(9) c c Close Data Files c close(31) c close(41) c close(51) c close(61) c close(71) c stop end c subroutine tridiagonal(il, iu, bb, dd, aa, cc) double precision aa, bb, cc, dd dimension aa(iu), bb(iu), dd(iu), cc(iu) lp=il+l do i = lp, iu r=bb(i)/dd(i-l) cc(i) = cc(i) - r*cc(i-l) enddo c c Back substitution c cc(iu)=cc(iu)/dd(iu) do i = lp, iu j = iu - i + i l cc(j)=(cc(j)-aa(j)*cc(j+l))/dd(j) enddo return end c 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 Num-ber\". 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\". Pro-ceedings 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 Sub-mitted, August 1993. [18] R. S. Seth, D. W. Francis, and C. P. J. Bennington. \"The Effect of Mechanical Treat-ment 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. "@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "1996-11"@en ; edm:isShownAt "10.14288/1.0080821"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Mechanical Engineering"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Fibre motion in shear flow"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/4545"@en .