The Orientation of Semi-Dilute Rigid Fibre Suspensions in a Linearly Contracting Channel by Paul Joseph Krochak B.A.Sc., The University of British Columbia, 2004 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy In The Faculty of Graduate Studies (Mechanical Engineering) The University Of British Columbia (Vancouver) May, 2008 © Paul Joseph Krochak 2008 Abstract This thesis has investigated the effects of long range hydrodynamic fibre-fibre in- teractions on the orientation state of a semi-dilute, rigid fibre suspension flowing through a linear contracting channel under laminar flow conditions. The effects of fibre-fibre interactions have been modeled mathematically, the governing equations solved numerically and the predicted results compared with experimental observa- tions. The theoretical model is based on the assumption that the orientation state of the suspension can be completely described by a probability distribution function and that fibre-fibre interactions are random in nature, thus giving rise to a diffusion-type process. The orientation distribution evolves spatially according to a Fokker-Plank type equation using closure equations for the rotary diffusion coefficient advanced by either (i) Folgar and Tucker (J. Reinforced Plast. Comp. 3 98–119 1984) or (ii) Koch (Phys. Fluids 7(8) 2086–2088 1995). Each of these two closure models for the rotary diffusion coefficient contain an unknown empirical constant that must be determined from experiments. These were fit to experimental data along the central streamline of the contraction as a function of fibre concentration. In both cases the diffusion coefficient was found to first increase with increasing suspension concentration up ii Abstract to a maximum, and then decrease with concentration above this point. This non- monotonic behavior was attributed to fibre clumping or flocculation, a mechanism not considered in the relationships for the rotary diffusion coefficient. The theoretical model was then extended to predict fibre orientation over the en- tire plane of the contracting channel and the two-way momentum coupling between the fluid and fibre phases was investigated numerically. The results show that the structure of the flow field within the contraction is significantly altered when the fibre phase is considered, demonstrating the non-negligible effect of the momentum exchange between the two phases. Comparison is made between the predicted ori- entation state of the suspension with experimental observations over the contraction plane. Good agreement is found between the model predictions and the experimen- tal observations, except in a small region near the solid boundaries where the two differed significantly. These near wall discrepancies are attributed to an inability to correctly handle the wall boundary conditions in the fibre orientation model. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 Description of Fibre Orientation . . . . . . . . . . . . . . . . . . . . 6 2.2 Equations of Motion for the Fluid Phase . . . . . . . . . . . . . . . . 13 2.3 Numerical Investigations . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4 Experimental Investigations . . . . . . . . . . . . . . . . . . . . . . . 17 2.5 Previous Headbox Studies . . . . . . . . . . . . . . . . . . . . . . . . 18 2.6 Fibre Orientation Near Solid Boundaries . . . . . . . . . . . . . . . . 20 2.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 iv Table of Contents 3 Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4 Estimating Dr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.2 Experimental Detail . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2.1 Numerical Estimates of Dr . . . . . . . . . . . . . . . . . . . 38 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.3.1 Comparison of Rotary Diffusion Models . . . . . . . . . . . . 45 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5 Validating Dr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.1 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.2 2D Orientation Results . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.2.1 Theoretical Predictions . . . . . . . . . . . . . . . . . . . . . 55 5.2.2 Experimental Observations . . . . . . . . . . . . . . . . . . . 57 5.3 Discussion of Fibre Orientation Results . . . . . . . . . . . . . . . . 60 5.4 Effect of 2 Way Coupling . . . . . . . . . . . . . . . . . . . . . . . . 62 5.4.1 Flow Field . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.4.2 Orientation Distribution . . . . . . . . . . . . . . . . . . . . . 65 6 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 68 6.1 General Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 6.2 Diffusion Coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.2.1 Folgar-Tucker Model . . . . . . . . . . . . . . . . . . . . . . . 70 6.2.2 Koch Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 v Table of Contents 6.3 Comparison of Diffusion Models . . . . . . . . . . . . . . . . . . . . 71 6.4 2D Fibre Orientation . . . . . . . . . . . . . . . . . . . . . . . . . . 71 6.5 Two Way Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 7 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 7.1 Mechanical Fibre-Fibre Interactions . . . . . . . . . . . . . . . . . . 74 7.2 Rotary Diffusion Coefficient . . . . . . . . . . . . . . . . . . . . . . . 75 7.3 Near Wall Behavior of Fibres . . . . . . . . . . . . . . . . . . . . . . 75 7.4 Two way coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Appendices A Experimental System . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 B Matlab Program Code . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 B.1 Fibre Tracking and Analysis . . . . . . . . . . . . . . . . . . . . . . . 87 B.1.1 Fibre Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . 87 B.1.2 Analysis of Fibre Observations . . . . . . . . . . . . . . . . . 93 B.2 Numerical Prediction of Fibre Orientation . . . . . . . . . . . . . . . 97 B.2.1 Fokker-Plank Equation Solvers . . . . . . . . . . . . . . . . . 97 B.2.2 Computation of the Fibre Stress . . . . . . . . . . . . . . . . 106 vi List of Figures 1.1 The channel geometry and fibre orientation angles used in this study. 2 4.1 The orientation of a fibre with respect to flow in the contraction. The angle φ is the angle of the projection of the fibre in the xy-plane and θ is the angle between the fibre and the z-axis. The angle γ is the angle of the projection in the xz-plane . . . . . . . . . . . . . . . . . . . . . 29 4.2 An image of the experimental test section. . . . . . . . . . . . . . . . 33 4.3 Two consecutive unprocessed images of fibre observation separated by 200ms. The upper image pair sequentially precedes the bottom image pair. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.4 An example of an original image (left) and a post-processed image with white fibres on a black background (right). The fibre concentration is nL3 = 8, the fibre length is 5mm, and the aspect ratio is 50. . . . . . 36 4.5 The region bounded by the red lines represents the analysis region. All fibres observed in this area are used for the 1D analysis and de- termination of CI and λ. . . . . . . . . . . . . . . . . . . . . . . . . . 38 vii List of Figures 4.6 Development of the orientation distribution along the contraction cen- terline for the φ projection. The experimental data are given at x L = 0(◦), x L = 0.33(), x L = 0.66() and x L = 1.0(). The con- centration shown here is nL3 = 8 and the solution shown here was obtained with Equation 4.2. . . . . . . . . . . . . . . . . . . . . . . . 40 4.7 Plot of the orientation variance σ2 along the contraction length for nL3 = 8(), nL3 = 16(), nL3 = 24(◦), and nL3 = 32(). The uncertainty of the data are represented by the error bars which have been computed at the 95% confidence interval. . . . . . . . . . . . . . 41 4.8 Plot of CI(−−) and λ(−−) vs nL3. . . . . . . . . . . . . . . . . . 42 4.9 Shown is a comparison of λ nL 3 ln2r (−−) with CI(−−) vs nL3. . . . . 44 4.10 Development of the orientation distribution along the contraction cen- terline for the φ projection. The experimental data are given at x L = 0(◦), x L = 0.33(), x L = 0.66() and x L = 1.0(). The best fit of Equation 4.4 is given as the solid lines and the predictions given as the hatched lines. The left column is the solution using Equation 4.2, and the right column is the solution using Equation 2.11. . . . . 46 5.1 Plot of the L2 norm of the relative change in the flow field between successive iterations. The concentration shown here is nL3 = 8. . . . 53 5.2 Plot of the L2 norm of the relative change in solutions along the central streamline between solutions obtained with mesh sizes of 736 x 500 and 1397 x 1000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 viii List of Figures 5.3 Numerical predictions of the mean fibre orientation angle. Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom). . . . . . . 56 5.4 Numerical predictions of the orientation variance, σ2. Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom). . . . . . . 58 5.5 Experimental observations of the mean fibre orientation angle. Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom) . . . . 59 5.6 Experimental observation of the orientation variance σ2. Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom) . . . . . . 60 5.7 Contour plots of numerically predicted velocity magnitude. Shown are the pure solvent (Top), nL3 = 8 (Mid-Top), nL3 = 16 (Mid-Bottom) and nL3 = 24 (Bottom). . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.8 Comparison of the orientation distribution functions obtained with and without the two way coupling at the end of the contraction. Shown on the right hand side is a close up near the peak of the distribution functions. Shown are nL3 = 8 (Top), nL3 = 16 (Middle), nL3 = 24 (Bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 A.1 A schematic diagram of the experimental system. . . . . . . . . . . . 83 A.2 A dimensioned sketch of the experimental section. . . . . . . . . . . . 84 ix List of Tables 4.1 χ2 Values for Each Model . . . . . . . . . . . . . . . . . . . . . . . . 47 5.1 Difference Between Coupled and Uncoupled Orientation Distributions 66 x Acknowledgments To Katie, for all her love and support. And, to my parents who have loved, encouraged and supported me my entire life, without whom this would never have been possible. I would also like to acknowledge my supervisors James Olson and Mark Martinez for their help with this project in addition to their support and inspiration over the years. Thank you very much for giving me this opportunity. i Chapter 1 Introduction In this work, the effect of long range hydrodynamic fibre-fibre interactions on the orientation state of rigid fibres in a viscous fluid flowing through a linear contracting channel are investigated. The case considered here is for flow with dilute to semi- dilute suspensions where the Reynolds numbers based on the length of the fibre is asymptotically small and based on the inlet channel width is order one. The monodispersed fibres studied here are cylindrical in shape with a length and diameter of 5mm and 100μm respectively; the aspect ratio of the fibres, r, that is the ratio of length to diameter is 50. The concentration limits of interest in this work are defined such that 1 ≤ nL3 ≤ r, where n is the number of fibres per unit volume and L is the fibre length. The contraction consists of two rigid sloping walls separated by 50.8mm at the entry and 5.08mm at the exit and the length of the contraction is 130mm, see Figure 1. Controlling the orientation state of fibres in contraction-type flows is one of the prin- cipal objectives in many science and industrial processes. For example, during the electrospinning of polymer based nanofibres, a polymer suspension is fed through a small contraction called a spinneret. A strong electric field is then applied which 1 Chapter 1. Introduction Figure 1.1: The channel geometry and fibre orientation angles used in this study. draws the suspension out through the contraction. This serves to align the polymer molecules in the flow direction. The motivation behind this thesis however is based on a similar but much larger-scale industrial process, namely papermaking. During papermaking, a flocculated fibre suspension is fluidized by turbulence created locally from a sudden change in geometry in an apparatus called a headbox. The fluidized fibre suspension is then formed into a plane liquid jet through an elongational flow in a contracting channel and then distributed on a moving permeable band or wire from which the water is drained and the paper is formed. Despite the longstanding use of headboxes, there are still many unanswered fundamental questions on the relationship between the operating state and the final fibre orientation distribution. One of the essential scientific difficulties is the incomplete understanding of the long- 2 Chapter 1. Introduction range hydrodynamic interactions between fibres in the suspension. A hydrodynamic fibre-fibre interaction is said to occur when one fibre is swept within the disturbance flow field of another fibre causing the fibre to undergo a small change in orientation (and vice-versa). This effect propagates throughout the suspension resulting in de- viations from the expected orientation state of the suspension. Another difficulty in understanding these flows is that momentum is exchanged between the fluid and fibre phases as the suspension flows. Both the magnitude and direction of this mo- mentum transfer are highly dependent on the orientation state of the suspension, a complication which effectively couples the structure of the flow with the orientation state of the suspension. Moreover, the momentum exchange between the two phases can affect the suspension rheology, viscous head loss and energy consumption within the headbox. What complicates both of these issues is that papermaking suspensions tend to mechanically clump or flocculate during processing. The ability to manipulate and control fibre orientation in the headbox is critical to forming higher quality paper. One of the key elements towards this is the ability to predict the effect of different system parameters on the orientation state. This means that a highly accurate and efficient numerical model must be available so that ini- tial predictions can be made before larger-scale prototypes are designed and tested. In this work, the effects of long range hydrodynamic fibre-fibre interactions on the orientation state of dilute to semi-dilute rigid fibre suspensions flowing through a linear contracting channel under laminar flow conditions is investigated numerically and experimentally. The effects of fibre-fibre interactions are modeled mathemat- ically by assuming that the orientation state of the suspension can be completely 3 Chapter 1. Introduction described by a probability distribution function and that fibre-fibre interactions are random in nature, thus giving rise to a diffusion-type process. The spatial evolution of the orientation distribution obeys a Fokker-Plank type equation where two dif- ferent relationships for the rotary diffusion coefficient are examined; one is a simple model that is easy to perform numerical predictions with but contains very little physical meaning, the other a rigorous model which encompasses the underlying physics of the problem but comes at a much higher computational cost. The two models are compared with experimental observations of fibre orientation along the channel length. Also studied is the effect of the two way momentum coupling on the orientation state of the suspension and on the structure of the flow field. This is accomplished by including the additional stress arising from the fibre phase in the fluid momentum equatons. The ultimate goal of this work is to improve the current state of fibre orientation modelling in the headbox. Modelling fibre motion is a formidable task. Perhaps the greatest obstacle relating to this problem is the extremely large computational domain required to solve the prob- lem regardless of the specific modelling technique being used. Furthermore, even the best available models require some apriori information that can only be determined from experiments. Approximation methods have been developed to circumvent the computational domain issue, but again, these methods are only approximate and also rely on additional apriori information about the system; information that can only be gained through experimentation or through direct computations, which again would require the enormous computational resources. A review of the literature pertinent to this study is presented in §2. Included in 4 Chapter 1. Introduction this section is a detailed development of the theoretical basis for the fibre orientation model and for the momentum exchange between the fluid and fibre phases. The objectives of this thesis are given in §3. In §4 the details on the experimental and numerical procedures used to measure the orientation distribution and to estimate Dr are presented. Included in this section are the the estimates for CI and λ, that is, the unknown closure parameters for the rotary diffusion coefficient. Also included within this section is a characterization of the experimentally observed orientation distribution along the contraction and a comparison of the two models for the rotary diffusion coefficient. In §5, the theoretical model is extended to predict fibre orientation over the entire xy-plane of the contracting channel. The effect of the two-way momentum coupling is investigated numerically and comparison is made between the predicted orientation state of the suspension and experimental observations over the xy-plane. Discussion on the effects of the two-way coupling and discrepancies between the numerical predictions and experimental observations of fibre orientation over the xy-plane are discussed. 5 Chapter 2 Literature Review In this section a brief overview of the literature regarding the evolution of the orientation-state of flowing fibre suspensions is given. In §2.1, we begin the dis- cussion by defining the traditional means of describing fibre orientation, that is, the Fokker-Plank equation and discuss the inherent closure problem associated with its use. In §2.2, we continue the discussion by describing the effect of the fibres on the carrier fluid and describe an additional stress imposed on the fluid. In §2.3 a de- scription of recent simulations in the general literature is given. This is followed by a review of experimental works given in §2.4. In §2.5 studies relevant to pulp and paper are presented. Finally in §2.6, we discuss one open question which limits the utility of these works, that is, the effect of the wall on the orientation distribution. 2.1 Description of Fibre Orientation The earliest study on the motion of a single isolated fibre dates back to Jefferey (1922) where he derived an equation of motion for a single rigid ellipsoid in a viscous unbounded Newtonian fluid. Using a no-slip boundary condition at the fibre’s surface 6 Chapter 2. Literature Review and matching the velocity field in the region near the fibre to the bulk flow field of the suspending medium, Jefferey solved for the flow field around the particle. A formula for the angular velocity vector of the fibre was then derived through a torque balance about the centre of the fibre. Jefferey’s derivation showed that a fibre rotates in one of a family of closed orbits around the vorticity axis termed the Jefferey orbits. Bretherton (1962) showed that Jefferey’s equations apply to cylindrical particles if the particle aspect ratio is replaced by an effective aspect ratio. What is clear from these studies is that fibers rotate continuously in response to shear within the flow and that there is no preferred orientation distribution. For cases above this dilute limit, fibre-fibre interactions create deviations from the expected orientation state. As such, modelling the orientation state of non-dilute suspensions must be handled differently. With shear flows, a number of authors (Altan et al. (1989), Folgar & Tucker (1984), Koch (1995), Leal & Hinch (1971), Lin & Zhang (2002), Parsheh et al. (2005)), describe the evolution of fibre orientation by a probability density function Ψ that obeys a Fokker-Planck relationship, that is DΨ Dt = 1 sin(θ) ∂ ∂θ (Drsin(θ) ∂Ψ ∂θ − sin(θ)θ̇Ψ) + 1 sin(θ) ∂ ∂φ (Dr 1 sin(θ) ∂Ψ ∂φ − sin(θ)φ̇Ψ) (2.1) where D Dt is the total derivative with fluid velocity field components; and u and v are the velocity components in the x and y directions respectively. The fibre orientation 7 Chapter 2. Literature Review vector, p, is defined as a unit vector pointing parallel to to the major axis of a fibre which, in three dimensional Euclidean space, is defined as p = ⎡ ⎢⎢⎢⎢⎣ cosφ sin θ sinφ sin θ cos θ ⎤ ⎥⎥⎥⎥⎦ (2.2) where φ and θ are the azimuth and zenith orientation angles respectively (see Figure 1.1) with corresponding unit vectors given by θ̂ = ⎡ ⎢⎢⎢⎢⎣ cosφ cos θ sinφ cos θ − sin θ ⎤ ⎥⎥⎥⎥⎦ and φ̂ = ⎡ ⎢⎢⎢⎢⎣ − sinφ cosφ 0 ⎤ ⎥⎥⎥⎥⎦ The rotational velocity can then be expressed as ṗ = φ̇ sin θφ̂ + θ̇θ̂ (2.3) Olson et al (2004) have derived the following expressions for the rotational velocities in terms of orientation angles φ̇ = 12 Lf 3 sin θ ∫ Lf/2 −Lf/2 lφ̂ · u(y + lp, t)dl (2.4) 8 Chapter 2. Literature Review and θ̇ = 12 Lf 3 ∫ Lf/2 −Lf/2 lθ̂ · u(y + lp, t)dl (2.5) In application to headbox flows where fibre length is small in comparison to the scale of the headbox, it is reasonable to locally linearize the flow field by taking a Taylor series expansion of the velocity field around the fibre centre, i.e., u(y + lp) ≈ u(y) + ∂u(y) ∂y lp (2.6) and substitute it into Equations 2.4 and 2.5 above. This results in the following simplified equations, θ̇ = θ̂ ( ∂u(y) ∂y p ) (2.7) and φ̇ = φ̂ sin(θ) ( ∂u(y) ∂y p ) (2.8) or expressing φ̇ in terms of its component partial derivatives for planar flow where u = u(u(x, y), v(x, y)), we have φ̇ = 1 2 ( ∂v ∂y − ∂u ∂x ) sin(2φ)− ∂u ∂y sin2(φ) + ∂v ∂x cos2(φ) (2.9) and 9 Chapter 2. Literature Review θ̇ = ∂u ∂x cos2(φ) cos(θ) sin(θ) + ∂u ∂y sin(φ) cos(φ) sin(θ) cos(θ) + ∂v ∂x cos(φ) sin(φ) sin2(θ) + ∂v ∂y sin2(φ) sin2(θ) (2.10) These equations are identical to those derived by Jefferey (1922) for fibres of large aspect ratio. The use of an orientation distribution function to describe the orientation state of a fibre suspension was further supported by Doi & Edwards (1986), who used stochastic arguments to predict the orientation state of rigid polymer suspensions. They argued that any large molecule in a suspension will experience a deterministic force and/or torque due to the flow field, but will also undergo random fluctuations in orientation due to Brownian rotary motion. Hence, they modeled the orientation state of polymer suspensions using a Fokker-Plank equation with a rotary diffusivity term, that is, Dr in Equation 2.1, which results from the effects of Brownian diffusion on the orientation state of the suspension. Folgar & Tucker (1984) used the Fokker-Plank equation to model the orientation state of rigid cylindrical particles in a semi-dilute suspension, however the rotary diffusivity term in their model was not due to Brownian effects, but rather was assumed to be the result of random hydrodynamic interactions with neighboring fibres; that is, hydrodynamic fibre-fibre interactions which cause fibres to deviate from their Jefferey orbits. They proposed, through dimensional analysis, a simple 10 Chapter 2. Literature Review relationship in which Dr is linearly proportional to the second invariant of the strain rate ‖α̇‖. For one dimensional flow in a linear contraction, Dr can be expressed as follows Dr = CI α̇ (2.11) where CI is traditionally called the interaction coefficient and is related to a number of suspension parameters such as concentration, aspect ratio, and fibre length; it is an empirical constant and must be determined from experiment. Jackson, Advani & Tucker (1986) used the Folgar-Tucker (FT) model to predict the 2-D orientation behavior of short fibres in compression moldings and showed that the model compares well with experiments once the interaction coefficient is determined. In studies by VerWeyst & Tucker (2002), the FT model for the rotary diffusion coefficient was used to investigate the effect of the two way coupling between the fluid and fibre phases for flows in complex geometries. Despite the fact that this model has been widely used with considerable success, Ranganathan & Advani (1991) discovered one major short-coming with the model. They used the Folgar-Tucker model to relate the interaction coefficient, CI , to the average interparticle spacing in non-dilute suspensions. They found a strong depen- dence between CI and the orientation state of the suspension. More specifically, by fitting the experimental data of Folgar (1983) for nL3 = 52.2 to Equation 5.1, where n is the number of fibres per unit volume and L is the fibre length, they found that CI was an order of magnitude larger when the suspension was in an isotropic orientation state, as opposed to an aligned orientation state. This phenomenon was 11 Chapter 2. Literature Review not however found to be significant when nL3 = 26. What these studies did show was an orientation state dependence of the interaction coefficient. Koch (1995) proposed a more physically complete form for the diffusion coefficient in three dimensional space and, like Ranganathan & Advani (1991), considered this parameter to vary as a function of the orientation state. With this in mind, he put forth the argument that the diffusion coefficient should be proportional to the product of the interaction rate and the average of the square of the orientation change per interaction. Koch proposed the following relationship Dr = nL3 α̇ ln2 r (λ1IE : 〈pppp〉 : E+ λ2E : 〈pppppp〉 : E) (2.12) where E is the fluid strain rate tensor, defined as E = 1 2 (∇u+∇Tu) (2.13) and r is the aspect ratio of the fibre, that is, the ratio of fibre length to diameter. The remaining term to be defined in 2.12 is the fourth order moment of the orientation distribution function Ψ, which is often referred to as the fourth order orientation tensor and is defined as 〈pppp〉 = ∫ ∫ pipjpkplΨ(φ)sin(θ)dθdφ (2.14) The first term on the right hand side in Equation 2.12 represents the isotropic rotary 12 Chapter 2. Literature Review diffusion coefficient that is proportional to the average of the square of the force per unit length on fibres. The second is an anisotropic correction term representing diffusion in the directions parallel to the force per unit length exerted by the other fibres. λ1 and λ2 are empirical constants which must be determined by experiment. Koch (1995) found that for simple shear λ1 = 3.16 × 10−3 and λ2 = 1.13 × 10−1; these values however may need to be verified for other flow fields, such as for the flow considered in this thesis. The main drawback to this model is that it introduces a non-linearity into Equation 2.1. Koch’s form for Dr has not been widely adopted, perhaps due to it’s mathematical and computational complexity and we find that the Folgar-Tucker model is found in most studies of laminar suspension flows, (e.g. Cintra & Tucker (1995), Lipscomb & Denn (1988), Ranganathan & Advani (1991), VerWeyst & Tucker (1999, 2002)). 2.2 Equations of Motion for the Fluid Phase The fact that the momentum is exchanged between the fluid and fibre phases in a flowing suspension has been known for some time now. Batchelor (1970) developed a general constitutive equation for the bulk stress in a suspension of rigid, inertia- less particles of arbitrary shape in a Newtonian fluid. This was accomplished by representing a single particle in suspension as a distribution of Stokeslets1 over a line enclosed by the particle body and then determining expressions for the resul- 1A Stokeslet is defined as a singularity in a Stokes flow which represents the resultant effect of a force applied to a fluid at that point. 13 Chapter 2. Literature Review tant force required to sustain translational motion and the resultant couple required to sustain rotational motion. These expressions were approximated in powers of , where = (log4r)−1 1. In a similar study, Cox (1970) determined the force per unit length exerted by the fluid on a cylindrical fibre in the limit as the Reynolds number tends to zero. He accomplished this slender body theory and generalized the results to an arbitrary, non-homogeneous flow field. He further showed that, to a first approximation the drag tensor is independent of fibre position and flow field. Khayat & Cox (1989) extended this result to small yet finite Reynolds number. Dinh & Armstrong (1984) extended Batchelor’s theory to account for the orientation state of elongated particles and its effect on the bulk stress within the suspension. This was accomplished by assuming the existence of a known orientation distribution function, that is, Ψ in Equation 2.1. By locally linearizing the flow field around the particle they were able to equate Batchelor’s constitutive equation to a new con- stitutive equation involving Ψ via a fourth order orientation moment tensor. By doing so they coupled the fibre orientation equations to the fluid momentum equa- tions. Bibbo, Dinh & Armstrong (2001) validated the work of Dinh & Armstrong (1984) by experimentally measuring the shear stress within a semi-concentrated sus- pension undergoing simple shear flow. Their measurements demonstrated a strong coupling between fibre alignment and the resulting shear stress within the suspen- sion and that the Dinh-Armstrong model was suitable for predicting the additional fibre stress. Shaqfeh, & Fredrickson (1990) used the Dinh-Armstrong model to derive an asymptotic form for the effective viscosity of a semi-dilute suspension of rods. The resulting relationship for the additional stress due to the fibre phase is as 14 Chapter 2. Literature Review follows τ fibre = μcr2E : 〈pppp〉 ln(1/c) + ln(ln(1/c)) + 1.439 (2.15) where c is the volume fraction of fibres within the suspension. With Equation 2.15 it becomes possible to simultaneously compute the fully coupled fibre orientation distribution and the resulting flow field. 2.3 Numerical Investigations More recently, researchers have been making great efforts to perform 3-D fibre ori- entation predictions for 2 and 3D flows inside complex geometries, as opposed to restricting their analysis to simple geometries only (e.g. Lin & Zhang (2002), Lip- scomb & Denn (1988), VerWeyst & Tucker (1999, 2002). The difficulty with this approach lies in the large computational domain required to resolve both the spatial and orientation domains when directly computing the orientation distribution of the suspension. For example, in order to model the orientation state of a three dimen- sional system, one would require a domain consisting of three spatial variables (e.g. the (x, y, z)-Cartesian coordinate system), and two orientation angle variables (e.g. the azimuthal and zenith angles in a polar system). If 50 grid points were to be used in each of the 5 variables, a total of 505 grid points would be required to solve the problem. This is on the order of 300 million nodes and would not nearly be a fine enough mesh for modelling many industrial suspension flows, particularly headbox 15 Chapter 2. Literature Review type flows. Despite the ongoing closure issues discussed in §2, there have still been many numeri- cal studies that incorporate the two-way coupling in conjunction with a Fokker-Plank equation to predict the orientation state of fibre suspensions in various complex ge- ometries. For the most part, these studies have shown good success. For example, Lipscomb & Denn (1988) used the tensor form of this model to predict the flow of dilute fibre suspensions through various axi-symmetric contractions using the finite element method (FEM). They showed that the suspension flow differs qualitatively from the flow of the suspending fluid alone. They also showed that their numerical results agreed well with experimental observation. VerWeyst & Tucker (1999) used a Galerkin finite element method to predict fibre orientation patterns in a variety of injection molded features using the Hele-Shaw flow model for non-dilute suspen- sion concentrations. They showed that their predictions compared very well with experiments when the detailed calculation can be reduced to two spatial dimensions and one orientation angle. They also found that 3-D calculations showed good qual- itative agreement with experimental observations. VerWeyst & Tucker (2002) used the fully coupled, Galerkin FEM solution to predict the flow of fibre suspensions through a variety of complex geometries which include axi-symmetric contractions, expansions and center-gated disks for a number of different suspension concentra- tions and Reynolds numbers. For the axi-symmetric expansion, they showed the formation of a large corner vortex when Re ≤ 5 and that this vortex grows in size with increasing fibre concentration. For Re ≥ 5, the size of the vortex decreased as the fibre concentration was increased. For the center-gated disk, they report that the 16 Chapter 2. Literature Review effect of the coupling was to displace streamlines toward the bottom surface as flow enters the disk region, and to align the fibres more rapidly in the radial direction. What is clear from these studies is that there is a significant difference between the coupled and uncoupled solution for the flow field and the orientation distribution. Lin, Gao, Zhou & Chan (2005) used the two way coupled Fokker-Plank formulation to predict the orientation state of a 2-D turbulent fibre suspension in a planar chan- nel flow. They ignored fibre-fibre interactions but used the rotary diffusion term to model the effect of turbulent fluctuations of fibre orientation. They showed that their numerical predictions compared well with experimental findings in 2-D. In studies by Paschkewitz et al. (2004), the drag reduction induced by the presence of rigid fibres in a turbulent flow was studied numerically. There work showed that drag reductions of up to 26% were theoretically possible for rigid, non-Brownian fibre suspensions in the semi-dilute concentration regime. 2.4 Experimental Investigations While there is an abundance of literature on the theoretical aspects of fibre orienta- tion, experimental investigations have been considerably fewer and far between. One of the earliest quantitative experimental studies was that by Anczurowski & Mason (1968) where they measured the periods of rotation of prolate spheroids in a Couette flow and then compared their results to the theoretical predictions of Jefferey (1922). These experiments were carried out for the projection of the azimuthal angle only, that is the fibre angle projected into the xy-plane. In a similar but more rigorous 17 Chapter 2. Literature Review study by Stover, Koch, & Cohen (1992), the orientations of rigid fibres in a semi- dilute, index-of-refraction-matched suspension were observed. These experiments were performed in the semi-dilute concentration regime and showed that even at the highest concentration investigated (i.e. nL3 = 45), fibres continue to rotate about the vorticity axis and remain highly aligned in the flow direction. Their measured orbital time constants were however found to be considerably different from those obtained by Anczurowski & Mason (1968), with the differences being attributed to the rotary diffusion resulting from long range hydrodynamic fibre-fibre interactions. In studies by Rahnama, Koch & Cohen (1995), an index of refraction matched system was used to directly measure the mean square orientational displacements in a dilute suspension of rigid fibres in an extensional flow. The mean square ori- entations were used to quantify the orientational diffusion resulting from fibre-fibre interactions, as well as to determine the effect of concentration on this diffusion-type process. Their results showed a clear increase in rotary diffusion with an increase in fibre concentration. 2.5 Previous Headbox Studies Within the pulp and paper literature several authors have investigated the orien- tation behavior of pulp fibres in headbox flows. In studies by Ullmar (1998) and Ullmar & Norman (1997), the fibre orientation distribution in the plane of the paper was measured at the headbox exit by digitally imaging nylon fibres in the suspension. These studies examined the effect of mean flow through the headbox, the headbox 18 Chapter 2. Literature Review contraction ratio, and fibre concentration on the orientation state of the suspension. In another study by Zhang (1998) the fibre orientation distribution was measured in both the plane of the paper and in the plane of contraction at several points along the axis of the headbox. Olson et al (2004) determined analytically the mean flow field along the centerline of a linear contraction, neglecting the two way coupling between the fluid and fibre phase, then used the Fokker-Plank equation to predict the 2-D orientation state for the resulting flow. They chose a global rotary diffusion coefficient which models the effect of turbulence on the orientation state by fitting their computations to the experimental data of Ullmar & Norman (1997) and Zhang (1998). Hyensjö et al. (2007) extended the work of Olson et al (2004) by using single phase CFD modeling to predict the flow field along streamlines of a linearly contracting channel. They used these results to compute the 2-D fibre orientation distribution along individual streamlines of the flow in the absence of fibre-fibre in- teractions and flow-fibre coupling. Parsheh et al. (2005) performed experimental studies to better understand the effect of turbulence on the orientation anisotropy of a highly dilute suspension of rigid rods in a planar contraction. They varied the turbulent intensity at the contraction inlet, then measured the orientation distribu- tion at a set of discrete points along the centerline of the contraction and finally compared these results to predictions based on a Fokker-Plank type equation. They showed that the rotational diffusion coefficient decays exponentially with respect to the local channel contraction ratio and depends on the inlet conditions. 19 Chapter 2. Literature Review 2.6 Fibre Orientation Near Solid Boundaries One major open question remains in the fibre orientation problem. That is, how do fibres behave near solid boundaries. This question has seen some research in the past years yet there have been no definitive answers. There are several possible orientation patterns for which a fibre may follow near a wall. The fibre could possibly align itself in the direction of the wall. The fibre could seemingly rotate, or flip about the wall. The fibre could shift away from the wall in response to the large flow gradients in this region or orient itself perpendicular to the wall. Stover, & Cohen (1990) were the first to rigorously address this issue. In their studies, the motion of rodlike particles in a low Reynolds number plane Poiseuille flow were observed experimentally and the effect of the wall on the period of fibre rotation was determined. They found that when particles with a high Jeffery orbit constant, (that is particles with a large period of rotation about the vorticity axis) came within a distance less than half a fibre length from the wall, an irreversible interaction between the fibre and wall occurred whereby the fibre was said to ”pole-vault” away from the wall to a distance of approximately half a fibre length. This interaction was accordingly named as it seemingly mimicked the flipping motion of pole-vaulter. Yet the fibre was said to never actually touch the wall but rather that the fibre tip would come within approximately one fibre diameter from the wall. They suggested the existence of some nonhydrodynamic force between the fibre and wall but could not determine the exact nature of this interaction. They observed quite a different behavior for fibres with a low Jeffery orbit constant that came within a half fibre length from the 20 Chapter 2. Literature Review wall. In this case, the fibre remained close to the wall indefinitely, however if the fibre had an increased period of rotation of approximately 50% it shifted its orientation away from the wall. In studies by Holm & Söderberg (2007) the influence of near wall shear stress on fibre orientation was investigated experimentally. They measured the resulting wall plane fibre orientation distribution of dilute suspensions in laminar flow down an inclined plane. Fibre length, aspect ratio and concentration were varied and images were collected along a set of wall parallel planes in an index-of-refraction matched system. They found that fibres were well aligned in the flow direction at distances ≥ 3mm from the wall, regardless of fibre length and aspect ratio. However, very close to the wall, that is at a distance of 1mm from the wall they found many fibres that were not aligned in the flow direction. For an aspect ratio of 40 there was a clear increased in the variance of the orientation distribution. For aspect ratios of 10, and 20, they found a bi-modal distributions, with one of the two modes at 0o relative to the flow direction and the other, more dominant mode at 90o relative to the stream direction. By increasing the fibre concentration, it was shown that the bi-modality of the distributions disapeared and the orientation distribution returned to having only one clear mean value at approximately 6o. 2.7 Summary From this body of works, we see that a complete description of the orientation state of a flowing suspension can be realized by solving a Fokker-Plank equation, that is, Equation 2.1 for the orientation distribution function, Ψ, where the effects of hy- 21 Chapter 2. Literature Review drodynamic fibre-fibre interactions are modeled with the rotary diffusion coefficient. Two models are available for the rotary diffusion coefficient, that is Equations 2.11 and 2.12 however both models contain unknown proportionality constants which must be determined before meaningful computations can be performed. Further- more the flow field is affected by the orientation state of the suspension creating a two way coupling between the fluid flow and the orientation state of the fibre phase. Experimental works in this area are far less evident, particularly those pertaining to contraction-type flows and the correct description of the near wall orientation behavior of the suspension is still uncertain. 22 Chapter 3 Problem Definition The description of the orientation state of rigid fibre suspensions in a variety of different flow conditions has received considerable attention. It is clear from the literature review that the focus of these works have been on the theoretical side and experimental works are less evident. From the above review, it should be emphasized that the equations describing the evolution of the orientation distribution function do not constitute a closed system. The success of a mathematical model in this case is dependent on the accuracy of the model for fibre-fibre interactions. Two models for the rotary diffusion resulting from fibre-fibre interactions are available; one a simple model that lacks physical meaning, the other a rigorous model that encompasses the physics of the problem but comes at a greater computational cost and has seen little use by other researchers. Furthermore, analysis of the effects of fibre-fibre interactions on the orientation state of suspensions subject to flow within linear contractions, that is, for flow within geometries typical of modern headboxes are non-existent. This leads to the first set of objectives for this research: 1. To experimentally measure the rotary diffusion coefficient as a function of fibre concentration for a linear contraction. 23 Chapter 3. Problem Definition 2. To determine the proportionality constants for each of the two models for the rotary diffusion coefficient. That is, the models of Folgar & Tucker (1984) and Koch (1995). 3. To compare the utility of these two models for predicting the suspension ori- entation state by comparison with experimental measurements. With the diffusion coefficient known, and a useful model at hand, it is possible to directly compute the orientation distribution function with good accuracy over a two dimensional spatial domain and without the need for orientation-tensor based closure approximations. This would prove to be very useful not only for predicting fibre orientation, but also for the further study of tensor-based closure approximations in contraction type flows. Also missing in the literature are investigations into the effects of the two way momentum coupling between the fluid and fibre phases on the suspension orientation state, and on the structure of the flow field for geometries similar to that of a headbox. Studies found in the literature have shown that the effect of the two way coupling is very significant in simpler flow geometries, however none have addressed its effect on suspension flow within a linearly contracting channel. Furthermore, the majority of the orientation analysis’ of linear contractions have been based on centerline calculations alone, an analysis that neglects the effects of 2 and 3 dimensional flow gradients on fibre alignment. This leads to the second set of objectives for this thesis: 1. To directly compute the 2D fully coupled planar orientation distribution func- tion over the entire plane of the contraction. 24 Chapter 3. Problem Definition 2. To investigate the effects of the two way coupling on both the flow field and orientation state of the suspension. Lastly, the numerical predictions of the suspension orientation state will be compared with 2D observations of fibre orientation over the entire xy-plane. 25 Chapter 4 Estimating Dr As described in §2 the traditional Eulerian approach to estimating the orientation distribution function involves solving a Fokker-Plank equation in conjunction with the Navier-Stokes equations. The utility of this approach is limited as the rotary dif- fusion coefficient, Dr must be determined empirically. We found two expressions in the literature, namely the Folgar-Tucker approach, as given in Equation 2.11 or the more complicated expression, that is Equation 2.12 derived by Koch. In this chapter we examine the usefulness of each expression by comparison with experimental data. In this chapter we compare data obtained along the central streamline. This repre- sents a limiting case of the model equations and as a result, is not computationally expensive to do so. This chapter is divided into a number of subsections. In §4.1 we re-derive the model equations and apply them to the central streamline. In §4.2, the experimental proce- dure is described along with the numerical details used to estimate Dr. The results are given in §4.3. 26 Chapter 4. Estimating Dr 4.1 Model Description The two general methods for modelling the behavior of fibres in suspension are the Lagrangian and Eulerian methods. In the Lagrangian method, the equations of motion for a single fibre are solved given a known velocity field and the process is repeated for each fibre within the suspension. While the Lagrangian method can be quite accurate, it is computationally intensive, particularly for non-dilute suspen- sions flowing within complex geometries. With the Eulerian method, the probability distribution of fibre orientation (and possibly position) is determined by solving a convection-diffusion equation, namely a Fokker-Plank type equation. With this ap- proach the mean flow field convects the fibres position and orientation while long range, hydrodynamic fibre-fibre interactions diffuses fibre orientation. This diffusion type process effectively results in a flux which opposes the local mean orientation state of the fibres. The following assumptions will be made for the Eulerian analysis outlined below: 1. The nominal fibre length is small enough that the flow field can be assumed to be linear along the length of the fibre. 2. The orientation state of the fibres at a point in space can be described by a probability distribution function, Ψ, where Ψ is defined such that the proba- bility of finding a fibre oriented between the angles φ and φ + ∂φ is Ψ(φ)∂φ. 3. The fibres are uniformly distributed in space throughout the flow field. 4. Fibre inertia and buoyancy are negligible. 27 Chapter 4. Estimating Dr 5. Fibres are rigid and large enough that Brownian motion can be neglected. In order to determine the orientation distribution function, Ψ, a steady-state Fokker- Plank equation is solved, that is, the steady-state form of Equation 2.1 u ∂Ψ ∂x +v ∂Ψ ∂y = 1 sinθ ∂ ∂θ (Drsinθ ∂Ψ ∂θ −sinθθ̇Ψ)+ 1 sinθ ∂ ∂φ (Dr 1 sinθ ∂Ψ ∂φ −sinθφ̇Ψ) (4.1) where u and v are the fluid velocity components in the x and y directions respectively. φ and θ are the azimuth and zenith orientation angles respectively (see Figure 4.1) and φ̇ and θ̇ are the rotational velocities for the 3D orientation vector, as given in Equations 2.9 and 2.10 respectively. Dr is the rotary diffusion coefficient which is given by either Equations 2.11 or 2.12. It should be noted that in this work we deal only with the first term from Equation 2.12, that is the isotropic component of the rotary diffusion coefficient. For applications to accelerating flows the second term in Equation 2.12 provides a lower order correction to the isotropic part of the diffusion coefficient. Therefore the marginal improvement in accuracy of the rotary diffusion coefficient does not justify the added computational cost of working with this term. Hence we define Koch’s rotary diffusion coefficient as follows Dr = nL3 α̇ ln2 r (λIE : 〈pppp〉 : E) (4.2) A simplification of the model equations can be realized by considering the flow of fibres along the central streamline only. A detailed derivation of this can be found 28 Chapter 4. Estimating Dr Figure 4.1: The orientation of a fibre with respect to flow in the contraction. The angle φ is the angle of the projection of the fibre in the xy-plane and θ is the angle between the fibre and the z-axis. The angle γ is the angle of the projection in the xz-plane in Olson et al (2004) and the references contained therein. If on the central streamline u(x) = u(0) 1− (1− 1 R ) ( x L ) , v(x) = 0 (4.3) where R is the contraction ratio of the headbox defined as the ratio of inlet to exit heights which for this analysis is fixed at a value of 10, Equation 2.1 reduces to 29 Chapter 4. Estimating Dr u ∂Ψ ∂x = 1 sinθ ∂ ∂θ (Drsinθ ∂Ψ ∂θ − sinθθ̇Ψ) + 1 sinθ ∂ ∂φ (Dr 1 sinθ ∂Ψ ∂φ − sinθφ̇Ψ) (4.4) Here the projected rotational velocities along the central streamline, that is φ̇ and θ̇ reduce to φ̇ = ∂u ∂x sin(2φ) (4.5) and θ̇ = ∂u ∂x sin(φ)cos(φ)cos2(θ) (4.6) The planar orientation distribution function, Ψp can then be extracted from the 3D orientation distribution function as follows. Ψp = ∫ π 0 Ψ(φ, θ)sin(θ)dθ (4.7) The boundary conditions on Ψ are as follows. Since one end of a fibre is indistin- guishable from the other end, Ψ must be periodic, that is Ψ(x, φ, θ) = Ψ(x, φ + π, π − θ) (4.8) 30 Chapter 4. Estimating Dr Ψ must also satisfy a normalization constraint ∫ π 2 −π 2 ∫ π 0 Ψ(x, φ, θ)sin(θ)dθdφ = 1 (4.9) At the inlet, we set Ψ equal to the experimentally observed inlet orientation distri- bution, that is Ψ(0, φ, θ) = Ψexp(0, φ, θ) (4.10) Equations 4.3 - 4.10 represent the system of equations which will be solved in order to estimate Ψ along the central streamline. 4.2 Experimental Detail The main objective of the experimental work is to measure the orientation distribu- tion functions along the central streamline. This data will be used in conjunction with Equations 4.3 - 4.10 to estimate the parameters CI and λ in Equations 2.11 and 4.2. A description of both the experimental device and the numerical procedure are given below. When dealing with non-dilute fibre suspensions it becomes increasingly difficult to distinguish individual fibres from one another, particularly as the fibre concentration increases. One way to deal with this problem is to implement an index of refraction matched system (e.g. Stover, Koch, & Cohen (1992), Yasuda et al (2004)), that is, a suspension is sought which consists of fluid and fibre phases with identical, or nearly 31 Chapter 4. Estimating Dr identical indexes of refraction. The result is that the fibre phase becomes indistin- guishable from the fluid phase when observed under white light; loosely speaking, the fibres become invisible. A small number of fibres (less than 1% of the total num- ber of fibres) are then coloured so that they can be visualized and their motion and orientation are observed and taken to represent the behaviour of all fibres within the suspension. For these experiments, cylindrical boroscilicate glass fibres (www.mosci.com) of di- mensions 5 × 0.1mm were employed as the fibre phase. The glass fibres have a density of approximately 2250 kg/m3. The index of refraction of these fibres was measured commercially and found to be 1.4719(±0.0005). Approximately 0.01% of the fibres in each suspension were silvered using Tollen’s solution, a mixture of 5 ml of 0.1 M AgNO3 with ∼ 10μL (5 drops) of 0.4 M NaOH . Before silvering, the fibres were washed in detergent, rinsed in alcohol and then in distilled water. After silvering the fibres were again washed in a similar manner to remove any loosely adsorbed AgNO3. We measure the motion and orientation of 4 different monodis- persed suspensions of concentration, nL3 = 8, 16, 24 and 32. The Newtonian fluid used in this was glycerine with a density of 1260 kg/m3 and a viscosity of 1.49 Pa s. Fibre settling was not observed over the time-scale of the experiment. The experiments were performed in a rectangular cross-sectional Plexiglas cell (50.8× 130mm) with approximately 1.5 litres of suspension. The contraction consists of two rigid sloping walls, separated by 50.8mm at the entry and 5.08mm at the exit; the contraction ratio is R = 10. The distance between the entry and the exit is 130mm, see Figure 4.2. 32 Chapter 4. Estimating Dr Figure 4.2: An image of the experimental test section. Up- and downstream of the channel are large reservoirs set at different heights to control the pressure drop over the cell. The flowrate was set at 0.0021m3/min. Based on the initial height of the channel, the inlet Reynolds number is estimated to be 0.5. During testing the suspension was recirculated between the two reservoirs with a Masterflex B/T variable speed peristaltic pump (www.coleparmer.com). The suspension was gently mixed prior to entering the channel in an effort to create a more uniform orientation distribution at the inlet while not disturbing the flow in the actual contracting section. During testing the suspension was recirculated between the two reservoirs. The piping for the recirculation loop was centred at the inlet and exit sides of the channel. The reason for this was to promote fibre trajectories through the middle of the channel and minimize any influence of the solid walls 33 Chapter 4. Estimating Dr on fibre orientation. Furthermore, the width of the channel is large relative to the length of the fibre so we neglect any influences of the wall boundary layers on the fibre orientation distribution. Our visualization system consiste5d of a progressive scan Basler A201b monochrome CCD camera (10 bit grey scale and 1008 × 1016 pixel spatial resolution with a framing rate of 5 frames per second), mounted with an F-mount Micro-NIKKOR 105 mm lens. A first surface mirror (Edmund Optics) was positioned at a 45◦ angle above the contraction in order to capture images from both the xy- and xz-planes simultaneously. The data obatained from observations of fibre orientation in the xz-plane were used only to initialize the orientation distribution function at the contraction inlet when measuring the diffusion coefficient. With the optical setup described the imaged volume is 130×65×5.1mm3 (length× width× elevation). For particle tracking, the lens aperture and focus, backlight intensity and camera exposure time were chosen so that the whole imaged volume was within the depth of field of the lens. The plexiglass cell was transilluminated using Schott-Fostec fibre optic dual backlight. The images were split between the two projection planes after which a mask was applied to isolate the actual flow region from the remainder of the image. The images were then saved using a Matrox Meteor II digital frame grabber, and the fibre tracking analysis was performed on the xy-plane. The orientation of each particle was calculated using an in-house fibre tracking algorithm written with Matlab v7.2 and consists of the following subroutines: 1. Particle identification - Particle edges are detected using a Prewitt edge de- 34 Chapter 4. Estimating Dr xy−plane xz−plane xy−plane xz−plane Figure 4.3: Two consecutive unprocessed images of fibre observation separated by 200ms. The upper image pair sequentially precedes the bottom image pair. 35 Chapter 4. Estimating Dr Figure 4.4: An example of an original image (left) and a post-processed image with white fibres on a black background (right). The fibre concentration is nL3 = 8, the fibre length is 5mm, and the aspect ratio is 50. 36 Chapter 4. Estimating Dr tection scheme with the original greyscale image converted to a binary image. Particle edge pixels are set to white (pixel value of 1) on a black background (pixel value 0). 2. Particles found in the image are dilated, the interior regions filled with white pixels, and then eroded back to their original size. The image now contains filled white particles on a black background (see Figure 4.4). 3. Any particle whose length and width are less than 10 pixels in size, or whose length and width are approximately equal are assumed to be artifacts in the image (e.g. air bubbles or non-fibrous debris) and are removed from the image. 4. Orientation Distribution - The orientation angle of each fibre is measured rel- ative to the horizontal by computing the arctangent of the ratio Δy Δx between the end points of a fibre. The associated fibre length and the position of the center of area is also recorded. 5. The flow region is partitioned into 5mm×5mm cells and the orientation angle of each fibre whose center of area lies within a particular cell is computed. In total we fit 26 cells in the horizontal direction along the 130mm length contraction centerline. There were on average approximately 5000 tracer fibre observations in each cell for each experiment. 37 Chapter 4. Estimating Dr 4.2.1 Numerical Estimates of Dr The diffusion coefficient is determined from observations of fibre orientation in a region near the central streamline. More specifically, the region of analysis for these particular experiments is defined by the area bounded by the straight streamlines that are exactly 30% of the distance from the centerline to the channel walls both above and below the centerline (see Figure 4.5). All fibre observations within this region are then binned into their corresponding spatial cells along the x-axis according to the position of the center of area of each observed tracer fibre (see experimental detail section for full description of binning tracer fibres into spatial cells). To reiterate, there are 26 cells along the central streamline, with a length of 5mm each. Analysis Region Figure 4.5: The region bounded by the red lines represents the analysis region. All fibres observed in this area are used for the 1D analysis and determination of CI and λ. In order to characterize the results, a method was developed to estimate Dr from these experimental measurements. An inverse solver was used to estimate the empir- ical parameters of the constitutive relationships, that is either Equation 2.11 or 4.2, by comparing the experimental data to Equation 4.4. The inverse solver consisted 38 Chapter 4. Estimating Dr of the following components. 1. Equation 2.1 was discretized using a second-order accurate centered difference scheme in both orientation angles φ and θ and a first-order accurate forward difference scheme in x. In this case the velocity field along the central streamline is assumed to be of the form u = (u(x), 0) with the boundary conditions given by Ψ(x, φ, θ) = Ψ(x, φ + π, π − θ) (4.11)∫ π 2 −π 2 ∫ π 0 Ψ(x, φ, θ)sin(θ)dθdφ = 1 (4.12) At the inlet, that is at x = 0, the orientation distribution is set equal to the orientation distribution observed experimentally. 2. We initialize the numerical procedure by assuming a trial value for CI or λ depending upon which rotary diffusion model is being used. 3. Equation 4.4 is solved using a Gauss-Seidel method. 4. When solving using Koch’s definition, the value of Dr is updated through 4.2 using the new estimate of Ψ. 5. Steps 3 and 4 are repeated until the L2-norm of the relative change in the solutions for Ψ and Dr between successive iterations diminish to a value smaller than 10−6. 6. A golden search routine is employed to minimize the L2-norm of the difference between the model and experimental prediction. 39 Chapter 4. Estimating Dr 4.3 Results To help understand the main findings of this section, it is instructive to discuss qual- itatively the evolution of Ψ along the length of the channel. We do so by displaying one representative result in Figure 4.6. Here we see the probability distribution at four different axial positions. Each distribution is centered around zero, that is they are aligned in the direction of the flow, and the spread in the distribution de- creases with increasing axial positions. These results are representative of all cases tested. −1.5 −1 −0.5 0 0.5 1 1.5 0 1 2 3 4 5 6 7 nL3 = 8 φ Ψ Figure 4.6: Development of the orientation distribution along the contraction center- line for the φ projection. The experimental data are given at x L = 0(◦), x L = 0.33(), x L = 0.66() and x L = 1.0(). The concentration shown here is nL3 = 8 and the solution shown here was obtained with Equation 4.2. To fully characterize the evolution of the orientation distribution we display the variance σ2 of the orientation distributions for all concentrations tested, see Figure 4.7. To help clarify this figure, a low value of σ2 corresponds to highly aligned orientation 40 Chapter 4. Estimating Dr 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 x σ2 Figure 4.7: Plot of the orientation variance σ2 along the contraction length for nL3 = 8(), nL3 = 16(), nL3 = 24(◦), and nL3 = 32(). The uncertainty of the data are represented by the error bars which have been computed at the 95% confidence interval. state, and a high value of σ2 corresponds to a more random, or isotropic orientation state. The first observation that can be made from Figure 4.7 is that the difference in σ2 at the outlet is small, but nonetheless unique for each concentration. Con- versely, at the inlet we see that the variance is significantly higher for the highest concentration. This implies that the orientation distribution of the suspension at the inlet of the channel at nL3 = 32 is more random than the other concentrations tested. Another interesting observation is that the highest concentration, that is, nL3 = 32, enters the contraction with the greatest degree of orientation isotropy, yet is the most aligned at the exit. We now turn to the main findings of this section where CI and λ are determined as a 41 Chapter 4. Estimating Dr function of fibre concentration. To reiterate, an optimization routine is performed to estimate CI and λ based on experimental data measured just prior to the contraction exit. The reason for measuring these parameters near the exit is based on the fact that the suspension becomes highly aligned in this region. As a result, the second order terms in Equation 4.4 play their most significant role in this region for fibre orientations parallel to the central axis (that is near φ = 0 and θ = π 2 ). Hence the accuracy of CI and λ are most critical near the exit and good agreement is still maintained throughout earlier stages of the contraction. Estimates using Equations 2.11 and 4.2 are presented in Figure 4.8. 5 10 15 20 25 30 35 3 4 5 x 10 −3 C I nL3 0 0.005 0.01 λ Figure 4.8: Plot of CI(−−) and λ(−−) vs nL3. As shown in this figure CI initially increases monotonically with concentration. This result was expected given that higher concentrations lead to an increase in the fre- quency of fibre-fibre interactions, hence greater diffusion of orientation. However, we see that for the highest concentration, nL3 = 32, there was a dramatic decrease in the 42 Chapter 4. Estimating Dr interaction coefficient. It is believed that at such a high concentration the orienta- tion model breaks down due to the assumption of purely hydrodynamic interactions. More specifically, it is believed that at higher concentrations the effects of fibre floc- culation and mechanical contact interactions between fibres has a hindering effect on fibre orientation. This observation further advances the argument that flocculation or mechanical entanglement is the cause for such a condition. Martinez et al (2001) have demonstrated that for papermaking fibres, fibre mobility diminishes through flocculation at nL3 = 31. This phenomenon is compounded with the aligning effects of the accelerating flow inside the contraction by which the suspension is forced into a highly aligned state and fibre mobility disappears almost completely. At this point it is worth comparing the above measurements of CI with those values published in the literature. For example, Folgar & Tucker (1984) measured values of CI in the range 0.003−0.0165 for fibre volume fractions, c, in the range 0.0004−0.16 and with a constant fibre aspect ratio of 16. This is comparable to our measurements, that is, CI in the range 0.0036− 0.0045 for volume fractions of 0.0025− 0.01 and a constant aspect ratio of 50. Of course the main difference is that our lowest value of CI corresponds to the highest concentration. Another probable cause for the differences, in addition to the fact that we have used a different flow geometry, is that Folgar & Tucker (1984) have based their measurements on a planar orientation model whereas the measurements made in this work are based on a more realistic, 3D orientation model. This was a mistake on the part of Folgar & Tucker (1984). It is also interesting to compare our results to those obtained by Phan-Thien et al. (2002) through direct simulation of fibre suspensions. They report simulated values 43 Chapter 4. Estimating Dr 5 10 15 20 25 30 35 3 4 5 x 10 −3 C I nL3 2 3 4 x 10 −3 λ n L3 /ln 2 r Figure 4.9: Shown is a comparison of λ nL 3 ln2r (−−) with CI(−−) vs nL3. of CI in the range 10 −2−10−3, which is again comparable to our measurements. One interesting aspect of their studies is that they base their correlations for CI on the parameter c · r as opposed to nL3, where r is the fibre aspect ratio, that is the ratio of length to diameter. The trend for λ is also interesting. We see that λ decreases monotonically with increasing fibre concentration, something which is unexpected since the effects of fibre concentration should be accounted for in the interaction rate term in the diffusion coefficient, that is, the nL3 part of Equation 4.2. This finding would suggest that the concentration dependence in Equation 4.2 is not quite correct, at least for a linearly contracting channel. What is interesting, however, is to plot the λ nL 3 ln2r term in Equation 4.2, alongside CI . This is shown in Figure 4.9. This plot shows several interesting points. The first observation is that the diffusion coefficient does scale with the velocity gradient ‖α̇‖ for this one-dimensional flow. 44 Chapter 4. Estimating Dr The second key observation is that the term in Equation 4.2 containing the fourth order orientation tensor effectively acts as a amplification factor to the diffusion coefficient for each concentration. 4.3.1 Comparison of Rotary Diffusion Models The evolution of Ψ is now characterized for the experimental conditions tested. An example of the development of both the predicted and observed orientation distribu- tion along the contraction centerline is shown in Figure 4.10 for a few axial positions. Included in these Figures is a comparison of the experimental data with the two dif- ferent equations for the rotary diffusion coefficient. The histograms of the orientation probability functions for each axial position were found to be smooth and clustered around a definite value near zero. This behavior is typical of all axial positions. In order to quantitatively evaluate the performance of the two models, a χ2 goodness of fit test is performed whereby each model is compared to the 1D experimentally determined orientation distribution function at the contraction exit. These results are summarized in Table 4.1 where a lower χ2 value indicates better agreement with the experimental distribution. The χ2 values shown in Table 4.1 indicate that Koch’s model better predicts the orientation state at the exit, however a visual inspection of Figure 4.10 suggests that the improvement over the simple Folgar-Tucker model is marginal. It should be emphasized that the computational requirements for solving Equation 4.4 using Equation 4.2 are far greater than those required to solve Equation 4.4 using Equation 2.11. For example, a grid size consisting of 50 x 200 x 200 nodes 45 Chapter 4. Estimating Dr −1.5 −1 −0.5 0 0.5 1 1.5 0 1 2 3 4 5 6 7 nL3 = 8 φ Ψ −1.5 −1 −0.5 0 0.5 1 1.5 0 1 2 3 4 5 6 7 φ Ψ nL3 = 8 −1.5 −1 −0.5 0 0.5 1 1.5 0 1 2 3 4 5 6 7 nL3 = 16 φ Ψ −1.5 −1 −0.5 0 0.5 1 1.5 0 1 2 3 4 5 6 7 φ Ψ nL3 = 16 −1.5 −1 −0.5 0 0.5 1 1.5 0 1 2 3 4 5 6 7 nL3 = 24 φ Ψ −1.5 −1 −0.5 0 0.5 1 1.5 0 1 2 3 4 5 6 7 φ Ψ nL3 = 24 −1.5 −1 −0.5 0 0.5 1 1.5 0 1 2 3 4 5 6 7 nL3 = 32 φ Ψ −1.5 −1 −0.5 0 0.5 1 1.5 0 1 2 3 4 5 6 7 φ Ψ nL3 = 32 Figure 4.10: Development of the orientation distribution along the contraction centerline for the φ projection. The experimental data are given at x L = 0(◦), x L = 0.33(), x L = 0.66() and x L = 1.0(). The best fit of Equation 4.4 is given as the solid lines and the predictions given as the hatched lines. The left column is the solution using Equation 4.2, and the right column is the solution using Equation 2.11. 46 Chapter 4. Estimating Dr for the x, φ, θ variables was required when using Equation 4.2, where as a grid size consisting of 50 x 100 x 100 nodes for the x, φ, θ variables was suitable when using Equation 2.11. In addition to the larger grid size size, several iterations (typically 4 or 5) were required when using Equation 4.2 since Dr had to be updated after each new solution was obtained to Equation 4.4. Table 4.1: χ2 Values for Each Model nL3 Koch Folgar-Tucker 8 29.5 63.5 16 29.5 63.5 24 28.8 62.1 32 28.8 62.1 4.4 Summary This chapter has presented a description of the model for predicting the orientation state of a fibre suspension flowing through a linearly contracting channel. A sim- plified model has also been presented for predicting the orientation state along the central streamline of the contraction. The experimental methods used to measure the orientation state of the suspension, and to measure the rotary diffusion coeffi- cient were given. Numerical estimates for the unknown closure parameters in the two models for the rotary diffusion coefficient have been presented, and the utility of the two models evaluated by comparison with experiments. It was shown that the interaction coefficient, that is, CI in Equation 2.11 first increases linearly with con- centration up to a point, after which it suddenly decreases when the concentration 47 Chapter 4. Estimating Dr is further increased. It is argued that this sudden decrease in CI indicates the point at which the model for the rotary diffusion coefficient breaks down and that this breaking point is the result of mechanical interactions and fibre flocculation. The relationship between λ in Equation 4.2 and concentration was also presented. These results show that λ does not remain constant but rather decreases monotonically as the concentration is increased. This finding would suggest that the scaling term in 4.2 is not quite correct for this particular type of suspension flow. A comparison of the two models for the rotary diffusion coefficient with experiments have shown that Koch’s model is better able to predict the orientation distribution at the contraction exit and it is believed that the improvement comes from the fourth order orientation tensor. The marginal improvement in accuracy gained through the use of Equation 4.2 comes at an extremely high computational cost, one which perhaps does not justify its usefulness. 48 Chapter 5 Validating Dr Estimates for the rotary diffusion coefficient were made in §4. These estimates were based on observations of fibre orientation along the central streamline of the contrac- tion. In this section, we use these estimates for Dr and CI to predict fibre orientation over the entire xy-plane of the contraction and include the two-way coupling between the fluid and fibre phases. We do so for concentrations of nL3 = 8, 16 and 24. The last data point, that is nL3 = 32 is not included in this analysis since it has been shown that the fibre orientation model breaks down at such a high concentration due to mechanical interactions and flocculation effects. This section includes a review of the model for predicting the planar orientation distribution over the xy-plane, the constitutive equations for the two way coupling, and the numerical formulation for computing Ψ. This is followed by a comparison of the predicted 2D orientation distribution with experiments and a discussion of the similarities and differences. In §5.4 the effects of the two way coupling on the resulting flow field and orientation distribution are investigated. 49 Chapter 5. Validating Dr 5.1 Model Formulation The assumptions made in this analysis are identical to those made in the 1D analysis. In order to predict the planar orientation distribution function, Ψ, a steady-state, 2D Fokker-Plank equation is solved, that is u ∂Ψ ∂x + v ∂Ψ ∂y = Dr ∂2Ψ ∂φ2 + ∂(φ̇Ψ) ∂φ (5.1) The orientation velocity, φ̇, was shown in §2 to be independent of θ and therefore assumes that same form as that given by Equation 2.9, that is φ̇ = 1 2 ( ∂v ∂y − ∂u ∂x ) sin(2φ)− ∂u ∂y sin2(φ) + ∂v ∂x cos2(φ) (5.2) The boundary conditions for Ψ are as follows: 1. Periodic boundary conditions are enforced with respect to the variable φ which results in the condition Ψ(x, y, φ) = Ψ(x, y, φ + π) (5.3) 2. Ψ must satisfy a normalization constraint ∫ π 2 −π 2 Ψ(x, y, φ)dφ = 1 (5.4) 50 Chapter 5. Validating Dr 3. At the contraction inlet, Ψ is set to equal to the experimentally observed ori- entation distribution, that is Ψ(0, y, φ) = Ψexp(0, y, φ) (5.5) The fluid flow is described using Cauchy’s equations of motion for an incompressible fluid, that is ∇ · u = 0 (5.6) ρ ( ∂u ∂t + u · ∇u ) = −∇P +∇ · τ (5.7) where ρ is the fluid density, P is the pressure and τ is the stress tensor, that is the sum of both the Newtonian fluid and fibre contributions τ = μ(∇u+∇uT ) + τ fibre (5.8) where τ fibre is given by Equation 2.15. The magnitude of the velocity distribution at the inlet is given a value identical to that used in the experiments. In order to simplify the problem and reduce computational time and space, Equation 5.1 is solved along individual streamlines of the flow under steady-state conditions. This reduces Equation 5.1 and 5.2 to a quasi-1D problem which is solved as fol- 51 Chapter 5. Validating Dr lows. 1. For each streamline, the values of u, v, ∂u ∂x , ∂v ∂x , ∂u ∂y , ∂v ∂y are used to directly com- pute φ̇(x, y) according to Equation 2.9 at each streamline position, s(x, y). 2. Equation 5.1 is then mapped onto the streamlines according to the following transformations ∂Ψ ∂x = ∂Ψ ∂s ∂s ∂x (5.9) ∂Ψ ∂y = ∂Ψ ∂s ∂s ∂y (5.10) This leads to the steady-state, quasi 1D form of 5.1 u ∂Ψ ∂s ∂s ∂x + v ∂Ψ ∂s ∂s ∂y = Dr ∂2Ψ ∂φ2 + ∂(φ̇(s(x, y))Ψ) ∂φ (5.11) 3. Equation 5.11 is discretized using forward differences with respect to the spatial streamline variable, s, and centered differences with respect to the orientation angle φ. 4. The flow field is computed in this device using a 2D commercial software pack- age (FLUENT) with the gradient of the additional fibre stress term treated as a momentum source. An iterative procedure is used whereby the flow field is initially determined for the pure fluid alone (i.e. water with no fibres present), 52 Chapter 5. Validating Dr after which Equations 5.1 and 5.2 are solved using the initial flow field data. The gradient of the fibre stress is then computed based on the initial solution for Ψ, which is then fed back into the fluid flow equations, which are then solved again to produce a new flow field, which is then used again in solving equation 5.1 for a new update of Ψ. 5. The process is repeated until the L2-norm of the change in the flow field between successive iterations is smaller than a critical value, 10−6. A total of four iterations were required for all concentrations considered. An example of the convergence of the flow field is shown in Figure 5.1 below and is typical of all concentrations. 1 1.5 2 2.5 3 3.5 4 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 iteration ||U i − U i− 1| | nL3 = 8 Figure 5.1: Plot of the L2 norm of the relative change in the flow field between successive iterations. The concentration shown here is nL3 = 8. Equation 5.1 was solved over the contraction using a mesh consisting of 736 spatial nodes and uniform grid in φ consisting of 500 nodes. When mapped onto streamlines, there were typically 60 spatial nodes along each streamline. Mesh independent solu- 53 Chapter 5. Validating Dr tions were obtained for Ψ. Figure 5.2 shows an example of the mesh independence whereby the L2 norm of the solution difference over the contraction is given for solu- tions obtained with the 736 x 500 mesh and another solution obtained with a 1397 x 1000 mesh. The difference in mesh solutions over the channel is seen to be small with the greatest difference occurring near the inlet and away from the central streamline region. The greatest difference occurs in a small region around the lower inlet of the contraction domain. This is a consequence of a very sharply peaked experimentally observed orientation distribution at this part of the inlet which is difficult to smooth with the numerical solver. This difference does considerably diminish further down the channel and is not believed to significantly alter the numerical predictions. 0.28 0.26 0.24 0.22 0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 Figure 5.2: Plot of the L2 norm of the relative change in solutions along the central streamline between solutions obtained with mesh sizes of 736 x 500 and 1397 x 1000. 5.2 2D Orientation Results The 2D orientation state results are presented below. This includes both the numer- ical predictions and the experimental observations over the xy-plane. 54 Chapter 5. Validating Dr 5.2.1 Theoretical Predictions In order to better visualize and interpret the orientation distribution over the en- tire 2D contraction plane we use the mean, φ̄, and variance, σ2, of the orientation distribution function. The mean value of orientation corresponds to the location of the peak in the orientation distribution, while the variance can be interpreted as a measure of the degree of orientation isotropy. That is, a small value of orientation variance corresponds to a suspension that is highly aligned in the mean direction and vice-versa. Each are defined in the traditional manner as follows φ̄ = ∫ π 2 −π 2 φΨ(x, y, φ)dφ (5.12) σ2 = ∫ π 2 −π 2 (φ− φ̄)2Ψ(x, y, φ)dφ (5.13) Figure 5.3 shows the local mean orientation as contour plots over the xy-plane. The mean value at the inlet is zero, largely a consequence of the experimental inlet condition. Away from the inlet, the mean orientation value tends to coincide with the streamline direction throughout the remainder of the contraction, and is essentially independent of concentration. The negligible effect of fibre concentration on the mean orientation angle is somewhat expected since fibre-fibre interactions are assumed to randomize fibre orientation but should not affect the direction of mean orientation. There is also a small region near the solid boundaries of the contraction, which shall be referred to as the wall layer region, where the mean fibre orientation angle exceeds 55 Chapter 5. Validating Dr 0.18 rad while the contraction angle is approximately 0.17 rad. This result suggests that near the walls, the model predicts a shift in fibre orientation away from the wall. 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1 -0.12 -0.14 -0.16 -0.18 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1 -0.12 -0.14 -0.16 -0.18 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1 -0.12 -0.14 -0.16 -0.18 Figure 5.3: Numerical predictions of the mean fibre orientation angle. Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom). Figure 5.4 shows contour plots of the model predictions of orientation variance, σ2, over the xy-plane. The variance is seen to be at its largest at the inlet, which is again a consequence of the inlet boundary conditions, and tends to decrease towards the contraction exit. This reflects the fact that fibres align in response to the large flow gradients near the contraction exit. Near the walls however we see an increase in 56 Chapter 5. Validating Dr the isotropy of fibre orientation, that is a high value of variance, suggesting that the walls effectively randomize fibre orientation. It is also interesting to note the large development region near the inlet where flow gradients are much smaller. It can be interpreted that in this region, fibre-fibre interactions have a much more dominant effect and help maintain the orientation isotropy of the suspension. The effect of increasing fibre concentration is reflected in the increased size of the development region. This is a consequence of the fact that increasing the suspension concentration leads to an increase in the frequency of fibre-fibre interactions, and ultimately a greater diffusion of orientation and a more isotropic orientation state. This effect is most obvious in the earliest stages of the contraction where the flow gradients are still relatively small. However, the effects of the fibre-fibre interactions completely disappear towards the exit of the contraction. 5.2.2 Experimental Observations The procedure for statistical analysis of the 2D data was performed in an analogous manner to the 1D analysis. However, with the 2D data set we measure the mean fibre orientation, orientation variance, and orientation distribution at a set of 5mm x 5mm spatial grid cells filling the entire section of the channel. Figure 5.5 shows the mean fibre orientation observed experimentally. The first ob- servation is that the suspension is oriented in the positive direction in the lower half of the contraction and in the negative direction in the upper half of the contraction. This supports the model predictions that fibres tend to align in the streamline di- 57 Chapter 5. Validating Dr 0.25 0.24 0.22 0.21 0.20 0.18 0.17 0.16 0.14 0.13 0.12 0.10 0.09 0.08 0.06 0.05 0.04 0.02 0.01 0.25 0.24 0.22 0.21 0.20 0.18 0.17 0.16 0.14 0.13 0.12 0.10 0.09 0.08 0.06 0.05 0.04 0.02 0.01 0.25 0.24 0.22 0.21 0.20 0.18 0.17 0.16 0.14 0.13 0.12 0.10 0.09 0.08 0.06 0.05 0.04 0.02 0.01 Figure 5.4: Numerical predictions of the orientation variance, σ2. Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom). rection. It can also be seen that the mean orientation approaches the zero direction along upper and lower streamlines near the contraction exit, an observation that was also predicted numerically. While, there is good qualitative agreement between the experimental results and the model predictions away from the walls, it is clear that the model has severely overpredicted fibre orientation in the wall layer region. That is, the experiments show that fibres align with the wall in the wall layer whereas the model predicts a shift in fibre orientation away from the walls in this region. Figure 5.6 shows the experimental observations of orientation variance over the con- 58 Chapter 5. Validating Dr 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1 -0.12 -0.14 -0.16 -0.18 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1 -0.12 -0.14 -0.16 -0.18 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1 -0.12 -0.14 -0.16 -0.18 Figure 5.5: Experimental observations of the mean fibre orientation angle. Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom) traction. The trend for σ2 away from the walls is very similar to the model predic- tions. That is, we see the greatest values of σ2 near and around the inlet, and this effect propagates a finite distance down the channel where the flow gradients are lower. Near the exit, σ2 becomes very small indicating a highly aligned suspension. The effect of concentration is also similar to the model predictions, that is, increas- ing the fibre concentration increases the size of the development region. Near the walls however, the experiments differ considerably from the model predictions. Near the walls the experiments reveal a very low value of σ2, indicating a highly aligned 59 Chapter 5. Validating Dr orientation state near the wall. This contradicts the model predictions of an increase in orientation isotropy near the walls. 0.25 0.24 0.22 0.21 0.20 0.18 0.17 0.16 0.14 0.13 0.12 0.10 0.09 0.08 0.06 0.05 0.04 0.02 0.01 0.25 0.24 0.22 0.21 0.20 0.18 0.17 0.16 0.14 0.13 0.12 0.10 0.09 0.08 0.06 0.05 0.04 0.02 0.01 0.25 0.24 0.22 0.21 0.20 0.18 0.17 0.16 0.14 0.13 0.12 0.10 0.09 0.08 0.06 0.05 0.04 0.02 0.01 Figure 5.6: Experimental observation of the orientation variance σ2. Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom) 5.3 Discussion of Fibre Orientation Results As discussed in Chapter 2, the near wall behavior of fibres has received some attention over the years. For example, Hyensjö et al. (2007) modeled the influence of shear on fibre orientation in turbulent headbox flows. They showed that near the walls, the 60 Chapter 5. Validating Dr tendency is for the mean fibre orientation angle to shift away from the wall and for the orientation variance to increase. Their numerical predictions compared well, at least qualitatively with the experimental results of Asplund & Norman (2004). In the latter studies it was shown that the degree of fibre alignment is at its lowest very close to the wall and then increases towards the middle of the channel where it remains approximately constant with distance from the wall. That is, fibre alignment, or orientation anisotropy as they called it, decreases significantly in a small region near a solid boundary, suggesting the idea that solid boundaries have a randomizing effect on fibre orientation. In studies by Hämäläinen & Hämäläinen (2007), fibre orientation was numerically predicted for a turbulent suspension in a 2D linear contraction. They prescribed conditions on Ψ at the wall such that the mean orientation angle was set equal to the wall angle however no physical justification for this choice of boundary conditions were given. Studies of near wall fibre orientation by Stover, & Cohen (1990) showed that fibres may orient away from walls in laminar flow, but only under certain conditions. They found that when fibres with a low Jeffery orbit constant, that is highly aligned fibres, came within a half fibre length from the wall the fibre remained close to, and aligned with the wall indefinitely. But fibres with a larger period of rotation (approximately 50% larger) shifted away from the wall. The experimental observations in this work suggest that both possibilities are possible, that is, fibres might align in the wall direction or they might reorient away from the wall. It should be pointed out that 5mm length fibres were used in this study which is roughly the same length scale as the wall layer region. It is possible that for fibre lengths comparable to the length scale of the wall layer, fibres will align with the 61 Chapter 5. Validating Dr wall, where as smaller fibres might not. This conjecture was not invetigated in this work. Other more subtle differences between the model predictions and the experimental observations, such as differences in size of the development region for example, may partly be the result of comparing a planar orientation model, that is Equation 5.1 to the projection of 3D experimental fibre orientations onto the xy-plane. It might be interesting to solve the full 3D fibre orientation problem, that is, Equation 2.1, extract the planar orientation distribution functions, Ψp, as was done in the 1D results section. Of course, the difficulty with this method is that it is quite computationally intensive, particularly near the solid boundaries where the flow gradients become large in non-streamwise directions. 5.4 Effect of 2 Way Coupling In this section we discuss the effects of the two way coupling between the fluid and fibre phases. That is, we investigate the altered structure of the resulting flow field, as well as any changes to the orientation distribution function, Ψ. 5.4.1 Flow Field Figure 5.7 shows contour plots of the magnitude of the two way coupled flow field for the different concentrations under investigation. Shown also is a contour plot of the magnitude of the flow field for the pure fluid, that is, water with no fibres present. 62 Chapter 5. Validating Dr 0.095 0.09 0.085 0.08 0.075 0.07 0.065 0.06 0.055 0.05 0.045 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0.095 0.09 0.085 0.08 0.075 0.07 0.065 0.06 0.055 0.05 0.045 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0.095 0.09 0.085 0.08 0.075 0.07 0.065 0.06 0.055 0.05 0.045 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0.095 0.09 0.085 0.08 0.075 0.07 0.065 0.06 0.055 0.05 0.045 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 Figure 5.7: Contour plots of numerically predicted velocity magnitude. Shown are the pure solvent (Top), nL3 = 8 (Mid-Top), nL3 = 16 (Mid-Bottom) and nL3 = 24 (Bottom). 63 Chapter 5. Validating Dr The general observation is that the effect of the fibre phase is to redirect mass near the walls towards the middle of the contraction, effectively reducing the velocity magnitude along and around the centerline and creating something similar to a plug- type flow condition. The predicted flow field with the fibre phase is considerably different than the flow field without fibres, where the velocity is seen to be at a maximum along the centerline and decreases towards the walls, a phenomenon typical of Newtonian channel flow. However when the fibres are present the contours become somewhat constant across the channel height, that is in the y-direction, except for in a small region very close to the wall, a phenomenon typically observed with slow flows of a yield stress fluid, that is a plug region forms. This observation becomes more and more prevalent as we approach the end of the contraction. In the earlier stages of the contraction, the flow actually becomes slowest along the centerline. The effect of increasing the fibre concentration on this observation is a very subtle one, where the flow becomes marginally slower as the concentration increases. This sort of flow structure has been observed by other researchers when modelling fibre suspension flows. For example, Heath et al. (2007) experimentally studied the flow of a 0.4% mass concentration pulp suspension through a 1:5 axis symetric sudden expansion. They found that for a mean velocity, U ≈ 0.5, the suspension moved as a plug far after the expansion, with the exception of a narrow gap near the wall where flow gradients were found to exist. This is similar to what has been shown in this work when the two way coupling is included. Attempts to model the flow of fibre suspensions using non-Newtonian stress models have seen limited success. For example, Duffy (2003) reviewed many attempts at using non-Newtonian 64 Chapter 5. Validating Dr models such as shear thinning and thickening, Bingham plastic, and friction factor- Reynolds number method and showed that none of these models were suitable when describing fibre suspension flow. Duffy (2003) proposed several mechanistic based models, however these models relied on empirical correlations in order for them to be successful. It would certainly make for interesting future work to confirm the predictions made in this thesis with experimental measurements. 5.4.2 Orientation Distribution Here we look at the effect of the two way coupling on the orientation distribution function. We restrict the analysis to the central streamline only for the sake of clarity, but it should be stated that this case is representative of all other points in the contraction. Shown in Figure 5.8 is a comparison of the orientation distribu- tion functions obtained with and without the two way coupling at the end of the contraction. What is clear from Figure 5.8 is that the effect of the two way coupling on the ori- entation distribution function is very small. The differences are summarized quan- titatively in Table 5.1 below where we give the L2 norm of the difference between coupled and uncoupled solutions at the end of the contraction. Also given is the rel- ative difference between the coupled and uncoupled solutions. For all concentrations the relative change in the orientation distribution is less than 2% when the two way coupling is included. 65 Chapter 5. Validating Dr Table 5.1: Difference Between Coupled and Uncoupled Orientation Distributions nL3 ||ΨCoup −ΨNoCoup||2 ||ΨCoup −ΨNoCoup||2/||ΨNoCoup||2 8 0.3561 0.0135 16 0.3064 0.0119 24 0.1039 0.0045 66 Chapter 5. Validating Dr −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0 1 2 3 4 5 6 7 φ Ψ exit −0.1 −0.05 0 0.05 0.1 3 3.5 4 4.5 5 5.5 6 6.5 φ Ψ exit Coupled Uncoupled −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0 1 2 3 4 5 6 φ Ψ exit −0.1 −0.05 0 0.05 0.1 3 3.5 4 4.5 5 5.5 6 φ Ψ exit Coupled Uncoupled −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 φ Ψ exit −0.1 −0.05 0 0.05 0.1 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 φ Ψ exit Coupled Uncoupled Figure 5.8: Comparison of the orientation distribution functions obtained with and without the two way coupling at the end of the contraction. Shown on the right hand side is a close up near the peak of the distribution functions. Shown are nL3 = 8 (Top), nL3 = 16 (Middle), nL3 = 24 (Bottom). 67 Chapter 6 Summary and Conclusions 6.1 General Summary This thesis has investigated the effects of long range hydrodynamic fibre-fibre in- teractions on the orientation state of a semi-dilute, rigid fibre suspension flowing through a linear contracting channel under laminar conditions. The effects of fibre- fibre interactions have been modeled mathematically, the governing equations solved numerically and the predicted results compared with experimental observations. The orientation state of the suspension has been completely described by an orientation distribution function which evolves according to a Fokker-Plank type equation. The fibre-fibre interactions are assumed to be random in nature and effectively result in a flux which opposes the local mean orientation state of the fibres thus mimicking diffusion-type process. In the first part of this work, the diffusion coefficient was measured along the central streamline of the contraction and the utility of two different models for the rotary diffusion were compared. Each of the two models contains one unknown propor- tionality constant which was determined from the experimental measurements of 68 Chapter 6. Summary and Conclusions the orientation distribution along the channel centerline. This was accomplished by means of an inverse solver where the evolution equation for the orientation distri- bution was solved and the solution fit at the contraction exit for four different fibre concentrations. A comparison was then made between the orientation distributions predicted by each of the two models with those observed experimentally. In the second part of this work, the orientation distribution was solved in two di- mensions over the entire plane of the contraction using the measured values for the diffusion coefficient from part one. A two way momentum coupling between the fluid and fibre phases was assumed and its effect on both the suspension orientation state and structure of the flow field was investigated. 6.2 Diffusion Coefficient The rotary diffusion coefficient was measured along the contraction centerline as a function of fibre concentration. In either case it was found to first increase with increasing suspension concentration up to a maximum, and then decrease with con- centrations above this point. This non-monotonic behavior was attributed to fibre clumping or flocculation, a mechanism not considered in the relationships for the rotary diffusion coefficient. To be more specific, it is believed that at fibre con- centrations above a critical value, nL3 = 24 in this case, the assumption of purely hydrodynamic fibre-fibre interactions is no longer a valid one and that the effects of mechanical contacts between fibres play a more dominant role. The result is that fibre mobility is greatly reduced, and that this reduction in mobility, combined with 69 Chapter 6. Summary and Conclusions the aligning effects of the high flow gradients near the exit of the contraction, force the suspension into a highly aligned state that can not be undone. 6.2.1 Folgar-Tucker Model The unknown proportionality constant in the Folgar-Tucker model, that is, the inter- action coefficient, CI , is determined as a function of fibre concentration. The results showed that CI increases linearly with concentration up to nL 3 = 24 after which it suddenly decreases to a minimum value at nL3 = 32. This behavior indicates that CI is effectively a measure of the diffusion coefficient and it indicates a linearity in rotary diffusion with respect to concentration for the range of validity. 6.2.2 Koch Model The unknown proportionality constant in the Koch model, that is, λ, in Equation 4.2 is determined as a function of fibre concentration. It was shown that λ decreases non- linearly with increasing fibre concentration. This result was unexpected as Equation 4.2 already contains concentration dependency, that is, the nL3 part of Equation 4.2. This might suggest that the rotary diffusion coefficient does not scale with nL3, but perhaps with some other combination of n, L, and r. This finding is unfortunate as it would be highly desirable from a modelling standpoint to be able to determine a proportionality constant once for a particular flow, then be able to reliably predict the effect of varying other parameters without measuring a new proportionality constant. 70 Chapter 6. Summary and Conclusions 6.3 Comparison of Diffusion Models The utility of the two models for the rotary diffusion coefficient were evaluated along the central streamline of the contraction by comparison with experiment. It was shown that both models give good predictions of fibre orientation along the contrac- tion once the diffusion coefficient has been determined. It was shown quantitatively through computation of the χ2 goodness of fit values that Koch’s model for Dr gives a marginal improvement in predicting fibre orientation at the contraction exit, but this improvement comes at a computational cost that perhaps does not justify its use, at least for the case of highly aligned contraction flows as studied in this work. In general, it is concluded that the diffusion coefficient must be determined empir- ically for a given set of suspension parameters and flow field, at least with the two models investigated in this work. This is a situation that continues to hinder the fibre orientation problem. 6.4 2D Fibre Orientation The two dimensional fibre orientation distribution was computed over the entire plane of the contraction based on inlet conditions observed experimentally. A two way coupling was assumed between the fibre and fluid phases. It was shown that fi- bres preferentially align in the streamline direction and are essentially aligned along the centerline at the contraction exit. The model predictions compared well with experimental observations of fibre orientation over the contraction plane except near 71 Chapter 6. Summary and Conclusions the solid walls where the model predicted an increase in orientation isotropy and a shift in the mean fibre orientation angle in a direction away from the wall. Ex- perimental observations however showed no such behavior but instead showed that fibres align with the wall indefinitely. It is speculated that the discrepancy between the model predictions and experimental observations near the wall are based on the length of the fibres used in the experiments. More specifically, it is believed that the fibre is not affected by the larger flow gradients near the wall due to the fact that fibre length is comparable to the length scale of the wall boundary layer. It is concluded that a more detailed study of wall-fibre interactions is needed to determine sound relationships between fibre length, aspect ratio and possibly concentration and near wall fibre behavior. This is another situation that has in the past hindered, and will continue to hinder the fibre orientation problem. 6.5 Two Way Coupling The effect of the two way coupling on the orientation distribution and on the flow field has been investigated numerically. This was accomplished through an iterative procedure where the flow field was solved without the fibre phase present and the orientation equations solved based on the initial flow field. The fibre stress term was computed and used to recompute the flow field, this time with the fibre phase present, after which the orientation equations were solved again with the updated flow field. The process was repeated until convergence was reached. The results showed that the effect of the two way coupling on the 2D planar orientation distribution was 72 Chapter 6. Summary and Conclusions negligible. However the structure of the flow field was altered considerably. It was shown that the flow field changed from a parabolic-type profile, typical of Newtonian channel flows, to something which resembled a plug-type flow. The change in flow structure was attributed to a redirection of mass and momentum away from the walls towards the middle of the contraction, thus slowing the flow away from the walls and creating a more uniform flow distribution throughout the middle section of the channel. 73 Chapter 7 Future Work Several interesting findings have been made in this work. These include the obser- vation that fibres flocculate at concentrations above some critical value, the short- comings of the available models for the rotary diffusion coefficient, the near wall behavior of fibres, the boundary conditions on Ψ and the two-way interaction be- tween the suspension orientation state and the resulting flow field. It is recommended that these issues all be further studied in greater detail. 7.1 Mechanical Fibre-Fibre Interactions In order to study the effects of flocculation and mechanical entanglement on the orientation distribution function, a useful future study would investigate both the orientation state and concentration state of the suspension simultaneously and any relationships between the two. Furthermore, a stochastic model for the effects of mechanical fibre-fibre interactions would be particularly useful as this would allow researchers to model much higher concentration suspensions with an Eulerian model much like the one used in this work. Including the effect of variable fibre concentra- 74 Chapter 7. Future Work tions throughout the channel would also likely improve the predictions made in this work. 7.2 Rotary Diffusion Coefficient As was shown in this work, neither of the two models for the rotary diffusion coeffi- cient were able to model different suspension concentrations with one single closure parameter. That is CI in Equation 2.11, and λ in Equation 4.2 had to be determined whenever the suspension concentration changed. An interesting and valuable study into the correct scaling for the rotary diffusion coefficient could possibly eliminate this problem. Phan-Thien et al. (2002) addressed this issue in their work where they developed numerical correlations between the Jefferey orbit constant and the concentration parameter c · r as opposed to nL3, where c was the volume fraction of fibres and r the fibre aspect ratio. This may provide a starting point for such a study. 7.3 Near Wall Behavior of Fibres The behavior of fibres near solid boundaries still remains a big question mark in the fibre orientation problem. The few studies that have been carried out in this area, although extremely interesting, have not provided many solid answers. Future work should rigorously investigate both the behavior individual fibres and the suspension as a whole near solid boundaries. That is, boundary conditions would be determined 75 Chapter 7. Future Work for the orientation distribution function. This study could look into the effects of fibre length, aspect ratio, suspension concentration and possibly flow conditions on the near wall behavior of fibres and fibre suspensions. 7.4 Two way coupling The final recommendation for future work would be to measure the flow field associ- ated with the flow of non-dilute suspensions in contraction type flows and compare this with the numerical predictions made in this work. This would further prove the usefulness of the model for the two way coupling presented in this work and would allow researchers to work with more complex flow, such as turbulent flows for example, with improved confidence. 76 Bibliography Advani, S.G., Tucker, C.L., 1987 The Use of Tensors to Describe and Predict Fibre Orientation in Short fibre Composites J. Rheology. 31, 751–784. Advani, S.G., Tucker, C.L., 1990 Closure Approximations for Three- Dimensional Structure Tensors J. Rheology. 33, 367–387. Altan, M.C., Advani, S.G., Güceri, S.I., Pipes, R.B. 1989 On the De- scription of the Orientation State for Fibre Suspensions in Homogeneous Flows J. Rheology. 33, 1129–1155. Anczurowski, E., Mason, S.G. 1968 Particle Motions in Sheared Suspensions: Rotation of Rigid Spheroids and Cylinders Trans. Soc. Rheol. 12, 209–215. Asplund, G., Norman, B. 2004 Fibre Orientation Profile Over the Thickness of a Headbox Jet J. Pulp Paper Sci. 30(8), 217–221. Batchelor, G.K. 1970 Slender-body Theory for Particles of Arbitrary Cross- Section in Stokes Flow J. Fluid Mech. 44, 419–440 Bibbo, M.A., Dinh, S.M. & Armstrong, R.C. 1985 Shear Flow Froperties of Semi-Concentrated Fibre Suspension J. Rheol. 29(6) 905–929 77 Bibliography Bretherton, F.P. 1962 The Motion of Rigid Particles in a Shear Flow at Low Reynolds Number J. Fluid Mech. 14, 284 Chung, D.H., & Kwon, T.H. 2001 Improved Model of Orthotropic Closure Approximation for Flow Induced Fibre Orientation Polymer Composites 22, 636– 649 Cintra, J.S. & Tucker, C.L. 1995 Orthotropic Closure Approximations for Flow-Induced Fibre Orientation J. Rheol. 39, 1095–1122 Cox, R.G. 1970 The Motion of Long Slender Bodies in a Viscous Fluid J.Fluid Mech. 44, 791–810. Dinh, S.M. & Armstrong, R.C. 1984 A Rheological Equation of State for Semi- Concentrated Fibre Suspensions J. Rheol. 28(3), 207–227 Doi, M. & Edwards, S.F. 1984 The Theory of Polymer Dynamics Oxford Uni- versity Press New York Duffy, G. G. 2003 The Significance of Mechanistic-Based Models in Fibre Sus- pension Flow Nordic J. of Pulp Paper Sci. 18(1), 74–80. Folgar, F. 1983 Orientation behavior of fibres in concentrated suspensions Ph.D. thesis, University of Illinois at Urbana-Champaign, 1983 Folgar, F. & Tucker, C.L. 1984 Fibre Orientation Distribution in Concen- trated Suspensions: A Predictive Model J. Reinf. Plastics and Comp. 3 98–119 78 Bibliography Hämäläinen, T. & Hämäläinen, J. 2007 Modelling of Fibre Orientation in the Headbox Jet J. Pulp Paper Sci. 33(1), 49–53. Heath H.S., Olson, J.A., Buckley, K.R., Lapi, S., Ruth, T.J., & Mar- tinez, D,M. 2007 Visualization of the Flow of a Fibre Suspension Through a Sudden Expansion Using PET AIChE Journal 53(2), 327–334. Holm, R. & Söderberg, D. 2007 Modelling the Effect of Shear Flow on Fibre Orientation Anisotropy in a Planar Contraction Overview of forming literature, 1990-2000 Trans. 12th Fundamental Research Symposium, Oxford 1 431–558 Hyensjö, M., Dahlkild, A., Krochak, P., Olson, J. 2007 Shear Influence on Fibre Orientation Nordic Pulp and Paper Research Journal 22, 376–382. Martinez, D.M., Buckley, K., Jivan, S., Lindström, A., Thurigaswamy, R., Ruth, T.J. & Kerekes, R.J. 2001 Characterizing the Mobility of Fibres During Sedimentation Pulp Pap. Fund. Res. Symp. (Oxford) 1 225–254 Norman, B. & Söderberg, D. 2001 Overview of Forming Literature, 1990-2000 Trans. 12th Fundamental Research Symposium, Oxford 1 431–558 Jackson, W.C., Advani, S.G. & Tucker, C.L. 1985 Predicting the Orienta- tion of Short Fibres in Thin Compression Moldings J. Comp. Mater. 20 539–557 Jefferey, G.B. 1922 The Motion of Ellipsoidal Particles Immersed in a Visous Fluid Proc. Royal Soc. London Series A 102, 161–179 Khayat, R. E., & Cox, R. G. 1989 Inertia Effects on the Motion of Long Slender Bodies J.Fluid Mech. 209 435–462 79 Bibliography Koch, D.L. 1995 A Model for Orientational Diffusion in fibre Suspensions Phys. Fluids 7(8) 2086–2088 Leal, L.G. & Hinch, E.J. 1971 The Effect of Weak Brownian Rotations on Particles in Shear Flow Chem. Engng. Comm. 108, 381–401. Lin, J. & Zhang, L. 2002 Numerical Simulation of Orientation Distribution Func- tion of Cylindrical Particle Suspensions App. Math. Mech. 23, 906–912. Lin, J.Z., Gao, Z.Y., Zhou, K. & Chan, T.L. 2005 Mathematical Modeling of Turbulent Fibre Suspension and Successive Iteration Solution in the Channel Flow App. Math. Mod. 30, 1010–1020. Lipscomb, G.G. & Denn, M.M. 1988 Flow of Fibre Suspensions in Complex Geometries Journal of Non-Newtonian Fluid Mechanics 26 297–325 Olson, J.A., Frigaard, I., Chan, C., & Hämäläinen, J. P. 2004 Modeling a Turbulent Fibre Suspension Flowing in a Planar Contraction: The One Dimensional Headbox Int. J. Multiphase Flow 30 51–66 Parsheh, M., Brown, M.L., Aidun, C.K. 2005 On the Orientation of Stiff Fibres Suspended in a Turbulent Flow in a Planar Contraction J. Fluid Mech. 545, 245–269. Paschkewitz, J.S., Yves, D., Dimitropoulos, C.D., Shaqfeh, E.S., Parviz, M. 2004 Numerical Simulation of Turbulent Drag Reduction Using Rigid Fibres J. Fluid Mech. 518, 281–317. 80 Bibliography Phan-Thien, N., Fan, X.J., Tanner, R.I., Zheng, R. 2002 FolgarTucker constant for a Fibre Suspension in a Newtonian Fluid J.Non-Newt. Fluid Mech. 103, 251–260. Rahnama, M., Koch, D.L. & Shaqfeh, E.S.G. 1995 The Effect of Hydrody- namic Interactions on the Orientation Distribution in a Fibre Suspension Subject to Simple Shear Flow Phys Fluids 7(3) 487–506 Rahnama, M., Koch, D.L. & Cohen, C. 1995 Observations of Fibre Ori- entation in Suspensions Subject to Planar Extensional Flows Phys. Fluids 7(8) 1811–1817 Ranganathan, S. & Advani, S.G. 1991 Fibre-Fibre Interactions in Homoge- neous Flows of Nondilute Suspensions Journal of Rheology 35(8) 1499–1522 Shaqfeh, E.S.G. & Fredrickson, G.H. 1990 The Hydrodynamic Stress in a Suspension of Rods Phys Fluids A 2(1) 7–24 Stover, C.A., & Cohen, C. 1990 The Motion of Rodlike Particles in the Pressure-Driven Flow Between Two Flat Plates Rheologica Acta 29 192–203 Stover, C.A., Koch D.L.& Cohen, C. 1992 Observation of Fibre Orientations in Simple Shear Flow of Semi Dilute Suspensions J. Fluid Mech. 238 277–296 Ullmar, M. & Norman, B. 1997 Observation of Fibre Orientation in a Head- box Nozzle at Low Consistency TAPPI Proceedings, Engineering and Papermaking conference pg 685 81 Bibliography Ullmar, M. 1998 On Fibre Alignment Mechanisms in a Headbox Nozzle Licentiate Thesis, Royal Institute of Technology Stockholm, Sweden VerWeyst, B.E. & Tucker, C.L. 1999 Fibre Orientation in 3-D Injection Molded Features: Prediction and Experiment International Polymer Processing 4 409–420 VerWeyst, B.E. & Tucker, C.L. 2002 Fibre Suspensions in Complex Geome- tries: Flow/fibre Coupling Canadian Journal of Chemical Engineering 80 1093– 1106 Yasuda, K., Kyuto, T. & Mori, N. 2004 An Experimental Study of Flow- Induced Fibre Orientation and Concentration Distributions in a Concentrated Sus- pension Flow Through a Slit Channel Containing a Cylinder Rheologica Acta (2004) 43 137-145 Zhang, X. 1998 Fibre Orientation in a Headbox Masters Thesis, The University of British Columbia Vancouver, Canada 82 Appendix A Experimental System A schematic diagram of the experimental system is shown below, see Figure A.1. Figure A.1: A schematic diagram of the experimental system. The experimental flow loop is constructed entirely of plexiglass with fully smooth walls, 5mm thick. Upstream of the experimental section is a short rectangular channel section where the flow is allowed to become fully developed. A dimensioned 83 Appendix A. Experimental System sketch of the experimental section is shown below, see figure A.2. Figure A.2: A dimensioned sketch of the experimental section. The suspension is continously recirculated with a Masterflex B/T variable speed peristaltic pump (www.coleparmer.com). The specifications for this pump are given below. • Tubing Size: 91 • Max Flow Rate: 37 Litres/min • Max Speed: 350 r/min • Input Voltage: 90 - 130 Vrms • Motor horsepower: 0.5 hp 84 Appendix A. Experimental System The index of refraction matched suspension consists of a boroscilicate glass fibre phase suspended in a glycerine fluid phase. The physical data for the suspending fluid, glycerine, at 20 0C is given below: • Density: 1260 kg/m3 • Viscosity: 1.49 kg/ms • Index of Refraction: 1.474 The physical data for the boroscilicate fibres used in the experiments is given be- low: • Density: 2230 kg/m3 • Index of Refraction: 1.472 Image sequences are captured using a progressive scan Basler A201b monochrome CCD camera mounted with an F-mount Micro-NIKKOR 105 mm lens. The camera lens is positioned 130 cm from the experimental section. This particular camera distance was chosen so that: (a) the entire contraction could be captured in a single image, (b) the diameter of an imaged fibre was exactly two pixels, and (c) the relative change in optical path due to the differences of index of refraction as light travels from the suspension to the camera lens is negligible. The specifications for the CCD camera and imaging system are given below. • 10 bit grey scale camera with 1008× 1016 pixel spatial resolution • Max framing rate of 5 frames per second 85 Appendix A. Experimental System • F-mount Micro-NIKKOR 105 mm lens • Schott-Fostec fibre optic dual backlight (tungsten halogen) 86 Appendix B Matlab Program Code In this section the Matlab codes used to solve various aspects of this work are pre- sented. This includes the programs used to track fibres within the experimental images of the flowing suspensions and the programs used to analyze the raw data generated by the fibre tracking algorithms. Also included in this section is the code used to perform the Eularian numerical predictions of fibre orientation. That is, the numerical solvers used to solve Equations 4.4 and 5.1 are given B.1 Fibre Tracking and Analysis B.1.1 Fibre Tracking The programs used to track fibres are called trackFibres.m, get2DThreshImage 3.m, idFibresE cor.m, and FindEndPoints.m. trackFibres.m is the main driver program. It calls the functions get2DThreshImage 3.m, a function used to accept an original image and return a black and white image consisting of white fibres on a black background. This image is then used by the function idFibresE cor.m, a function 87 Appendix B. Matlab Program Code that returns the values of the centre of area, and two end points of each observed fibre in the black and white image. This function also filters images for unwanted debris and crossing fibres. The values of the fibre end points is determined by the function FindEndPoints.m called within idFibresE cor.m. These values are then saved in a seperate file. The process is repeated for every image in a sequence. trackFibers.m %%% Filename: trackFibers.m %%% Written by: Paul J. Krochak %%% Last Modified: 22-Feb-2008 clear all home % Search lenghts dx1 = 5; dy1 = 5; dx2 = 5; dy2 = 5; dt = 200e-3; seqNr = 1; load roiMsk_jan_10_07 %This is a manually defined mask to isolate the flow region for seqNr = 1:10, %Define sequence file location movPath = [’F:\CNum = 0\Jan 17_07\’]; movName = [’Sequence’,num2str(seqNr),’.avi’]; readName = [movPath,movName]; mov = aviread(readName); [m,nImages] = size(mov); nImSet = 2; %Number of images through which a fiber will try to be tracked imOrg = mov(1).cdata; [roiMsk1,xi1,yi1] = roipoly(imOrg); [roiMsk2,xi2,yi2] = roipoly(imOrg); for vw = 1:1, % side view and top view handled seperately if vw == 1, % side view 88 Appendix B. Matlab Program Code roiMsk = roiMsk1; xi = xi1; yi = yi1; else %top view roiMsk = roiMsk2; xi = xi2; yi = yi2; end imNmbr = 0; setNmbr = 0; while imNmbr < nImages-1, setNmbr = setNmbr+1; imNmbr = imNmbr+1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % identify fibres from the first image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if vw == 1, imOrg = mov(imNmbr).cdata; imRoi = imOrg.*uint8(roiMsk); imRoi = imRoi(:,1:508); im1 = get2DThreshImage_3(imRoi); else imOrg = mov(imNmbr).cdata; imRoi = imOrg.*uint8(roiMsk); imRoi = imRoi(:,520:1008); im1 = get2DThreshImage_3(imRoi); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Save Eularian Statistical Data From First Image %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [fibSP(imNmbr),fibCOG(imNmbr),fibEP(imNmbr),NumOfFibres1] =... idFibresE_cor(im1); end %%% End while imNbr loop savename = [’C:\Documents and Settings\Administrator\My Documents\... 3D Fiber Tracking\Results\CN = 32\Apr_4_07\seq_’,num2str(seqNr),... ’view_’,num2str(vw),’_Ori’]; save(savename,’fibSP’,’fibCOG’,’fibEP’,’xi’,’yi’) clear fibSP fibCOG fibEP xi yi 89 Appendix B. Matlab Program Code end %End vw loop....vw = view = [top,bottom] end % End seqNr loop get2DThreshImage 3.m % Filename: get2DThreshImage_3.m % Created by: Paul Krochak % Last modified: 20-Feb-2008 function M = get2DThreshImage_3(imRoi) %imRoi is the original image consisting of flow region only s = size(imRoi); % Perform thresholding and edge-detection BWs = edge(imRoi, ’prewitt’, (graythresh(imRoi) * 0.15)); % Define structure elements for dilation of fibre edges se90 = strel(’line’, 3, 90); se0 = strel(’line’, 2, 0); se45 = strel(’line’, 2, 45); se60 = strel(’line’, 2, 60); BWsdil = imdilate(BWs, [se90 se0 se45 se60]); % Fill fibres with white pixels BWdfill = imfill(BWsdil, ’holes’); BWnobord = imclearborder(BWdfill, 8); % Erode fibres back to original size seD = strel(’diamond’,1); BWfinal = imerode(BWnobord,seD); BWfinal = imerode(BWnobord,seD); M = BWfinal; [A,NumOfFibres] = bwlabel(M,4); eps = 10; eps2 = 5; %%%% Filter Debris from image 90 Appendix B. Matlab Program Code if (NumOfFibres), for j = 1 : NumOfFibres; [x,y] = find(A == j); if ((abs(max(x)-min(x))<= eps) & (abs(max(y)-min(y))<=eps))... |(abs(max(x)-min(x))-(max(y)-min(y)) <=eps2), M(x,y) = 0; else end end end get2DThreshImage 3.m function [sp,cog,ep,NumOfFibres] = idFibresE_cor(im) % labels fibres from tresholded image and returns coordimates, number of % fibres found and end points % Written by: Paul Krochak % Last modified: 20-Feb-2008 fibNr = 0; [A,NumFibres] = bwlabel(im,8); NumFibres if (NumFibres), for j = 1 : NumFibres, [x,y] = find(A == j); %%%%%%%Check for x-ing fibers: If cogVal is true (1) then COG is %% on the fibre and fibre is assumed not to be crossed with another %% If detected fibre length is greater than 55 pixels it is likely %% two fibres crossing [sp0,cog0,ep0] = FindEndPoints(x,y); L0 = sqrt((ep0(1)-sp0(1))^2+(ep0(2)-sp0(2))^2); cogVal = A(cog0(1),cog0(2)); if cogVal & (L0 <= 55) fibNr = fibNr+1; %fibre(j).coord = [x,y]; [sp(1).coord(fibNr,:),cog(1).coord(fibNr,:),ep(1).coord(fibNr,:)]... = FindEndPoints(x,y); 91 Appendix B. Matlab Program Code end end else sp(1).coord = []; cog(1).coord = []; ep(1).coord = []; end NumOfFibres = fibNr; FindEndPoints.m function [sp,cog,ep] = FindEndPoints(x,y) % Finds 2d-end points and center of gravity of fibres % Written by: Paul Krochak % Last modified: 20-Feb-2008 [minX,minIndX] = min(x); [minY,minIndY] = min(y); [maxX,maxIndX] = max(x); [maxY,maxIndY] = max(y); %Determine orientation of fibre and compute end points if x(maxIndY)-x(minIndY) < 0 sp = [x(maxIndX),y(minIndY)]; cog = [x(minIndX)+round((x(maxIndX)-x(minIndX))/2),... y(maxIndY)+round((y(minIndY)-y(maxIndY))/2)]; ep = [x(minIndX),y(maxIndY)]; else sp = [x(minIndX),y(minIndY)]; cog = [x(minIndX)+round((x(maxIndX)-x(minIndX))/2),... y(minIndY)+round((y(maxIndY)-y(minIndY))/2)]; ep = [x(maxIndX),y(maxIndY)]; end 92 Appendix B. Matlab Program Code B.1.2 Analysis of Fibre Observations The data obtained from the fibre tracking algorithm is converted to position and orientation data with the code presented below. Included in this is the routine used to measure the fibre orientation distribution along the central streamline only, however a similar program is used to measure the orientation distribution across the entire xy-plane. The routine called analyzeOri1D Str.m is the main driver program for doing so. This program goes through the fibre data from each image, that is the data representing the centre of area and two end points of each observed fibre as measured with the code from §B.1.1, and computes the orientation of each fibre. These orientations are binned into the correct spatial cell based on the centre of area of the fibre through the helper function binSpace.m. The orientation distribution function is then computed for each spatial cell. It should be pointed out that we only present the analysis of one data set. That is, the code shown below is for the 10 sequences obtained on Jan 10, 2007, however the exact same procedure is used for all other data sets obtained in this work. analyzeOri1D Str.m % File: analyzeOri1D_Str.m % Written by: Paul J. Krochak % Last modified: Dec 20 2006 clear all home xN = round(1016/(7.8154*5)); %Number of spatial points in x-direction x = linspace(round(1016/26),1016,26); 93 Appendix B. Matlab Program Code DX = abs(x(2)-x(1)); %%%% Initialize spactial distribution to empty struct for k=1:xN spCount(k).phi = []; end % Structs used to store position and orientation data oriCL.loc = []; oriCL.phi = []; eps = DX/2; %%%%%%Window Size fibInd = 0; vw = 1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% ANALYZE JANUARY 9 DATA %%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% loadPath = [’C:\Documents and Settings\Paul Krochak\My Documents\... Fiber Tracking\3D Fiber Tracking\Results\CNUM = 8\Jan_9_07\’]; loadName = [’C:\Documents and Settings\Paul Krochak\My Documents\... Fiber Tracking\3D Fiber Tracking\Results\CNUM = 8\Jan_9_07\seq_1view_1_Ori.mat’]; load(loadName); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% Determine CenterLine and centrLn direction (phiCL) %%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% y0 = round(min(xi)) + round(0.5*(max(xi)-min(xi))); x0 = 0.5*(yi(1)+yi(2)); ye = round(xi(4)) + round(0.5*(xi(3)-xi(4))); xe = 0.5*(yi(3)+yi(4)); phiCL = atan((ye-y0)/(xe-x0)); y1 = xi(1);x1 = yi(1); y2 = xi(2);x2 = yi(2); y3 = xi(3);x3 = yi(3); y4 = xi(4);x4 = yi(4); yCL = y0 + (y0-ye)/(x0-xe)*x; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% Define Analysis Region %%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% y0_t = y0 + 0.3*(y2-y0); ye_t = ye + 0.3*(y3-ye); 94 Appendix B. Matlab Program Code y0_b = y0 + 0.3*(y1-y0); ye_b = ye + 0.3*(y4-ye); kt = ye_t + (y0_t-ye_t)/(x2-x3)*x; %%% Top perimeter of Window kb = ye_b + (y0_b-ye_b)/(x1-x4)*x; %%% Bot perimeter of Window for seqNr = 1:10, %10 sets of sequences obtained for this experimental set loadSeq = [’seq_’,num2str(seqNr),’view_’,num2str(vw),’_Ori’,’.mat’]; loadName = [loadPath,loadSeq]; load(loadName); [n,nIms] = size(fibCOG); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Go through each Image %%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for n=1:nIms, clear cog_n cog_n = fibCOG(n).coord; for i=1:26 xc = x(i); % Find fibres in region of interest in i’th spatial bin expr = (abs(cog_n(:,1)-xc) <= eps) & ((kt(i) >= cog_n(:,2))... & (cog_n(:,2) >=kb(i))); Ind = find(expr); if Ind for j=1:length(Ind) fibInd = fibInd+1; dx = fibSP(n).coord(Ind(j),1) - fibEP(n).coord(Ind(j),1); dy = fibSP(n).coord(Ind(j),2) - fibEP(n).coord(Ind(j),2); % Correct for 3 pixel shift in end point definition if dy < 0 dy = dy+3; elseif dy > 0 dy = dy-3; end %save orientation and position data in struct oriCL(fibInd).phi = atan(dy/dx)+phiCL; oriCL(fibInd).loc = fibCOG(n).coord(Ind(j),:); end else 95 Appendix B. Matlab Program Code end end end end spCount = binSpace(oriCL,spCount,yi); binSpace.m function [distr] = binSpace(ori,oldDistr,yi) % Filename: binSpace.m % Written by: Paul J Krochak % Last modified: 20-Feb-2008 %%% Resolution of 7.8154 pixels per 5mm x_L = 0:5*7.8154:1016; %discretized centerline locations (pixel dimensionality). if ~isempty(ori(1).loc) for j=1:length(ori), x1(j) = ori(j).loc(:,1); end x_f = abs(x1 - max(yi)); distr = oldDistr; eps = 3.9077*5; %Halfway between 2 spactial grid points for i=1:length(x_L) clear Ind Ind = find(abs(x_f - x_L(i)) <= eps); if Ind, distr(i).phi = [distr(i).phi,ori(Ind).phi]; end end % Else is envoked if no fibres were found else distr = []; end 96 Appendix B. Matlab Program Code B.2 Numerical Prediction of Fibre Orientation In this section the Matlab code used to solve the Fokker-Plank equations describing the evolution of the orientation distribution function, Ψ, are presented. Included in this section is the code used to solve Equation 4.4, that is, the code used to predict the 3D fibre orientation distribution along the central streamline. Also shown is the code used to solve Equation 5.1 to predict the 2D planar orientation distribution over the xy-plane. This is followed by the code used to compute the additional fibre stress term described by Equation 2.15. B.2.1 Fokker-Plank Equation Solvers The two routines used to solve Equations 4.4 and 5.1 are solvePsi3CL.m and solveP- siStr.m respecively. The code shown here for solving Equation 4.4 is for the case when Equation 4.2 is used to define the diffusion coefficient. The helper function, geta4 CL.m, used in this code computes and returns the fourth order orientation tensor. solvePsi3CL.m % Filename: solvePsi3CL.m % Written by: Paul J Krochak % Last modified: 20-Feb-2008 clear all home 97 Appendix B. Matlab Program Code CN = 8; %CN = nL^3 LAM = 0.005:0.0002:0.008; r = 50; L = 0.005; R = 10; Nphi = 200; phi = linspace(-pi/2,pi/2,Nphi); phi = linspace(-pi/2,pi/2,Nphi); dphi = phi(2)-phi(1); Ntheta = 200; theta = linspace(0,pi,Ntheta); dth = theta(2)-theta(1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% Define Velocity Field %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nnod = 50; u0 = 1.0; Len = 1.0; x = linspace(0,1,Nnod); x = x’; u = u0./(1-(1-1/R)*(x/Len)); dudx = u0*(1-1/R)./((1-(1-1/R)*(x/Len)).^2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%% Load Experimental Data for phi projection %%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Psi_exp_p = dlmread(’PsiP_Exp_CL_Crowd_8.dat’); [m n] = size(Psi_exp_p); NnodE = m; NphiE = n; phiE = linspace(-pi/2,pi/2,NphiE); IND = find(isnan(Psi_exp_p)); Psi_exp_p(IND) = 0; %%%%%%%%% Interpolate Psi_exp to refine mesh %%%%%%%%%%%%%%%%%% Psi_exp_p = interpPsi(phiE,Psi_exp_p,phi); initPsi_p = Psi_exp_p(4,:); clear Psi_exp_p phiE NphiE NnodE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 98 Appendix B. Matlab Program Code %%%%%%%%% Load Experimental Data for gamma projection %%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Psi_exp_g = dlmread(’PsiG_Exp_CL_Crowd_8.dat’); [m n] = size(Psi_exp_g); NnodE = m; NthE = n; thE = linspace(0,pi,NthE); IND = find(isnan(Psi_exp_g)); Psi_exp_g(IND) = 0; %%%%%%%%% Interpolate Psi_exp to refine mesh %%%%%%%%%%%%%%%%%% Psi_exp_g = interpPsi(thE,Psi_exp_g,theta); initPsi_g = Psi_exp_g(4,:); clear Psi_exp_g phiE NthE initPsi = initPsi_p’*initPsi_g; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Re-nomralize inlet distribution %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for kk=1:Ntheta, t = theta(kk); sumPsi(kk) = sin(t)*(initPsi(1,kk) + 2*sum(initPsi(2:(Nphi-1),kk))... +initPsi(end,kk)); end intPsi = dphi*dth/4*(sumPsi(1) + 2*sum(sumPsi(2:(Ntheta-1)))+sumPsi(end)); initPsi(:,:) = abs(1/intPsi*initPsi(:,:)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RES = []; Psi = 1/(2*pi)*ones(Nnod,Nphi,Ntheta); Psi(1,:,:) = initPsi; Psi0 = norm(Psi(:),2); load Psi3DExp_CN_8_lambda_0.0065.mat Psi Psi_old = Psi; for iL =1:length(LAM), lambda = LAM(iL) tol = 1e-4; 99 Appendix B. Matlab Program Code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Define Diffusion Coefficient %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% a4 = geta4_CL(Psi,phi,theta,Nnod,Nphi,Ntheta); Dr = CN/(log(r)^2)*lambda*a4.*dudx; Dr_old = Dr; Dr_err = 1e6; while Dr_err > tol, Psi_old = Psi; tol = 1e-4; itemax = 100000; ite = 0; err = 1e6; while (ite < itemax) & (err > tol), ite = ite + 1; %Enforce BC at inlet Psi(1,:,:) = initPsi; % Psi(1,:,:) = 1/(2*pi); for j=1:Nnod-1, dx = x(j+1)-x(j); %Enforce BC’s in theta for k=2:Ntheta-1 %Enforce BC’s in phi Psi(j+1,Nphi-1,k) = Psi(j+1,Nphi,k); for i=2:Nphi-1, p = phi(i); t = theta(k); phidot = -sin(p)*cos(p)*dudx(j); thetadot = dudx(j)*cos(p)^2*sin(t)*cos(t); dphidot = -cos(2*p)*dudx(j); dthetadot = dudx(j)*cos(p)^2*cos(2*t); cVel = u(j)/dx; Psi(j+1,i,k) = (cVel*Psi(j,i,k) + Dr(j)/(sin(t)^2*dphi^2)*... (Psi(j+1,i-1,k)+Psi(j+1,i+1,k))-0.5*(phidot*1/dphi*... (Psi(j+1,i+1,k)-Psi(j+1,i-1,k)))... +Dr(j)/(dth^2)*(Psi(j+1,i,k+1)+Psi(j+1,i,k-1))... +0.5*Dr(j)/dth*cot(t)*(Psi(j+1,i,k+1)-Psi(j+1,i,k-1))... 100 Appendix B. Matlab Program Code -0.5/dth*thetadot*(Psi(j+1,i,k+1)-Psi(j+1,i,k-1)))... *1/(cVel + 2*Dr(j)*1/(sin(t)^2*dphi^2)+2*Dr(j)*... 1/(dth^2) + dphidot + cot(t)*thetadot + dthetadot); end end end Psi = normPsi(phi,theta,Psi,Nnod,Nphi,Ntheta); delta = Psi-Psi_old; err = norm(abs(delta(:)),2)/Psi0; Psi_old = Psi; end a4 = geta4_CL(Psi,phi,theta,Nnod,Nphi,Ntheta); Dr = CN/(log(r)^2)*lambda*a4.*dudx; Dr_err = norm(abs(Dr_old-Dr),2) Dr_old = Dr; end %End Dr While loop saveStr = [’save Psi3DExp_CN_8_lambda_’,num2str(lambda),’.mat Psi’]; eval(saveStr); end geta4 CL.m function a4 = geta4_CL(Psi,phi,theta,Nnod,Nphi,Ntheta) % Filename: geta4_CL.m % Written by: Paul J Krochak % Last modified: 20-Feb-2008 %%%% Only one component ever gets used along CL a4 = zeros(Nnod,1); dp = abs(phi(2)-phi(1)); dt = abs(theta(2)-theta(1)); A = zeros(Nphi,Ntheta); for j=1:Nnod, for ii=1:Nphi, for kk=1:Ntheta, 101 Appendix B. Matlab Program Code A(ii,kk) = Psi(j,ii,kk); end end for i=2:Nphi, for k=1:Ntheta, t = theta(k); sumPsi(k) = sin(t)*(A(1,k)*(cos(1)*sin(t))^4+... 2*sum(A(2:(Nphi-1),k).*(cos(phi(2:(Nphi-1)))’*sin(t)).^4)+... A(end,k)*(cos(phi(end))*sin(t))^4); end a4(j,1) = dp*dt/4*(sumPsi(1) + 2*sum(sumPsi(2:(Ntheta-1)))+sumPsi(end)); clear sumPsi end end solvePsiStr.m % Filename: solvePsiStr.m % Written by: Paul J Krochak % Last modified: 20-Feb-2008 clear all home CN = 8; %CN = nL^3 CI = 0.0037; r = 50; L = 0.005; Nphi = 500; phi = linspace(-pi/2,pi/2,Nphi); dphi = phi(2)-phi(1); nStrLns = 13; %Number of streamlines for nSt = 1:nStrLns, %load flow field data for nSt’th streamline strName = [’TecStrDat_’,num2str(nSt),’.dat’] fieldDat = dlmread(strName); x = fieldDat(:,1); 102 Appendix B. Matlab Program Code y = fieldDat(:,2); u = fieldDat(:,3); v = fieldDat(:,4); dudx = fieldDat(:,5); dvdx = fieldDat(:,6); dudy = fieldDat(:,7); dvdy = fieldDat(:,8); modV = sqrt(diff(x).^2+diff(y).^2); Nnod = length(u); IF = find(isnan(fieldDat)) clear fieldDat %%%% Compute Streamline Parameter s (Arc Length) s(1) = 0; for i=2:Nnod, s(i) = sum(modV(1:i-1)); end dsdx = diff(s)./diff(x’); dsdy = diff(s)./diff(y’); Psi = 1/pi*ones(Nnod,Nphi); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Load Experimental Data and set Inlet Condition to Experiment %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if (nSt==1) | (nSt==2) | (nSt==3), Psi_exp = dlmread(’Psi_Exp_SL_2_Crowd_8.dat’); [m n] = size(Psi_exp); NnodE = m; NphiE = n; phiE = linspace(-pi/2,pi/2,NphiE); IND = find(isnan(Psi_exp)); Psi_exp(IND) = 0; %%%%%%%%% Interpolate Psi_exp to refine mesh %%%%%%%%%%%%%%%%%% Psi_exp = interpPsi(phiE,Psi_exp,phi); initPsi = Psi_exp(4,:); clear Psi_exp phiE elseif (nSt==4) | (nSt==5), Psi_exp = dlmread(’Psi_Exp_SL_1_Crowd_8.dat’); [m n] = size(Psi_exp); NnodE = m; 103 Appendix B. Matlab Program Code NphiE = n; phiE = linspace(-pi/2,pi/2,NphiE); IND = find(isnan(Psi_exp)); Psi_exp(IND) = 0; %%%%%%%%% Interpolate Psi_exp to refine mesh %%%%%%%%%%%%%%%%%% Psi_exp = interpPsi(phiE,Psi_exp,phi); initPsi = Psi_exp(4,:); clear Psi_exp phiE elseif (nSt==6) | (nSt==7) | (nSt==8), Psi_exp = dlmread(’Psi_Exp_CL_Crowd_8.dat’); [m n] = size(Psi_exp); NnodE = m; NphiE = n; phiE = linspace(-pi/2,pi/2,NphiE); IND = find(isnan(Psi_exp)); Psi_exp(IND) = 0; %%%%%%%%% Interpolate Psi_exp to refine mesh %%%%%%%%%%%%%%%%%% Psi_exp = interpPsi(phiE,Psi_exp,phi); initPsi = Psi_exp(4,:); clear Psi_exp phiE elseif (nSt==9) | (nSt==10), Psi_exp = dlmread(’Psi_Exp_SL_p1_Crowd_8.dat’); [m n] = size(Psi_exp); NnodE = m; NphiE = n; phiE = linspace(-pi/2,pi/2,NphiE); IND = find(isnan(Psi_exp)); Psi_exp(IND) = 0; %%%%%%%%% Interpolate Psi_exp to refine mesh %%%%%%%%%%%%%%%%%% Psi_exp = interpPsi(phiE,Psi_exp,phi); initPsi = Psi_exp(4,:); clear Psi_exp phiE elseif (nSt==11) | (nSt==12) | (nSt==13), Psi_exp = dlmread(’Psi_Exp_SL_p2_Crowd_8.dat’); [m n] = size(Psi_exp); NnodE = m; NphiE = n; phiE = linspace(-pi/2,pi/2,NphiE); 104 Appendix B. Matlab Program Code IND = find(isnan(Psi_exp)); Psi_exp(IND) = 0; %%%%%%%%% Interpolate Psi_exp to refine mesh %%%%%%%%%%%%%%%%%% Psi_exp = interpPsi(phiE,Psi_exp,phi); initPsi = Psi_exp(4,:); clear Psi_exp phiE end Psi_old = Psi; if nSt==7, v = zeros(size(v)); dsdy = zeros(size(dsdy)); dvdx = zeros(size(dvdx)); dvdy = zeros(size(dvdy)); dudy = zeros(size(dudy)); end Dr = CI*sqrt(dudx.^2+dudy.^2+dvdx.^2+dvdy.^2); Dr_err = 1e6; tol = 1e-6; itemax = 100000; ite = 0; err = 1e6; while (ite < itemax) & (err > tol), ite = ite + 1; %Enforce BC at inlet Psi(1,:) = initPsi; for j=1:Nnod-1, ds = s(j+1)-s(j); %Enforce BC’s in phi Psi(j+1,1) = Psi(j+1,Nphi); for i=2:Nphi-1, p = phi(i); phidot = -sin(p)*cos(p)*dudx(j)-sin(p)^2*dudy(j)+... cos(p)^2*dvdx(j)+sin(p)*cos(p)*dvdy(j); dphidot = -cos(2*p)*dudx(j)-sin(2*p)*dudy(j)-sin(2*p)*... dvdx(j)+cos(2*p)*dvdy(j); Psi(j+1,i) = ((u(j)*dsdx(j)+v(j)*dsdy(j))*1/ds*Psi(j,i)+... 105 Appendix B. Matlab Program Code Dr(j)/(dphi^2)*(Psi(j+1,i-1)+Psi(j+1,i+1))-... 0.5*phidot*1/dphi*(Psi(j+1,i+1)-Psi(j+1,i-1)))*... 1/((u(j)*dsdx(j)+v(j)*dsdy(j))*1/ds + ... 2*Dr(j)*1/(dphi^2) + dphidot); end end Psi = normPsi(phi,Psi,Nnod,Nphi); delta = Psi-Psi_old; err = norm(abs(delta(:)),2); Psi_old = Psi; end saveStr = [’save Psi_1D_Str_’,num2str(nSt),’.dat Psi’]; eval(saveStr); clear Psi Dr s x y u v dudx dudy dvdx dvdy dsdx dsdy end B.2.2 Computation of the Fibre Stress The Matlab code used to compute the fibre stress term and it’s gradient is given below. This routine first computes the relationship given by Equation 2.15, then it computes the gradient of this term which is treated as a momentum source term in the fluid momentum equations. The fibre stress is computed along each streamline for which a solution was obtained to Equation 5.1 then it’s gradient is interpolated onto the mesh used by Fluent to solve for the flow field. The resulting file containing this data is later loaded into the Fluent environment to solve the fluid flow equations. This routine calls the helper function get fa4.m which computes the following fa4 = E :< pppp > (B.1) 106 Appendix B. Matlab Program Code get phiStressStr.m % Filename: get_phiStressStr.m % Written by: Paul J Krochak % Last modified: 20-Feb-2008 % clear all home % Nphi = 500; phi = linspace(-pi/2,pi/2,Nphi); dphi = abs(phi(2)-phi(1)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% Compute Fiber Orientation Momentum Effect %%%%%% %%%%%%%%%%%%%%%% For Each StreamLine %%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Define Apparent Viscosity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mu = 0.00103; rho = 998.2; n = 6.4000e+007; %number density L = 0.005; %fiber length d = 1e-4; c = pi/4*n*L*d^2; %volume fraction r = 50; mu_f = mu*c*r^2/(log(1/c)+log(log(1/c))+1.439); nStrLns = 13; %number of streamlines M = []; X = []; Y = []; for n=1:nStrLns, %Load flow field data for n’th streamline strName = [’TecStrDat_’,num2str(n),’.dat’] fieldDat = dlmread(strName); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 107 Appendix B. Matlab Program Code x = fieldDat(:,1); y = fieldDat(:,2); u = fieldDat(:,3); v = fieldDat(:,4); dudx = fieldDat(:,5); dvdx = fieldDat(:,6); dudy = fieldDat(:,7); dvdy = fieldDat(:,8); modV = sqrt(diff(x).^2+diff(y).^2); Nnod = length(u); clear fieldDat X = [X;x]; Y = [Y;y]; s(1) = 0; for i=2:Nnod, s(i) = sum(modV(1:i-1)); end dsdx = diff(s)./diff(x’); dsdy = diff(s)./diff(y’); if n==7, v = zeros(size(v)); dsdy = zeros(size(dsdy)); dvdx = zeros(size(dvdx)); dvdy = zeros(size(dvdy)); dudy = zeros(size(dudy)); end % Load Psi for n’th streamline loadStr = [’Psi_1D_Str_’,num2str(n),’.dat’]; Psi = dlmread(loadStr); [Nnod,Nphi] = size(Psi); phi = linspace(-pi/2,pi/2,Nphi); dphi = abs(phi(2)-phi(1)); a4 = geta4(Psi,phi,Nnod,Nphi); f_a4 = getf_A4(a4,Nnod,dudx,dvdx,dudy,dvdy); f_x = diff(f_a4(:,1,1)’)./diff(s).*dsdx + diff(f_a4(:,1,2)’)./diff(s).*dsdy; f_y = diff(f_a4(:,2,1)’)./diff(s).*dsdx + diff(f_a4(:,2,2)’)./diff(s).*dsdy; M = [M;mu_f*[f_x(1),f_y(1)];mu_f*[f_x’,f_y’]]; clear Psi Dr s x y u v dudx dudy dvdx dvdy dsdx dsdy 108 Appendix B. Matlab Program Code end % Load Fluent Mesh load 2dExpRectMesh.mat Nnod = length(coords); Nelems = length(elems); fxI = griddata(X,Y,M(:,1),coords(:,2),coords(:,3)); fyI = griddata(X,Y,M(:,2),coords(:,2),coords(:,3)); IX = find(isnan(fxI)); fxI(IX) = 0; IY = find(isnan(fyI)); fyI(IY) = 0; fid = fopen(’C:\Fluent.Inc\fluent6.2.16\Sims\phiStress.dat’,’w’); % For each element average the fibre stress and save in a file for Fluent for n=1:Nelems, fx_0 = fxI(elems(n,1)); fx_1 = fxI(elems(n,2)); fx_2 = fxI(elems(n,3)); fx_3 = fxI(elems(n,4)); fy_0 = fyI(elems(n,1)); fy_1 = fyI(elems(n,2)); fy_2 = fyI(elems(n,3)); fy_3 = fyI(elems(n,4)); fx(n) = 1/4*(fx_0+fx_1+fx_2+fx_3); fy(n) = 1/4*(fy_0+fy_1+fy_2+fy_3); fprintf(fid,’\n %f, %f’,fx(n),fy(n)); end fclose(’all’); get fa4.m function f_a4 = getf_A4(a4,jmax,dudx,dvdx,dudy,dvdy) % Helper function used to return f_a4 = E:<pppp> in % the orientation tensor model equation 109 Appendix B. Matlab Program Code f_a4 = zeros(jmax,2,2); f_a4(:,1,1) = a4(:,1,1,1,1).*dudx+a4(:,1,1,2,1).*1/2.*... (dudy+dvdx)+a4(:,1,1,1,2).*1/2.*(dudy+dvdx)+a4(:,1,1,2,2).*dvdx; f_a4(:,1,2) = a4(:,1,2,1,1).*dudx+a4(:,1,2,2,1).*1/2.*... (dudy+dvdx)+a4(:,1,2,1,2).*1/2.*(dudy+dvdx)+a4(:,1,2,2,2).*dvdx; f_a4(:,2,1) = a4(:,2,1,1,1).*dudx+a4(:,2,1,2,1).*1/2.*... (dudy+dvdx)+a4(:,2,1,1,2).*1/2.*(dudy+dvdx)+a4(:,2,1,2,2).*dvdx; f_a4(:,2,2) = a4(:,2,2,1,1).*dudx+a4(:,2,2,2,1).*1/2.*... (dudy+dvdx)+a4(:,2,2,1,2).*1/2.*(dudy+dvdx)+a4(:,2,2,2,2).*dvdx; 110
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The orientation state of semi-dilute rigid fibre suspensions...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The orientation state of semi-dilute rigid fibre suspensions in a linearly contracting channel Krochak, Paul Joseph 2008
pdf
Page Metadata
Item Metadata
Title | The orientation state of semi-dilute rigid fibre suspensions in a linearly contracting channel |
Creator |
Krochak, Paul Joseph |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | This work investigates the effects of long range hydrodynamic fibre-fibre interactions on the orientation state of a semi-dilute, rigid fibre suspension flowing through a linear contracting channel under laminar flow conditions. The effects of fibre-fibre interactions are modeled mathematically, the governing equations solved numerically and the predicted results compared with experimental observations. The theoretical model is based on the assumption that the orientation state of the suspension can be completely described by a probability distribution function and that fibre-fibre interactions are random in nature, thus giving rise to a diffusion-type process. The orientation distribution evolves spatially according to a Fokker-Plank type equation using closure equations for the rotary diffusion coefficient advanced by either (i) Folgar and Tucker (J. Reinforced Plast. Comp. 3 98–119 1984) or (ii) Koch (Phys. Fluids 7(8) 2086–2088 1995). Each of these two closure models for the rotary diffusion coefficient contains an unknown empirical constant that must be determined from experiments. These were fit to experimental data along the central streamline of the contraction as a function of fibre concentration. The diffusion coefficient was found to first increase with increasing suspension concentration up to a maximum, and then decrease with concentration above this point. This non-monotonic behavior was attributed to fibre flocculation, a mechanism not considered in the relationships for the rotary diffusion coefficient. The theoretical model is then extended to predict fibre orientation over the entire plane of the contraction and the two-way momentum coupling between the fluid and fibre phases were investigated numerically. The results show that the structure of the flow field within the contraction is significantly altered when the fibre phase is considered, demonstrating the non-negligible effect of the momentum exchange between the two phases. Comparison is made between the predicted orientation state of the suspension with experimental observations over the contraction plane. Good agreement was found between the model predictions and the experimental observations except in a small region near the solid boundaries. These near wall discrepancies were attributed to an inability to correctly handle the wall boundary conditions in the fibre orientation model. |
Extent | 1862140 bytes |
Subject |
Fibre suspensions Multiphase flow Fibre orientation Fluid mechanics |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-05-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0066399 |
URI | http://hdl.handle.net/2429/842 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2008-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2008_spring_krochak_paul.pdf [ 1.78MB ]
- Metadata
- JSON: 24-1.0066399.json
- JSON-LD: 24-1.0066399-ld.json
- RDF/XML (Pretty): 24-1.0066399-rdf.xml
- RDF/JSON: 24-1.0066399-rdf.json
- Turtle: 24-1.0066399-turtle.txt
- N-Triples: 24-1.0066399-rdf-ntriples.txt
- Original Record: 24-1.0066399-source.json
- Full Text
- 24-1.0066399-fulltext.txt
- Citation
- 24-1.0066399.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0066399/manifest