UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Fibre motion in shear flow 1996

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata


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

Full Text

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 <j) = 180° for Trials 3(a) to 3 (c) 58 3.20 Fibre Motion for Trials 3(a) to 3 (c) 59 3.21 Fibre Motion for Flexible Trials 61 vi 3.22 Stiffness vs Reynolds Number 62 3.23 Fibre Orbits for Series 3 63 3.24 Fibre Motion for Trials 4(a), 4(b), and 4(c) 65 3.25 a) TG/(TG)rigid versus B (Eq. 3.42), b) TG/{TG)rigid versus B (Eq. 3.43) . . . 69 E.26 Rigid 10 Element Fibre, T G vs Rep 89 vii List of Tables 1.1 Variation of Dry Relative Density for Douglas Fir trees in British Columbia . . . 5 2.2 Fibre and Flow Parameters 39 2.3 Dimensionless Groups 40 3.4 Constants and Parameters 46 3.5 Constant Values used In Series 1 48 3.6 Results for Trial 1(a) 49 3.7 Results for Trial 1(b) 51 3.8 Constant Values used In Series 3 and 4 54 3.9 Series 3 Elastic Parameters 55 3.10 Results for Trial 3(a) 56 3.11 Constant Values used In Series 4 64 3.12 Results for a flexible 10 element fibre and the two calculated bending ratios . . . 68 3.13 Results from Series 3 and the two calculated bending ratios 68 3.14 Variation of fibre parameters for softwood pulp in British Columbia 71 E.15 Results from Series 2 . 88 viii LIST OF S Y M B O L S A.ca cross sectional area of fibre wall Af cross sectional area of fluid portion of element B fibre bending number E Young's modulus CTOt rotational constant CtT translational constant ph hydrodynamic force ps spring force G shear rate I> Ics second moment of area of fibre cross-section h element mass moment of inertia If mass moment of inertia of fluid portion of element Ifw mass moment of inertia of fibre wall portion of element N number of elements Re Reynolds Number 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 <p is the orientation angle of the ellipsoid as it rotates in a clockwise direction from the y axis as shown in Fig. 1.5. Chapter 1. Introduction 15 Figure 1.5: Orientation angle of Rotating Particles When integrated with time, an equation for the angular displacement is formed, <t> = 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-^<l (1.5) that is, the product of the relative velocity of the ellipsoid by its length, ul must be smaller than the kinematic viscosity of the fluid, v = r\\p. Chapter 1. Introduction 16 70 60 50 * r e 40 30 20 10 , , i , , i i i • , , , 1 , , , ' , 1 . . . . 1 . . 20 40 60 80 100 Cylindrical Aspect Ratio, rc 120 Figure 1.6: Equivalent Ellipsoidal Aspect Ratio, r* versus Cylindrical Aspect Ratio, rc. This figure was adapted from data in [22], [9], and [14]. Bretherton [3] subsequently extended JefFery's work to describe the motion of any axisym- metric particle by using the concept of an equivalent ellipsoidal axis ratio. In 1950, Mason [13] stated that pulp fibres may be regarded as cylindrical particles. Trevelyan and Mason [21] showed that JefFery's equation could describe the motion of a cylindrical particle by substituting an equivalent ellipsoidal axis ratio, r* for re in JefFery's equation. The value of r* was calculated by comparing measured periods of rotation from experiments with cylindri- cal rods of aspect rc with the periods obtained from JefFery's equation. Fig. 1.6 exhibits the variation of rc with r*. The difference in shape between the rods used in the experiments and the cigar shaped ellipsoids upon which the theory is based results in differing rates of rotation. Chapter 1. Introduction 17 The assumption of Stokes flow for the motion of a suspended pulp fibre allows Jeffery's equation to be used to describe its motion. This assumption is supported by the small size of a typical softwood fibre (length = 4 mm). It can also be assumed that the fibre will move at a low velocity relative to the surrounding fluid due to the similarity in densities. The fibre wall has a specific gravity of approximately 1.5. Depending on the thickness of the fibre wall the fibre will have an overall specific gravity « 1.1 when immersed in water. Jeffery's equation assumes that, apart from the disturbance produced in the immediate vicinity of the particle, the fluid motion is steady and varies in space on a scale which is large compared with the dimensions of the particle. The assumption of an isolated fibre in a Newtonian fluid also serves as the basis for the computational model developed in this thesis. For this reason, Jeffery's equation, using r*, will serve as a benchmark for the rigid fibre model. 1.4.2 Flexible Fibre Motion The modelling of flexible fibre motion in a suspension is much more difficult than that of rigid fibre motion. Even a theoretical model of an isolated fibre will be limited by the fact that the fibre shape, and hence forces acting on it, is not constant. Much of the literature has focussed on experiments conducted with pulp and other types of fibres in order to describe the characteristics of flexible fibre motion. A large body of work on the motion of particles in suspension was created by Mason and others after 1950 at the Pulp and Paper Research Institute of Canada (PAPRICAN). Much of it was published as a series of articles in the Journal of Colloid Science [21, 1, 6, 7]. This work has particular application to the present research due to its focus on the theory and classification of pulp fibre motion. Chapter 1. Introduction 18 Mason [14] showed that the motion of a single wood fibre is more complex than the model in Jeffery's equation due to the fibre flexibility. Forgacs, Robertson, and Mason [8] divided the range of motion of fibres into the four "orbit classes" shown in Fig. 1.7. The fibre orbits are described below: • Type I - Rigid Rotation The fibre remains rigid throughout the rotation and its period and rate can be described by Jeffery's equations. • Type II - Springy Rotation The fibre begins unstressed in a horizontal position. As it rotates, a compressive force acts on the fibre due to the fluid. This force reaches a maximum at 45° and the fibre responds by bending. The deformed rotation continues until the fibre approaches 135°. At this point the fibre straightens and continues its rotation. • Type III - Flexible Rotation As the flexibilty of the fibre increases, one or both ends of the fibre bend inward soon after it has rotated from the horizontal position. In a loop turn both ends bend in to create an S-shape. This shape persists in the fibre as it rotates before straightening along the x axis near 180°. A snake turn occurs when the forward end of the fibre bends upward and undulates down the fibre until it is facing backward. This turn is more common than the loop turn in actual shear flows due to asymmetries in the fibre. • Type IV - Complex Orbit When the fibre is sufficiently flexible it will continue to rotate in a snake turn motion without straightening at any point in the rotation. Chapter 1. Introduction 19 1 t I — \ / — / — x II - v ( r - III — s s J - — — < _ c c - — IV O C L _ c Figure 1.7: Typical rotational orbits of fibres in a one dimensional shear field as shown on the top left of the figure. The fibre flexibility increases from Classes I to IV. Adapted from Fig. 14 in [8]. The authors [8] also identified factors in actual pulp fibres which account for their behaviour. In particular, the nonuniformity of fibre stiffness was noted. There are sections along the length of the fibre which are more flexible than others. There often exist localized areas of weakness which are caused by fractures in the cell wall. Chapter 1. Introduction 20 Onset of Buckling The ability of a fibre to resist deforming in a shear flow depends on its stiffness and the flow parameters: viscosity and shear rate. Forgacs and Mason [6] put forward a theory for the deformation of cylindrical particles rotating in a velocity field. Mason had stated in an earlier paper [13] that wood fibres could be regarded as cylindrical in shape. Equations were developed which calculated the critical value of the ratio E/Grj above which the shear-induced axial compression causes the particle to buckle. The theory was based on an approximate theory of the forces acting on a thin rod in laminar shear developed by Burgers [4]. For a rigid particle rotating in a one dimensional shear, Burgers showed that the total axial force, F acting on the central cross section of the fibre is approximately, F = , ,  TGf M (1.6) ln(2rc) - 1.75 V ; (1.7) where I is half the fibre length, rc is the cylindrical aspect ratio (length to diameter ratio), M = sine)) cos<j>, and </» is the angle of the rod measured clockwise from the y axis. The sign of F depends upon whether M is negative or positive, and determines whether the fibre is under compression or tension. The maximum compressive force is reached at d> = —45°. By considering Euler's classical equation for the shape of a solid rod suffering small 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<o Figure 2.14: Schematic Diagram showing nonslip condition. Oi + Oj = 2ctij The use of finite difference solutions to differential equations introduces numerical errors which cause the model to violate precise continuity. In particular, an angular non-slip condtion may be broken when the model bends. The geometric form of the condition is: where Q ^ J is the angle of from the fixed axes. This condition requires that nieghbouring elements " ro l l " wi th respect to each other so that their contact point does not slip. A diagram of this condition can be seen in F i g . 2.14. The use of spheres in the figure implies that the contact point remains at a constant distance aD from the centres of the elements. If longer elements were used in the model the condition would change to reflect the new geometry. Another possibility in the future development of the model would be to use a coordinate system which guarantees continuity of the fibre automatically. Bi + Oj = 2ai3 (2.25) Chapter 2. Numerical Procedure 36 In the present model the correction is applied by adding the contribution of a torque, where kg is a constant that is adjusted according to the degree of correction desired. In practice the value of the bending constant kf, was also used for kg. The rotational acceleration induced by tg results in an additional contribution to the angular velocity of each element over the time step 6t. Average velocities (u, v, and u) over the time step St are calculated from the mean of the velocities at the start and end of the time step. New translational and rotational displacements are found by adding the contribution of the velocities over the time step. These values are recorded for each element and the time step is advanced. New equations of motion are then formed by calculating the hydrodynamic and fibre structural forces. The motion of the whole fibre may be followed by continuing this process over a series of time steps. The computational algorithm followed in the application of the procedure is detailed in the next section. h[(0j + Oi) - 2aij]/Crot (2.26) Chapter 2. Numerical Procedure 37 2.3 Application 2.3.1 Non-dimensional Groups In order to implement the numerical procedure outlined in this chapter, a computer code was written in the Fortran programming language. First the equations were expressed in dimensionless form. It is the fluid forces which are responsible for the motion of the fibre in this model. As such, a20r\G was chosen as the unit of force in the simulation. The units of length and time were selected as aa and 1/G respectively. Other dimensional units used in the simulations are identified in Appendix C as well as the dimensionless equations of motion. These equations are as follows, % = ^ ; (^ + E ^ + £ n (2-27) ^ = ^ ( E ^ + E n (2-28) % = 7^( r' l* + E r 6 * + E / t * X ^ + Er9*) (2-29) where * refers to the dimensionless value. C t 7 . and CTOt are the translational and rotational constants which relate the inertial and hydrodynamic forces in the model. Their formulation is shown below. C t T = 2 7 r ^ ^ (2.30) Crot = ^ (2.31) Chapter 2. Numerical Procedure 38 Ctr is proportional to Rep, the particle Reynolds number defined as Rep =  a °y G . This dimensionless group is unlike the traditional Reynolds number since it is the inertia of the element that is considered in its formulation; not that of the fluid surrounding it. Fluid inertia! forces are neglected because of the assumption of Stokes flow. The element density, pe includes the contribution of the fibre wall and the fluid immobilized inside the walls. The fluid density thus has an effect on the inertia of the element. The term particle Reynolds number will be used in this thesis to emphasize the use of the particle's inertia as opposed to that of the fluid. Crot is related to Rep in a more complicated manner than Ctr because of the form of the element's mass moment of inertia Ie which is, U = (3at2 + lal) + ^ (3at2 + 4a*) (2.32) where rrifw is the mass of the fibre wall and mw is the the mass of the fluid immobilized inside the hollow fibre. As can be seen from the equations for J e , the inner radius a{ of an element has an effect on its inertial properties. It does so by changing the wall thickness of the fibre as well as the amount of immobilized fluid. The wall thickness is represented nondimensionally in the model through the ratio of inner to outer radius in the fibre, a* — ai/aa (2.33) The elastic properties of the fibre are defined by the product EI, where E and I refer to the Young's Modulus and second moment of area respectively. Changing the wall thickness ratio, a*, affects J as shown in Appendix B. A list of the parameters and the forces which they influence is given in Table 2.2. Chapter 2. Numerical Procedure 39 Force Model Parameter Inertia! Fibre Outer Radius Inner Radius Fibre Wall Density Fluid Fluid Density Hydrodynamic Fibre Outer Radius Fluid Shear Rate Viscosity Elastic Fibre Outer Radius Inner Radius Elastic Modulus Table 2.2: The fibre and fluid parameters which define the model. Note that the fluid density which effects the fibre inertia refers to the fluid immobilized inside the fibre walls. Since the fibres are immersed in water it is assumed that they fill with water. A dimensionless group involving E appears in the process of forming the equations. The elastic coefficient, kel = E/Gri (2.34) relates the elastic force to the hydrodynamic force. This group is found in the dimensionless forms of the bending and stretching constants, A£ and A:*. These constants determine the rigidity of the fibre. They are functions of the Young's modulus, and the cross-sectional properties (ie. area, second moment of area). The forms of these constants are: *: = EAC3/a20r,G (2.35) = f(kei,a*2) (2.36) k*b = EIca/atvG (2.37) = / ( W 4 ) (2.38) Chapter 2. Numerical Procedure 40 Group Meaning Definition Rep Inertial / Hydrodynamic V Elastic / Hydrodynamic E/r/G Aspect Ratio rc = N a* Wall Thickness Ratio a* — a.i/a0 Table 2.3: Dimensionless groups and the parameters involved in their formulation The dynamic behaviour of a fibre is influenced strongly by its shape. Higher aspect ratio fibres rotate at a slower rate in a shear flow than lower aspect ratio fibres. This geometric factor is represented by rc, the cylindrical aspect ratio which is equal to the length of the fibre divided by its width. Since the elements used in the model have equal length and dimameter the aspect ratio will simply be equal to N, the number of elements. Chapter 2. Numerical Procedure 41 2.3.2 Computational Algorithm This section outlines the algorithm that was followed to perform simulations. It closely follows the code listing given in Appendix F. 1. PreUminary Declarations Before beginning calculations, all variables other than counters are declared as double precision. The number of elements, N is set as a parameter. Array declarations are then based on this parameter. The number of iterations, along with the dimensionless time step, are chosen to simulate the desired total time for the fibre motion. The largest time step allowable is chosen in order to reduce computational time. The issues of calculation time and time step limitations will be discussed in the Results & Discussion chapter. 2. Parameters The parameters which define the fibre (wall density, inner and outer radius, and Young's modulus) and the fluid model (water density, shear rate and viscosity) models are set according to the trial being performed. These values are then used to calculate the cross- sectional area, second moment of area, mass, mass moment of inertia, density, bending and spring constant, and the translational and rotational constants. 3. Initial Conditions Before entering the computational loop, the initial conditions of the fibre model are set. In all cases, the fibre is positioned initially with its centreline aligned with the y axis. The importance of the initial position will be discussed in the next chapter. The vertical placement of the fibre is defined by the location of the centre of the first element. Re- maining elements are connected above the first. Internal stresses are set to zero as the fibre is straight and therefore undeformed. The velocity and rotation of the elements are initialized to zero. Chapter 2. Numerical Procedure 42 4. Data Files Before calculations begin, the data files are created and opened. The primary file stores the position of each element at selected time intervals. Other files may be opened to store data on the forces, torques, et cetera acting on the fibre. If the model is continuing the results from a previous trial, a file will be read which contains the data from the last time step of the previous trial. 5. Computational Loop (a) At t = 0, the shear flow begins to act on the motionless fibre and the computational loop begins. A series of calculations are performed which depend on values from the previous time step (or initial conditions, at t = 0). The normal directions (Eq. 2.10), spring forces (Eq. 2.9), and bending torques (Eq. 2.12) are determined from the position of the elements relative to each other. The hydrodynamic force and torque are calculated by Eqs. 2.20 and 2.21 from the relative velocity between the elements and the fluid. (b) Equations of motion (Eqs. 2.22 and 2.23) are then formed. The accelerations are due to the forces and torques calculated above. The velocities from the previous time step are saved. New velocities are found by adding the contribution of the accelerations over the time step, u(i) = u0(i) + u(i)St (2.39) where u0(i) is the velocity of the element from the previous time step. (c) The velocity v} due to the tangential force ft is found by solving (N — 1) simultaneous equations. These equations are constructed by finding ti* from Eqs. 2.22 and 2.23 and substituting it into Eq. 2.24. A matrix is then formed from the coefficients of u*. The development of this matrix is found in Appendix D. The tridiagonal structure of the matrix decreases the computations required in the inversion subroutine tridiagonal. Chapter 2. Numerical Procedure 43 The contributions of /* to the translational and angular velocity of each element are then added. A further adjustments are made to the angular velocity so that the model does not break the geometric condition in Eq. 2.25. A torque, t9 is calculated by applying Eq. 2.26 to the angular values from the beginning of the time step. An additional rotational velocity, u9 due to the torque is found. The angular velocities are then updated to include this contribution. When the velocities have been found at the new time step, the displacement of each element over the time step may be calculated. This is done by considering the average velocity over the time step. For the x direction this calculated as, x(i) = x(i) + 0.5 [uQ(i) + u(i)] St (2.40) where u0(i) is the velocity at the beginning of the time step and u(i) is the velocity at the end of the time step. The new position of each element may then be recorded in a data file. 6. Continue Loop or End At this point, the program returns to the beginning of the loop or finishes. The time step is advanced, and a new series of calculations and matrix manipulations are begun. Once the loop is finished, the values of all variables used in the program are stored in a final data file. This file may be used if the simulation is to be continued. By repeating the computational loop over an adequate number of time steps, the motion of the elements and hence the fibre may be followed. Chapter 3 Results and Discussion 3.1 Introduction A series of trials have been conducted to simulate rotational and translational motion of a single fibre using the numerical model developed in the previous chapter. A preliminary series of tests was conducted with a rigid 10 element model to compare the present calculations with previously published results. These trials were also run to determine the behaviour of the model over a range of values for the inertial constants CtT and C r ot- Subsequent trials were performed with a 60 element model. The second series of trials was conducted over a range of values for the elastic constant kei. Results were determined for fibres with and without bending deformation. A final series was completed where the fibre model had localized weaknesses (kinks) in the bending stiffness along its length. The results of these trials as well particular aspects of the model's behaviour will be discussed in this chapter. Motion is induced in a fibre by placing it in a Couette flow where it begins at rest with zero velocity as described earlier. After time zero, the shear flow exerts a force on the fibre which responds by rotating in a clockwise direction. The subsequent motion of the fibre is determined by its initial position. A preliminary test was conducted with the fibre centred on the origin along the y-axis. However the remaining tests were run with the fibre aligned with the y-axis extending upward from the origin. This resulted in both rotational and translational motion. A schematic diagram of the two initial positions is found in Figure 3.15. The recurring feature in all trials is the presence of a shearing flow which causes the fibre model to rotate. It follows that the fundamental determinant of fibre motion is G, the shear rate of the flow field and that the dimensionless time unit is t* = tG. 44 Chapter 3. Results and Discussion 4 5 Figure 3.15: The in i t ia l position of fibres used in the trials. The fibre in a) wi l l rotate about the origin when a shear is exerted on i t . The fibre in b ) extending up from the origin, w i l l translate and rotate. The response of the fibre to the shear can be measured by T , the period of rotat ion. For this reason, many of the results discussed in this chapter wi l l be related to TG, the dimensionless t ime required to perform one period of rotat ion. The computat ional t ime for the trials was determined by the number of time steps required to calculate the period. Most trials were executed for the time required to complete 360° of rotat ion from the starting posit ion. C o n s t a n t s a n d P a r a m e t e r s The behaviour of the fibre model is determined by the forces and torques exerted on i t . In the previous chapter, it was shown that three forces are involved: inert ial , hydrodynamic, and elastic. The influence of the these forces is determined by the value of the the dimensionless groups shown i n Table 2.3 in the previous chapter. These groups in turn take their values from the parameters listed in Table 2.2. Chapter 3. Results and Discussion 46 Parameter Values Used in trials Fluid Density, pf 1000 kg/m3 Viscosity, 77 1Q~3 kg/ms Shear Rate, G 10, 1000 s'1 Outer Radius, aQ 10,100 /xm Fibre Wall Density, pfw 10- 104 kg/m3 Young's Modulus, E 100 - 2 X 106 Pa Aspect Ratio, N 10, 60 Table 3.4: The constants and values chosen for the parameters. The fibre wall density, pfw was kept constant at IbWkg/m3 in the trials with the 60 element model. The values chosen for the parameters depended on the fibre and flow properties being modelled. Several parameters were kept constant for most trials. The fluid properties were those of water at 20° Celsius so that the fluid density, pf and viscosity, 77 were set at 1000kg/m3 and 10~3kg/ms respectively for all trials. Again, it should be emphasized that the fluid density plays no role in the hydrodynamic force exerted on the fibre. It solely affects the fibre's inertia due to the fluid immobilized inside the fibre walls. The fluid velocity field was then determined by one parameter G, the shear rate. A fibre wall density, pfw of 1500%/m3 was used in most trials. The fibre therefore had a specific gravity, S of between 1.0 and 1.5, depending on the wall thicknes. The overall dimension of the fibre was set by the outer radius, aQ. Cylindrical elements had a length and diameter of 2aa matching those of spherical elements. The fibre length is then equal to 2a0iV where N is the number of elements in the fibre model. Table 3.4 summarizes the constants and parameter ranges used in the trials. A similar table will be shown at the beginning of each section giving the parameters used in that particular trial. Chapter 3. Results and Discussion 47 Format of Results Data was collected from the model in a Lagrangian fashion. The main data file recorded the position of the fibre. Since the fibre's position is represented by the midpoint of each element, the results displayed here appear half an element shorter than the actual fibre model at each end. Other data files recorded the forces acting in the fibre during the trial. At the end of each trial, a final data file was created which stored all variable values so that the model could be restarted if needed. Each data file was divided into separate zones for each specified time step. In theory, the number of zones is limited by the size of the dimensionless time step, St and the length of the trial. In practice, a zone was written to the data file every 0.1 t* in order to keep the file at a manageable size. This resolution was found to be high enough to observe the behaviour of the model in sufficient detail. The data representations in this chapter were created using Tecplot. This program processes ASCII data files to produce binary plot files allowing the data to be graphed in the manner desired by the user. This software was also used to create animated results of the fibre motion. Tecplot limits each plot file to a maximum of one hundred zones. This limitation required that plot files to be formed by taking zones from a portion of, or at regular intervals from the parent data file for each trial. Results from each trial were analyzed to determine the model's dynamic behaviour. For rigid fibres, the angle <j> formed by the fibre with the y axis was plotted against time and the period TG was determined. The period 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 <p = 180° are exhibited in Fig. 3.19. The time taken to reach <p = 180° is listed above each plot. The distinguishing feature of Trials 3(b) and 3(c) was the bending that occurred in the model. The amount of deflection in Trial 3(b) amounted to 0.184 of a fibrelength and increased to 0.242 of a fibrelength in Trial 3(c). These deflections are due to the increase in the component of fluid forces exerted along the length of the fibre as it rotates from the horizontal to the vertical. Chapter 3. Results and Discussion 59 (a) (b) (c) 1.0 0.5 0.0 t =128 t =138 64 65 66 67 68 69 64 65 66 67 68 69 Figure 3.20: The motion of the fibre for the three stiffest trials: (a) kei = 2 X 106, (b) kei = 106, and (c) kel = 5 X 105. The axes are measured in fibre lengths. Fig. 3.20 plots the fibre centrelines for Trials 3(a), 3(b), and 3(c) as they pass through the "tumbling" region (100° < <p < 260°) between t* = 128 and 140. The fibres are recorded at intervals of Ai* = 1 to show the detail of the motion. The motion of Trial 3(b) could be classified as "springy" in the terminology of Forgacs and Mason [6]. A region of bent motion occurs between t* — 132 and 134. The majority of its motion was identical to Trial 3(a). The rotational behaviour of Trial 3(c) could be classified as a "loop turn". Bending increased between t* = 128 and 133. The fibre rotated as a rigid body between t* = 133 and 134 before recovering its straight shape at t* = 13. Due to the flexibility in Trial 3(c), a smaller fibre profile is exposed to the fluid flow in the x direction. This may have significance for the hydrodynamic drag of a flexible fibre. Chapter 3. Results and Discussion 60 Very Flexible Motion - Trials 3(d), 3(e) The elastic constant was reduced to kei = 105 and 104 for Trials 3(d) and 3(e) respectively. Fig. 3.21 plots the fibre centrelines at ten evenly spaced locations over Ai* = 18 during the tumbling region. The dimensionless time at the beginning and end of the region are shown above the fibres. The results from Trials 3(d) and 3(e) differ considerably from the trials in Series 3. Most notably, deformation was more prominent in the fibres and was not confined to the tumbling region but occurred over a wide range of the motion. This made it dimcult to identify the exact angle of rotation of the fibre. In Fig. 3.21 the fibre estimated to represent the (p = 180° position is identified by the time displayed above it. The periods for Trials 3(d) and 3(e) were thus estimated as T G = 210 and 148 respectively compared to 266.6 for Trial 3(a). The motion of the fibre in Trial 3(d) is similar to that seen in Trial 3(c) and could be described as a more compact loop turn. The reduction in the period between the trials is likely due to the increased deformation. In Fig. 3.21(a), it can be seen that the fibre ends were pushed further apart into the flow field by the buckling action around (p — 90°. This induced the tumbling motion in Trial 3(d) that is not seen in Fig. 3.20(c) until the fibre has rotated past <p = 100°. Chapter 3. Results and Discussion 61 (a) 1.0 0.5 0.0 t =98 t =110 0 t =116 49 50 51 52 53 54 55 56 57 58 (b) 1.0 0.5 0.0 t=62 t' = 74 t' = 80 31 32 33 34 35 36 37 38 39 40 Figure 3.21: The motion of the fibre for the two most flexible trials: (a) ke\ = 105 and (b) kei = 104. The axes are measured in fibre lengths. The motion of Trial 3(e) seen in Fig. 3.21(b) was dominated by fibre deformation. Almost no rigid body motion was visible. The fibre ends bent inwards and migrated toward the middle of the fibre. After passing the midpoint of the fibre at t* = 70 the fibre ends moved apart until the fibre assumed a straight geometry at approximately t* = 88. At this point the pattern seen between t* = 0 and 62 is repeated. The high flexibility of the model in this trial would likely result in a more complex asymmet- ric orbit for an actual wood fibre with the same properties. However, the symmetry inherent in the formulation of the model causes it to perform in this manner. This symmetry is difficult to create even in highly idealized experimental situations. This represents one of the simpli- fications of the model which was addressed in the next series of trials. Another weakness of the model is the free draining approximation. The close proximity of the various parts of the fibre in Trial 3(e) would result in some manner of hydrodynamic interaction. The model in its present form does not take this into account. Chapter 3. Results and Discussion 62 Conclusions for Series 3 Five trials were conducted in Series 3 to determine the range of motion for a 60 element fibre model with stiffness varying from kei = 2 X 106 to 104. The rigid 60 element fibre model, Trial 3(a), showed an increase in the tendency to align with the flow direction over that seen with the 10 element model. This behaviour follows that seen in experimental and theoretical work. The period of rotation for Trial 3(a) was 266.6. A slightly decrease was seen in Trials 3(b and 3(c) where T G = 266.2 and 264.2 respectively. Rigid body motion predominated in these trials. When the stiffness was reduced in Trials 3(d) and 3(e) the period dropped to T G = 210 and 148 respectively. This resulted from the deformation of the fibre which caused the tumbling motion to occur at an earlier point than in the stiffer trials. The period calculated from JefFery's equation was 250.23 using r* — 39.8. The variation of the period with flexibility is plotted in Fig. 3.22. 250 200 r 150 -o •o TG • Results form Series 3 - • Jeffery's equation 100 h 50 h 10s Elastic Constant (E /nG) Figure 3.22: Stiffness vs Reynolds Number for the five trials in Series 3. Chapter 3. Results and Discussion 63 i i i i i Li •j/ 3(a) = - - 3(b) 3(c) 3(d) t 3(e) Figure 3.23: Fibre Orbits for Series 3. The axes are marked in tenths of a fibrelength. It was observed in the trials that the fibre profile exposed to the flow direction during the tumbling region was reduced for the flexible fibres. This phenomena is shown here by plotting the orbits traced by the fibre tips. The orbits were found by having the point of view translate at the expected velocity for the fibre's centroid (ie. the average fluid velocity). The orbits shown in Fig. 3.23 are shifted slightly to the left due to the translational lag of the fibre from zero velocity. The maximum profile of the fibre exposed to the flow direction as a percentage of fibrelength decreases from 100% for the rigid fibre in Trial 3(a) to 90.2, 57.7, 32.4, and 17.1% for Trials 3(b) to 3(e) respectively. This would suggest that the most flexible fibres have the lowest profile. However, if the profile is averaged over a whole rotation it is likely that the most flexible fibres will have a higher profile due to the deformation that exists throughout the rotation. Stiffer fibres would have a much smaller profile due to the large proportion of time they spend aligned with the flow. Another consequence of the orbits is the smaller effective area occupied by flexible fibres while rotating. This could reduce the possibility of interaction between fibres in a suspension of flexible fibres in a simple shear flow compared with rigid fibres at the same concentration. Chapter 3. Results and Discussion 64 3.3.2 Localised Wall Weakness - Series 4 Series 4 consisted of five trials where a localized weakness (kink) in the bending stiffness was placed at five evenly spaced locations along the fibre. These trials were conducted to address the deficiencies of the symmetric formulation of the model. The conditions for these trials were the same as Trial 3(c). The parameters used in this series are listed in Table 3.11. The fibre kinks were modelled by reducing the bending stiffness at one of 5 locations along the fibre for each trial. The bending stiffness, k* at the kink was reduced to 20% of the global value. Element Radius, a0 10 \xm Shear Rate, G 1000 s~l Wall Density, pjw 1500 kg/md E/VG 5 x 105 Table 3.11: These values were kept constant in the trials run in Series 4 The results for the trials are shown in Fig. 3.24 for the tumbling region between t* = 128 and 140. It can be seen from the figure that Trials 4(c) and 3(c) yielded similar results. The period of all the trials in Series 4 was unchanged from Trial 3(c). Figs. 3.24 shows that a kink caused the model to reproduce a "snake turn" as described in the literature as opposed to the symmetric "loop turn" seen in Trial 3(c). The motions in Trials 4(a) and 4(e) were identical but 180° of phase. The same was true for Trials 4(b) and 4(d). The kink in Trial 4(a) allowed one end of the fibre to buckle first but the bend did not occur at the kink as was the case in Trial 4(b). The importance of these trials lies in demonstrating that the model can simulate the actual motion seen in wood fibres by applying a local modification to the bending stiffness. Further tests could be run at the lower bending stiffnesses used in Trials 3(d) and 3(e) to determine if the periods of the rotation found in these trials are effected by a kink. Chapter 3. Results and Discussion 65 Figure 3.24: Fibre motion over the tumbling region for the kink located at: a) 1/6, b) 1/3, c) 1/2 d) 2/3, and e) 5/6 from the left end of the fibre. The arrow points to the location of the kink. The axes are measured in fibrelengths. Chapter 3. Results and Discussion 66 3.4 Discussion To this point, results have been presented separately for fibres of different aspect ratios with rigid or flexible motion. A single dimensionless number is required which can identify the range of motion in a fibre for any fluid or fibre conditions. A simple dimensional analysis is used in the development of this number below. Fibre Bending Number, B First, consider the viscous forces exerted on the fibre due to the shear flow, where n is the viscosity, a c is the outer radius of the fibre,and G is the shear rate. The viscous force per unit length on each element is then, ^ oc VGa0 Now, consider the deflection associated with the buckling or bending of the fibre where it is considered a beam or column, wL* 0 oc — — EI where 8 is the deflection, w is the load per unit length, L is the length of the fibre, and EI is the stiffness. For a critical deflection, (^) the load per unit length is, EI W K V Chapter 3. Results and Discussion 67 A dimensionless bending ratio, B can then be formed, B = ^ (3.4!) nGa0 = &&T? (3'42) Hence, a fibre is stiff (will not deform) if B exceeds a critical value. A more exact ratio for B may be formed by considering the theory of Forgacs and Mason [13]. This theory, developed in the Introduction, considers the case of a fibre where the viscous forces due to the fluid vary along its length. Eq. 1.8 is repeated below. E % 2r 4 >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 <ft dt* Appendix C. Dimensionless Parameters and Equations of Motion Hydrodynamic Terms Fh - -6irr)a0 (u - V) Fh* = -67T77a0 (tt - Gy) ( ^ j g ) = -6TT (u-Gy) (—p;) = -6TT (ti* - y*) Sirr/a3, (UJ — fi) rph* -8*r,al(u + G/2)(—) -8* ( « + f ) ( ^ ) = -8?r («* + 1/2) Appendix C. Dimensionless Parameters and Equations of Motion 84 Structural Terms F' = -ks (I - 2a0) F3* = -k*s {I* - 2) 7r . a 2 . k, = 2 ^ ( 1 - ^ ) 2 a0r/(j a0 ~ 2 77G U a2 j T 6 = -ftt AO Th* = -k* AO H = §*:ci-$ v - * E n "<i Appendix D Velocity Boundary Condition Generalized Coordinates dr{ d6{ d-Ti dOj _ + 0 o _ X n , . = _ + a o _ X 7 l . i Generalized Non-Dimensionalized Coordinates dr* dOi dr* d0j dt* dt* 3 dt* dt* 3 Component of velocity perpendicular to 71^+1 (* is assumed) Ui X nyi — Vi X nxi — oji = v.i+i X nyi — Vi+\ X nxi + tv̂ +i Additional Velocities due to unknown force ft = FtX St Ui = Ui + {fi X nyi- / * _ ! X 7l1/i_i)/C'trans Vi = Vi - (/,• X nXi - / / _ ! X nXi-x)ICtrans "i = Ui-(fi+fU)/Crot 85 Appendix D. Velocity Boundary Condition 86 Application of the Condition at Interface i —> 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.


Citation Scheme:


Usage Statistics

Country Views Downloads
China 1 0
City Views Downloads
Beijing 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}


Share to:


Related Items