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 eﬀects of long range hydrodynamic ﬁbre-ﬁbre interactions on the orientation state of a semi-dilute, rigid ﬁbre suspension ﬂowing through a linear contracting channel under laminar ﬂow conditions. The eﬀects of ﬁbre-ﬁbre interactions have been 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 ﬁbre-ﬁbre interactions are random in nature, thus giving rise to a diﬀusion-type process. The orientation distribution evolves spatially according to a Fokker-Plank type equation using closure equations for the rotary diﬀusion coeﬃcient 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 diﬀusion coeﬃcient contain an unknown empirical constant that must be determined from experiments. These were ﬁt to experimental data along the central streamline of the contraction as a function of ﬁbre concentration. In both cases the diﬀusion coeﬃcient was found to ﬁrst increase with increasing suspension concentration up ii Abstract to a maximum, and then decrease with concentration above this point. This nonmonotonic behavior was attributed to ﬁbre clumping or ﬂocculation, a mechanism not considered in the relationships for the rotary diﬀusion coeﬃcient. The theoretical model was then extended to predict ﬁbre orientation over the entire plane of the contracting channel and the two-way momentum coupling between the ﬂuid and ﬁbre phases was investigated numerically. The results show that the structure of the ﬂow ﬁeld within the contraction is signiﬁcantly altered when the ﬁbre phase is considered, demonstrating the non-negligible eﬀect 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 is found between the model predictions and the experimental observations, except in a small region near the solid boundaries where the two diﬀered signiﬁcantly. These near wall discrepancies are attributed to an inability to correctly handle the wall boundary conditions in the ﬁbre 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 Deﬁnition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4 Estimating Dr 4.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.2 Experimental Detail . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 . . . . . . . . . . . . . . . . . . . 38 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.2.1 4.3 Results 4.3.1 Numerical Estimates of Dr . . . . . . . . . . . . 45 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5 Validating Dr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.4 Comparison of Rotary Diﬀusion Models 5.1 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.2 2D Orientation Results 54 . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Theoretical Predictions . . . . . . . . . . . . . . . . . . . . . 5.2.2 Experimental Observations . . . . . . . . . . . . . . . . . . . 5.3 Discussion of Fibre Orientation Results 5.4 Eﬀect of 2 Way Coupling 55 57 . . . . . . . . . . . . . . . . 60 . . . . . . . . . . . . . . . . . . . . . . . . 62 5.4.1 Flow Field . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.4.2 Orientation Distribution . . . . . . . . . . . . . . . . . . . . . 65 6 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 68 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.2.1 Folgar-Tucker Model . . . . . . . . . . . . . . . . . . . . . . . 70 6.2.2 Koch Model 70 6.1 General Summary 6.2 Diﬀusion Coeﬃcient . . . . . . . . . . . . . . . . . . . . . . . . . . . v Table of Contents 6.3 Comparison of Diﬀusion 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 Diﬀusion Coeﬃcient . . . . . . . . . . . . . . . . . . . . . . . 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 B.2.2 Computation of the Fibre Stress . . . . . . . . . . . . . . . . . 97 . . . . . . . . . . . . . . . . 106 vi List of Figures 1.1 The channel geometry and ﬁbre orientation angles used in this study. 4.1 The orientation of a ﬁbre with respect to ﬂow in the contraction. The 2 angle φ is the angle of the projection of the ﬁbre in the xy-plane and θ is the angle between the ﬁbre 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 ﬁbre observation separated by 200 ms. The upper image pair sequentially precedes the bottom image pair. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 35 An example of an original image (left) and a post-processed image with white ﬁbres on a black background (right). The ﬁbre concentration is nL3 = 8, the ﬁbre length is 5mm, and the aspect ratio is 50. . . . . . 4.5 36 The region bounded by the red lines represents the analysis region. All ﬁbres observed in this area are used for the 1D analysis and determination of CI and λ. . . . . . . . . . . . . . . . . . . . . . . . . . 38 vii List of Figures 4.6 Development of the orientation distribution along the contraction centerline for the φ projection. x L = 0(◦), x L = 0.33( ), x L The experimental data are given at = 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. . . . . . . . . . . . . . . . . . . . . . . . 4.7 40 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% conﬁdence interval. . . . . . . . . . . . . . 41 4.8 Plot of CI (− −) and λ(− −) vs nL3 . . . . . . . . . . . . . . . . . . 42 4.9 nL −) with CI (− −) vs nL3 . . . . . Shown is a comparison of λ ln 2 r (− 3 44 4.10 Development of the orientation distribution along the contraction centerline for the φ projection. x L = 0(◦), x L = 0.33( ), x L The experimental data are given at = 0.66( ) and x L = 1.0( ). The best ﬁt 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. . . . . 5.1 Plot of the L2 norm of the relative change in the ﬂow ﬁeld between successive iterations. The concentration shown here is nL3 = 8. . . . 5.2 46 53 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 ﬁbre orientation angle. Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom). . . . . . . 5.4 Numerical predictions of the orientation variance, σ 2 . Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom). . . . . . . 5.5 59 Experimental observation of the orientation variance σ 2 . Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom) . . . . . . 5.7 58 Experimental observations of the mean ﬁbre orientation angle. Shown are nL3 = 8 (Top), nL3 = 16 (Middle) and nL3 = 24 (Bottom) . . . . 5.6 56 60 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). . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 63 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 Diﬀerence 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 eﬀect of long range hydrodynamic ﬁbre-ﬁbre interactions on the orientation state of rigid ﬁbres in a viscous ﬂuid ﬂowing through a linear contracting channel are investigated. The case considered here is for ﬂow with dilute to semidilute suspensions where the Reynolds numbers based on the length of the ﬁbre is asymptotically small and based on the inlet channel width is order one. The monodispersed ﬁbres studied here are cylindrical in shape with a length and diameter of 5 mm and 100 μm respectively; the aspect ratio of the ﬁbres, r, that is the ratio of length to diameter is 50. The concentration limits of interest in this work are deﬁned such that 1 ≤ nL3 ≤ r, where n is the number of ﬁbres per unit volume and L is the ﬁbre length. The contraction consists of two rigid sloping walls separated by 50.8 mm at the entry and 5.08 mm at the exit and the length of the contraction is 130 mm, see Figure 1. Controlling the orientation state of ﬁbres in contraction-type ﬂows is one of the principal objectives in many science and industrial processes. For example, during the electrospinning of polymer based nanoﬁbres, a polymer suspension is fed through a small contraction called a spinneret. A strong electric ﬁeld is then applied which 1 Chapter 1. Introduction Figure 1.1: The channel geometry and ﬁbre orientation angles used in this study. draws the suspension out through the contraction. This serves to align the polymer molecules in the ﬂow direction. The motivation behind this thesis however is based on a similar but much larger-scale industrial process, namely papermaking. During papermaking, a ﬂocculated ﬁbre suspension is ﬂuidized by turbulence created locally from a sudden change in geometry in an apparatus called a headbox. The ﬂuidized ﬁbre suspension is then formed into a plane liquid jet through an elongational ﬂow 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 ﬁnal ﬁbre orientation distribution. One of the essential scientiﬁc diﬃculties is the incomplete understanding of the long- 2 Chapter 1. Introduction range hydrodynamic interactions between ﬁbres in the suspension. A hydrodynamic ﬁbre-ﬁbre interaction is said to occur when one ﬁbre is swept within the disturbance ﬂow ﬁeld of another ﬁbre causing the ﬁbre to undergo a small change in orientation (and vice-versa). This eﬀect propagates throughout the suspension resulting in deviations from the expected orientation state of the suspension. Another diﬃculty in understanding these ﬂows is that momentum is exchanged between the ﬂuid and ﬁbre phases as the suspension ﬂows. Both the magnitude and direction of this momentum transfer are highly dependent on the orientation state of the suspension, a complication which eﬀectively couples the structure of the ﬂow with the orientation state of the suspension. Moreover, the momentum exchange between the two phases can aﬀect 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 ﬂocculate during processing. The ability to manipulate and control ﬁbre orientation in the headbox is critical to forming higher quality paper. One of the key elements towards this is the ability to predict the eﬀect of diﬀerent system parameters on the orientation state. This means that a highly accurate and eﬃcient numerical model must be available so that initial predictions can be made before larger-scale prototypes are designed and tested. In this work, the eﬀects of long range hydrodynamic ﬁbre-ﬁbre interactions on the orientation state of dilute to semi-dilute rigid ﬁbre suspensions ﬂowing through a linear contracting channel under laminar ﬂow conditions is investigated numerically and experimentally. The eﬀects of ﬁbre-ﬁbre interactions are modeled mathematically by assuming that the orientation state of the suspension can be completely 3 Chapter 1. Introduction described by a probability distribution function and that ﬁbre-ﬁbre interactions are random in nature, thus giving rise to a diﬀusion-type process. The spatial evolution of the orientation distribution obeys a Fokker-Plank type equation where two different relationships for the rotary diﬀusion coeﬃcient 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 ﬁbre orientation along the channel length. Also studied is the eﬀect of the two way momentum coupling on the orientation state of the suspension and on the structure of the ﬂow ﬁeld. This is accomplished by including the additional stress arising from the ﬁbre phase in the ﬂuid momentum equatons. The ultimate goal of this work is to improve the current state of ﬁbre orientation modelling in the headbox. Modelling ﬁbre motion is a formidable task. Perhaps the greatest obstacle relating to this problem is the extremely large computational domain required to solve the problem regardless of the speciﬁc 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 ﬁbre orientation model and for the momentum exchange between the ﬂuid and ﬁbre 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 diﬀusion coeﬃcient. 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 diﬀusion coeﬃcient. In §5, the theoretical model is extended to predict ﬁbre orientation over the entire xy-plane of the contracting channel. The eﬀect 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 eﬀects of the two-way coupling and discrepancies between the numerical predictions and experimental observations of ﬁbre 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 ﬂowing ﬁbre suspensions is given. In §2.1, we begin the discussion by deﬁning the traditional means of describing ﬁbre 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 eﬀect of the ﬁbres on the carrier ﬂuid and describe an additional stress imposed on the ﬂuid. In §2.3 a description 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 eﬀect of the wall on the orientation distribution. 2.1 Description of Fibre Orientation The earliest study on the motion of a single isolated ﬁbre dates back to Jeﬀerey (1922) where he derived an equation of motion for a single rigid ellipsoid in a viscous unbounded Newtonian ﬂuid. Using a no-slip boundary condition at the ﬁbre’s surface 6 Chapter 2. Literature Review and matching the velocity ﬁeld in the region near the ﬁbre to the bulk ﬂow ﬁeld of the suspending medium, Jeﬀerey solved for the ﬂow ﬁeld around the particle. A formula for the angular velocity vector of the ﬁbre was then derived through a torque balance about the centre of the ﬁbre. Jeﬀerey’s derivation showed that a ﬁbre rotates in one of a family of closed orbits around the vorticity axis termed the Jeﬀerey orbits. Bretherton (1962) showed that Jeﬀerey’s equations apply to cylindrical particles if the particle aspect ratio is replaced by an eﬀective aspect ratio. What is clear from these studies is that ﬁbers rotate continuously in response to shear within the ﬂow and that there is no preferred orientation distribution. For cases above this dilute limit, ﬁbre-ﬁbre interactions create deviations from the expected orientation state. As such, modelling the orientation state of non-dilute suspensions must be handled diﬀerently. With shear ﬂows, 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 ﬁbre orientation by a probability density function Ψ that obeys a Fokker-Planck relationship, that is 1 ∂Ψ ∂ ∂ 1 ∂Ψ DΨ ˙ + 1 ˙ = (Dr sin(θ) − sin(θ)θΨ) (Dr − sin(θ)φΨ) Dt sin(θ) ∂θ ∂θ sin(θ) ∂φ sin(θ) ∂φ (2.1) where D Dt is the total derivative with ﬂuid velocity ﬁeld components; and u and v are the velocity components in the x and y directions respectively. The ﬁbre orientation 7 Chapter 2. Literature Review vector, p, is deﬁned as a unit vector pointing parallel to to the major axis of a ﬁbre which, in three dimensional Euclidean space, is deﬁned as ⎡ ⎤ ⎢ cos φ sin θ ⎢ p=⎢ ⎢ 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 θ ⎢ − sin φ ⎥ ⎢ ⎥ ˆ = ⎢ cos φ ⎥ and φ ⎢ ⎥ ⎣ ⎦ 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ The rotational velocity can then be expressed as ˆ + θ˙θ ˆ p˙ = φ˙ sin θφ (2.3) Olson et al (2004) have derived the following expressions for the rotational velocities in terms of orientation angles φ˙ = 12 3 Lf sin θ Lf /2 −Lf /2 ˆ · u(y + lp, t)dl lφ (2.4) 8 Chapter 2. Literature Review and 12 θ˙ = Lf 3 Lf /2 −Lf /2 ˆ · u(y + lp, t)dl lθ (2.5) In application to headbox ﬂows where ﬁbre length is small in comparison to the scale of the headbox, it is reasonable to locally linearize the ﬂow ﬁeld by taking a Taylor series expansion of the velocity ﬁeld around the ﬁbre centre, i.e., u(y + lp) ≈ u(y) + ∂u(y) lp ∂y (2.6) and substitute it into Equations 2.4 and 2.5 above. This results in the following simpliﬁed equations, ∂u(y) θ˙ = θˆ p ∂y (2.7) and φ˙ = φˆ sin(θ) ∂u(y) p ∂y (2.8) or expressing φ˙ in terms of its component partial derivatives for planar ﬂow where u = u(u(x, y), v(x, y)), we have 1 φ˙ = 2 ∂v ∂u − ∂y ∂x sin(2φ) − ∂v ∂u 2 sin (φ) + cos2 (φ) ∂y ∂x (2.9) and 9 Chapter 2. Literature Review ∂u ∂u θ˙ = cos2 (φ) cos(θ) sin(θ) + sin(φ) cos(φ) sin(θ) cos(θ) ∂x ∂y ∂v ∂v sin2 (φ) sin2 (θ) + cos(φ) sin(φ) sin2 (θ) + ∂x ∂y (2.10) These equations are identical to those derived by Jeﬀerey (1922) for ﬁbres of large aspect ratio. The use of an orientation distribution function to describe the orientation state of a ﬁbre 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 ﬂow ﬁeld, but will also undergo random ﬂuctuations in orientation due to Brownian rotary motion. Hence, they modeled the orientation state of polymer suspensions using a Fokker-Plank equation with a rotary diﬀusivity term, that is, Dr in Equation 2.1, which results from the eﬀects of Brownian diﬀusion 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 diﬀusivity term in their model was not due to Brownian eﬀects, but rather was assumed to be the result of random hydrodynamic interactions with neighboring ﬁbres; that is, hydrodynamic ﬁbre-ﬁbre interactions which cause ﬁbres to deviate from their Jeﬀerey 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 ﬂow in a linear contraction, Dr can be expressed as follows Dr = CI α˙ (2.11) where CI is traditionally called the interaction coeﬃcient and is related to a number of suspension parameters such as concentration, aspect ratio, and ﬁbre 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 ﬁbres in compression moldings and showed that the model compares well with experiments once the interaction coeﬃcient is determined. In studies by VerWeyst & Tucker (2002), the FT model for the rotary diﬀusion coeﬃcient was used to investigate the eﬀect of the two way coupling between the ﬂuid and ﬁbre phases for ﬂows 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 coeﬃcient, CI , to the average interparticle spacing in non-dilute suspensions. They found a strong dependence between CI and the orientation state of the suspension. More speciﬁcally, by ﬁtting the experimental data of Folgar (1983) for nL3 = 52.2 to Equation 5.1, where n is the number of ﬁbres per unit volume and L is the ﬁbre 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 signiﬁcant when nL3 = 26. What these studies did show was an orientation state dependence of the interaction coeﬃcient. Koch (1995) proposed a more physically complete form for the diﬀusion coeﬃcient 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 diﬀusion coeﬃcient 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 (λ1 IE : pppp : E + λ2 E : pppppp : E) α˙ ln2 r (2.12) where E is the ﬂuid strain rate tensor, deﬁned as 1 E = (∇u + ∇T u) 2 (2.13) and r is the aspect ratio of the ﬁbre, that is, the ratio of ﬁbre length to diameter. The remaining term to be deﬁned 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 deﬁned as pppp = pi pj pk pl Ψ(φ)sin(θ)dθdφ (2.14) The ﬁrst term on the right hand side in Equation 2.12 represents the isotropic rotary 12 Chapter 2. Literature Review diﬀusion coeﬃcient that is proportional to the average of the square of the force per unit length on ﬁbres. The second is an anisotropic correction term representing diﬀusion in the directions parallel to the force per unit length exerted by the other ﬁbres. λ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 veriﬁed for other ﬂow ﬁelds, such as for the ﬂow 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 ﬁnd that the Folgar-Tucker model is found in most studies of laminar suspension ﬂows, (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 ﬂuid and ﬁbre phases in a ﬂowing suspension has been known for some time now. Batchelor (1970) developed a general constitutive equation for the bulk stress in a suspension of rigid, inertialess particles of arbitrary shape in a Newtonian ﬂuid. 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 resul1 A Stokeslet is deﬁned as a singularity in a Stokes ﬂow which represents the resultant eﬀect of a force applied to a ﬂuid 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 ﬂuid on a cylindrical ﬁbre 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 ﬂow ﬁeld. He further showed that, to a ﬁrst approximation the drag tensor is independent of ﬁbre position and ﬂow ﬁeld. Khayat & Cox (1989) extended this result to small yet ﬁnite Reynolds number. Dinh & Armstrong (1984) extended Batchelor’s theory to account for the orientation state of elongated particles and its eﬀect 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 ﬂow ﬁeld around the particle they were able to equate Batchelor’s constitutive equation to a new constitutive equation involving Ψ via a fourth order orientation moment tensor. By doing so they coupled the ﬁbre orientation equations to the ﬂuid momentum equations. Bibbo, Dinh & Armstrong (2001) validated the work of Dinh & Armstrong (1984) by experimentally measuring the shear stress within a semi-concentrated suspension undergoing simple shear ﬂow. Their measurements demonstrated a strong coupling between ﬁbre alignment and the resulting shear stress within the suspension and that the Dinh-Armstrong model was suitable for predicting the additional ﬁbre stress. Shaqfeh, & Fredrickson (1990) used the Dinh-Armstrong model to derive an asymptotic form for the eﬀective viscosity of a semi-dilute suspension of rods. The resulting relationship for the additional stress due to the ﬁbre phase is as 14 Chapter 2. Literature Review follows τ f ibre = μcr 2 E : pppp ln(1/c) + ln(ln(1/c)) + 1.439 (2.15) where c is the volume fraction of ﬁbres within the suspension. With Equation 2.15 it becomes possible to simultaneously compute the fully coupled ﬁbre orientation distribution and the resulting ﬂow ﬁeld. 2.3 Numerical Investigations More recently, researchers have been making great eﬀorts to perform 3-D ﬁbre orientation predictions for 2 and 3D ﬂows inside complex geometries, as opposed to restricting their analysis to simple geometries only (e.g. Lin & Zhang (2002), Lipscomb & Denn (1988), VerWeyst & Tucker (1999, 2002). The diﬃculty 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 dimensional 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 ﬁne enough mesh for modelling many industrial suspension ﬂows, particularly headbox 15 Chapter 2. Literature Review type ﬂows. Despite the ongoing closure issues discussed in §2, there have still been many numerical studies that incorporate the two-way coupling in conjunction with a Fokker-Plank equation to predict the orientation state of ﬁbre suspensions in various complex geometries. 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 ﬂow of dilute ﬁbre suspensions through various axi-symmetric contractions using the ﬁnite element method (FEM). They showed that the suspension ﬂow diﬀers qualitatively from the ﬂow of the suspending ﬂuid alone. They also showed that their numerical results agreed well with experimental observation. VerWeyst & Tucker (1999) used a Galerkin ﬁnite element method to predict ﬁbre orientation patterns in a variety of injection molded features using the Hele-Shaw ﬂow model for non-dilute suspension 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 qualitative agreement with experimental observations. VerWeyst & Tucker (2002) used the fully coupled, Galerkin FEM solution to predict the ﬂow of ﬁbre suspensions through a variety of complex geometries which include axi-symmetric contractions, expansions and center-gated disks for a number of diﬀerent suspension concentrations 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 ﬁbre concentration. For Re ≥ 5, the size of the vortex decreased as the ﬁbre concentration was increased. For the center-gated disk, they report that the 16 Chapter 2. Literature Review eﬀect of the coupling was to displace streamlines toward the bottom surface as ﬂow enters the disk region, and to align the ﬁbres more rapidly in the radial direction. What is clear from these studies is that there is a signiﬁcant diﬀerence between the coupled and uncoupled solution for the ﬂow ﬁeld 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 ﬁbre suspension in a planar channel ﬂow. They ignored ﬁbre-ﬁbre interactions but used the rotary diﬀusion term to model the eﬀect of turbulent ﬂuctuations of ﬁbre orientation. They showed that their numerical predictions compared well with experimental ﬁndings in 2-D. In studies by Paschkewitz et al. (2004), the drag reduction induced by the presence of rigid ﬁbres in a turbulent ﬂow was studied numerically. There work showed that drag reductions of up to 26% were theoretically possible for rigid, non-Brownian ﬁbre suspensions in the semi-dilute concentration regime. 2.4 Experimental Investigations While there is an abundance of literature on the theoretical aspects of ﬁbre orientation, 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 ﬂow and then compared their results to the theoretical predictions of Jeﬀerey (1922). These experiments were carried out for the projection of the azimuthal angle only, that is the ﬁbre 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 ﬁbres in a semidilute, 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), ﬁbres continue to rotate about the vorticity axis and remain highly aligned in the ﬂow direction. Their measured orbital time constants were however found to be considerably diﬀerent from those obtained by Anczurowski & Mason (1968), with the diﬀerences being attributed to the rotary diﬀusion resulting from long range hydrodynamic ﬁbre-ﬁbre 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 ﬁbres in an extensional ﬂow. The mean square orientations were used to quantify the orientational diﬀusion resulting from ﬁbre-ﬁbre interactions, as well as to determine the eﬀect of concentration on this diﬀusion-type process. Their results showed a clear increase in rotary diﬀusion with an increase in ﬁbre concentration. 2.5 Previous Headbox Studies Within the pulp and paper literature several authors have investigated the orientation behavior of pulp ﬁbres in headbox ﬂows. In studies by Ullmar (1998) and Ullmar & Norman (1997), the ﬁbre orientation distribution in the plane of the paper was measured at the headbox exit by digitally imaging nylon ﬁbres in the suspension. These studies examined the eﬀect of mean ﬂow through the headbox, the headbox 18 Chapter 2. Literature Review contraction ratio, and ﬁbre concentration on the orientation state of the suspension. In another study by Zhang (1998) the ﬁbre 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 ﬂow ﬁeld along the centerline of a linear contraction, neglecting the two way coupling between the ﬂuid and ﬁbre phase, then used the Fokker-Plank equation to predict the 2-D orientation state for the resulting ﬂow. They chose a global rotary diﬀusion coeﬃcient which models the eﬀect of turbulence on the orientation state by ﬁtting their computations to the experimental data of Ullmar & Norman (1997) and Zhang (1998). Hyensj¨o et al. (2007) extended the work of Olson et al (2004) by using single phase CFD modeling to predict the ﬂow ﬁeld along streamlines of a linearly contracting channel. They used these results to compute the 2-D ﬁbre orientation distribution along individual streamlines of the ﬂow in the absence of ﬁbre-ﬁbre interactions and ﬂow-ﬁbre coupling. Parsheh et al. (2005) performed experimental studies to better understand the eﬀect 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 distribution at a set of discrete points along the centerline of the contraction and ﬁnally compared these results to predictions based on a Fokker-Plank type equation. They showed that the rotational diﬀusion coeﬃcient 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 ﬁbre orientation problem. That is, how do ﬁbres behave near solid boundaries. This question has seen some research in the past years yet there have been no deﬁnitive answers. There are several possible orientation patterns for which a ﬁbre may follow near a wall. The ﬁbre could possibly align itself in the direction of the wall. The ﬁbre could seemingly rotate, or ﬂip about the wall. The ﬁbre could shift away from the wall in response to the large ﬂow gradients in this region or orient itself perpendicular to the wall. Stover, & Cohen (1990) were the ﬁrst to rigorously address this issue. In their studies, the motion of rodlike particles in a low Reynolds number plane Poiseuille ﬂow were observed experimentally and the eﬀect of the wall on the period of ﬁbre rotation was determined. They found that when particles with a high Jeﬀery orbit constant, (that is particles with a large period of rotation about the vorticity axis) came within a distance less than half a ﬁbre length from the wall, an irreversible interaction between the ﬁbre and wall occurred whereby the ﬁbre was said to ”pole-vault” away from the wall to a distance of approximately half a ﬁbre length. This interaction was accordingly named as it seemingly mimicked the ﬂipping motion of pole-vaulter. Yet the ﬁbre was said to never actually touch the wall but rather that the ﬁbre tip would come within approximately one ﬁbre diameter from the wall. They suggested the existence of some nonhydrodynamic force between the ﬁbre and wall but could not determine the exact nature of this interaction. They observed quite a diﬀerent behavior for ﬁbres with a low Jeﬀery orbit constant that came within a half ﬁbre length from the 20 Chapter 2. Literature Review wall. In this case, the ﬁbre remained close to the wall indeﬁnitely, however if the ﬁbre had an increased period of rotation of approximately 50% it shifted its orientation away from the wall. In studies by Holm & S¨oderberg (2007) the inﬂuence of near wall shear stress on ﬁbre orientation was investigated experimentally. They measured the resulting wall plane ﬁbre orientation distribution of dilute suspensions in laminar ﬂow 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 ﬁbres were well aligned in the ﬂow direction at distances ≥ 3 mm from the wall, regardless of ﬁbre length and aspect ratio. However, very close to the wall, that is at a distance of 1 mm from the wall they found many ﬁbres that were not aligned in the ﬂow 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 ﬂow direction and the other, more dominant mode at 90o relative to the stream direction. By increasing the ﬁbre 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 ﬂowing suspension can be realized by solving a Fokker-Plank equation, that is, Equation 2.1 for the orientation distribution function, Ψ, where the eﬀects of hy21 Chapter 2. Literature Review drodynamic ﬁbre-ﬁbre interactions are modeled with the rotary diﬀusion coeﬃcient. Two models are available for the rotary diﬀusion coeﬃcient, 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. Furthermore the ﬂow ﬁeld is aﬀected by the orientation state of the suspension creating a two way coupling between the ﬂuid ﬂow and the orientation state of the ﬁbre phase. Experimental works in this area are far less evident, particularly those pertaining to contraction-type ﬂows and the correct description of the near wall orientation behavior of the suspension is still uncertain. 22 Chapter 3 Problem Deﬁnition The description of the orientation state of rigid ﬁbre suspensions in a variety of diﬀerent ﬂow 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 ﬁbre-ﬁbre interactions. Two models for the rotary diﬀusion resulting from ﬁbre-ﬁbre 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 eﬀects of ﬁbre-ﬁbre interactions on the orientation state of suspensions subject to ﬂow within linear contractions, that is, for ﬂow within geometries typical of modern headboxes are non-existent. This leads to the ﬁrst set of objectives for this research: 1. To experimentally measure the rotary diﬀusion coeﬃcient as a function of ﬁbre 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 diﬀusion coeﬃcient. That is, the models of Folgar & Tucker (1984) and Koch (1995). 3. To compare the utility of these two models for predicting the suspension orientation state by comparison with experimental measurements. With the diﬀusion coeﬃcient 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 ﬁbre orientation, but also for the further study of tensor-based closure approximations in contraction type ﬂows. Also missing in the literature are investigations into the eﬀects of the two way momentum coupling between the ﬂuid and ﬁbre phases on the suspension orientation state, and on the structure of the ﬂow ﬁeld for geometries similar to that of a headbox. Studies found in the literature have shown that the eﬀect of the two way coupling is very signiﬁcant in simpler ﬂow geometries, however none have addressed its eﬀect on suspension ﬂow 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 eﬀects of 2 and 3 dimensional ﬂow gradients on ﬁbre alignment. This leads to the second set of objectives for this thesis: 1. To directly compute the 2D fully coupled planar orientation distribution function over the entire plane of the contraction. 24 Chapter 3. Problem Definition 2. To investigate the eﬀects of the two way coupling on both the ﬂow ﬁeld and orientation state of the suspension. Lastly, the numerical predictions of the suspension orientation state will be compared with 2D observations of ﬁbre 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 diffusion coeﬃcient, 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 represents 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 procedure 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 ﬁbres in suspension are the Lagrangian and Eulerian methods. In the Lagrangian method, the equations of motion for a single ﬁbre are solved given a known velocity ﬁeld and the process is repeated for each ﬁbre within the suspension. While the Lagrangian method can be quite accurate, it is computationally intensive, particularly for non-dilute suspensions ﬂowing within complex geometries. With the Eulerian method, the probability distribution of ﬁbre orientation (and possibly position) is determined by solving a convection-diﬀusion equation, namely a Fokker-Plank type equation. With this approach the mean ﬂow ﬁeld convects the ﬁbres position and orientation while long range, hydrodynamic ﬁbre-ﬁbre interactions diﬀuses ﬁbre orientation. This diﬀusion type process eﬀectively results in a ﬂux which opposes the local mean orientation state of the ﬁbres. The following assumptions will be made for the Eulerian analysis outlined below: 1. The nominal ﬁbre length is small enough that the ﬂow ﬁeld can be assumed to be linear along the length of the ﬁbre. 2. The orientation state of the ﬁbres at a point in space can be described by a probability distribution function, Ψ, where Ψ is deﬁned such that the probability of ﬁnding a ﬁbre oriented between the angles φ and φ + ∂φ is Ψ(φ)∂φ. 3. The ﬁbres are uniformly distributed in space throughout the ﬂow ﬁeld. 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 FokkerPlank equation is solved, that is, the steady-state form of Equation 2.1 u ∂Ψ 1 ∂ 1 ∂ 1 ∂Ψ ∂Ψ ∂Ψ ˙ ˙ +v = (Dr sinθ −sinθθΨ)+ (Dr −sinθφΨ) (4.1) ∂x ∂y sinθ ∂θ ∂θ sinθ ∂φ sinθ ∂φ where u and v are the ﬂuid 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 diﬀusion coeﬃcient which is given by either Equations 2.11 or 2.12. It should be noted that in this work we deal only with the ﬁrst term from Equation 2.12, that is the isotropic component of the rotary diﬀusion coeﬃcient. For applications to accelerating ﬂows the second term in Equation 2.12 provides a lower order correction to the isotropic part of the diﬀusion coeﬃcient. Therefore the marginal improvement in accuracy of the rotary diﬀusion coeﬃcient does not justify the added computational cost of working with this term. Hence we deﬁne Koch’s rotary diﬀusion coeﬃcient as follows Dr = nL3 (λIE : pppp : E) α˙ ln2 r (4.2) A simpliﬁcation of the model equations can be realized by considering the ﬂow of ﬁbres 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 ﬁbre with respect to ﬂow in the contraction. The angle φ is the angle of the projection of the ﬁbre in the xy-plane and θ is the angle between the ﬁbre 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 − R1 ) x L , v(x) = 0 (4.3) where R is the contraction ratio of the headbox deﬁned as the ratio of inlet to exit heights which for this analysis is ﬁxed at a value of 10, Equation 2.1 reduces to 29 Chapter 4. Estimating Dr u ∂Ψ ∂Ψ 1 ∂ ˙ + 1 ∂ (Dr 1 ∂Ψ − sinθφΨ) ˙ = (Dr sinθ − sinθθΨ) ∂x sinθ ∂θ ∂θ sinθ ∂φ sinθ ∂φ (4.4) Here the projected rotational velocities along the central streamline, that is φ˙ and θ˙ reduce to ∂u sin(2φ) φ˙ = ∂x (4.5) ∂u sin(φ)cos(φ)cos2 (θ) θ˙ = ∂x (4.6) and The planar orientation distribution function, Ψp can then be extracted from the 3D orientation distribution function as follows. π Ψp = Ψ(φ, θ)sin(θ)dθ (4.7) 0 The boundary conditions on Ψ are as follows. Since one end of a ﬁbre is indistinguishable 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 π Ψ(x, φ, θ)sin(θ)dθdφ = 1 (4.9) 0 At the inlet, we set Ψ equal to the experimentally observed inlet orientation distribution, 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 distribution 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 ﬁbre suspensions it becomes increasingly diﬃcult to distinguish individual ﬁbres from one another, particularly as the ﬁbre 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 ﬂuid and ﬁbre phases with identical, or nearly 31 Chapter 4. Estimating Dr identical indexes of refraction. The result is that the ﬁbre phase becomes indistinguishable from the ﬂuid phase when observed under white light; loosely speaking, the ﬁbres become invisible. A small number of ﬁbres (less than 1% of the total number of ﬁbres) are then coloured so that they can be visualized and their motion and orientation are observed and taken to represent the behaviour of all ﬁbres within the suspension. For these experiments, cylindrical boroscilicate glass ﬁbres (www.mosci.com) of dimensions 5 × 0.1 mm were employed as the ﬁbre phase. The glass ﬁbres have a density of approximately 2250 kg/m3 . The index of refraction of these ﬁbres was measured commercially and found to be 1.4719(±0.0005). Approximately 0.01% of the ﬁbres 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 ﬁbres were washed in detergent, rinsed in alcohol and then in distilled water. After silvering the ﬁbres were again washed in a similar manner to remove any loosely adsorbed AgNO3 . We measure the motion and orientation of 4 diﬀerent monodispersed suspensions of concentration, nL3 = 8, 16, 24 and 32. The Newtonian ﬂuid used in this was glycerine with a density of 1260 kg/m3 and a viscosity of 1.49 P a 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× 130 mm) with approximately 1.5 litres of suspension. The contraction consists of two rigid sloping walls, separated by 50.8 mm at the entry and 5.08 mm at the exit; the contraction ratio is R = 10. The distance between the entry and the exit is 130 mm, 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 diﬀerent heights to control the pressure drop over the cell. The ﬂowrate was set at 0.0021 m3/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 Masterﬂex B/T variable speed peristaltic pump (www.coleparmer.com). The suspension was gently mixed prior to entering the channel in an eﬀort to create a more uniform orientation distribution at the inlet while not disturbing the ﬂow 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 ﬁbre trajectories through the middle of the channel and minimize any inﬂuence of the solid walls 33 Chapter 4. Estimating Dr on ﬁbre orientation. Furthermore, the width of the channel is large relative to the length of the ﬁbre so we neglect any inﬂuences of the wall boundary layers on the ﬁbre 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 ﬁrst 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 ﬁbre orientation in the xz-plane were used only to initialize the orientation distribution function at the contraction inlet when measuring the diﬀusion coeﬃcient. With the optical setup described the imaged volume is 130×65×5.1 mm3 (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 ﬁeld of the lens. The plexiglass cell was transilluminated using Schott-Fostec ﬁbre optic dual backlight. The images were split between the two projection planes after which a mask was applied to isolate the actual ﬂow region from the remainder of the image. The images were then saved using a Matrox Meteor II digital frame grabber, and the ﬁbre tracking analysis was performed on the xy-plane. The orientation of each particle was calculated using an in-house ﬁbre 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 xy−plane xz−plane xz−plane Figure 4.3: Two consecutive unprocessed images of ﬁbre observation separated by 200 ms. 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 ﬁbres on a black background (right). The ﬁbre concentration is nL3 = 8, the ﬁbre 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 ﬁlled with white pixels, and then eroded back to their original size. The image now contains ﬁlled 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-ﬁbrous debris) and are removed from the image. 4. Orientation Distribution - The orientation angle of each ﬁbre is measured relative to the horizontal by computing the arctangent of the ratio Δy Δx between the end points of a ﬁbre. The associated ﬁbre length and the position of the center of area is also recorded. 5. The ﬂow region is partitioned into 5 mm × 5 mm cells and the orientation angle of each ﬁbre whose center of area lies within a particular cell is computed. In total we ﬁt 26 cells in the horizontal direction along the 130 mm length contraction centerline. There were on average approximately 5000 tracer ﬁbre observations in each cell for each experiment. 37 Chapter 4. Estimating Dr 4.2.1 Numerical Estimates of Dr The diﬀusion coeﬃcient is determined from observations of ﬁbre orientation in a region near the central streamline. More speciﬁcally, the region of analysis for these particular experiments is deﬁned 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 ﬁbre 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 ﬁbre (see experimental detail section for full description of binning tracer ﬁbres into spatial cells). To reiterate, there are 26 cells along the central streamline, with a length of 5 mm each. Analysis Region Figure 4.5: The region bounded by the red lines represents the analysis region. All ﬁbres 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 empirical 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 diﬀerence scheme in both orientation angles φ and θ and a ﬁrst-order accurate forward diﬀerence scheme in x. In this case the velocity ﬁeld along the central streamline is assumed to be of the form u = (u(x), 0) with the boundary conditions given by Ψ(x, φ, θ) = Ψ(x, φ + π, π − θ) π 2 − π2 (4.11) π Ψ(x, φ, θ)sin(θ)dθdφ = 1 (4.12) 0 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 diﬀusion model is being used. 3. Equation 4.4 is solved using a Gauss-Seidel method. 4. When solving using Koch’s deﬁnition, 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 diﬀerence between the model and experimental prediction. 39 Chapter 4. Estimating Dr 4.3 Results To help understand the main ﬁndings of this section, it is instructive to discuss qualitatively 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 diﬀerent axial positions. Each distribution is centered around zero, that is they are aligned in the direction of the ﬂow, and the spread in the distribution decreases with increasing axial positions. These results are representative of all cases tested. nL3 = 8 7 6 5 4 Ψ 3 2 1 0 −1.5 −1 −0.5 0 φ 0.5 1 1.5 Figure 4.6: Development of the orientation distribution along the contraction centerline for the φ projection. The experimental data are given at Lx = 0(◦), Lx = 0.33( ), x = 0.66( ) and Lx = 1.0( ). The concentration shown here is nL3 = 8 and the L 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 ﬁgure, a low value of σ 2 corresponds to highly aligned orientation 40 Chapter 4. Estimating Dr 0.08 0.07 0.06 0.05 σ2 0.04 0.03 0.02 0.01 0 0 0.1 0.2 0.3 0.4 0.5 x 0.6 0.7 0.8 0.9 1 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% conﬁdence interval. state, and a high value of σ 2 corresponds to a more random, or isotropic orientation state. The ﬁrst observation that can be made from Figure 4.7 is that the diﬀerence in σ 2 at the outlet is small, but nonetheless unique for each concentration. Conversely, at the inlet we see that the variance is signiﬁcantly 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 ﬁndings of this section where CI and λ are determined as a 41 Chapter 4. Estimating Dr function of ﬁbre 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 signiﬁcant role in this region for ﬁbre 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. −3 5 x 10 0.01 CI 4 3 0.005 λ 5 10 15 20 25 30 0 35 3 nL Figure 4.8: Plot of CI (− −) and λ(− −) vs nL3 . As shown in this ﬁgure CI initially increases monotonically with concentration. This result was expected given that higher concentrations lead to an increase in the frequency of ﬁbre-ﬁbre interactions, hence greater diﬀusion 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 coeﬃcient. It is believed that at such a high concentration the orientation model breaks down due to the assumption of purely hydrodynamic interactions. More speciﬁcally, it is believed that at higher concentrations the eﬀects of ﬁbre ﬂocculation and mechanical contact interactions between ﬁbres has a hindering eﬀect on ﬁbre orientation. This observation further advances the argument that ﬂocculation or mechanical entanglement is the cause for such a condition. Martinez et al (2001) have demonstrated that for papermaking ﬁbres, ﬁbre mobility diminishes through ﬂocculation at nL3 = 31. This phenomenon is compounded with the aligning eﬀects of the accelerating ﬂow inside the contraction by which the suspension is forced into a highly aligned state and ﬁbre 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 ﬁbre volume fractions, c, in the range 0.0004 − 0.16 and with a constant ﬁbre 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 diﬀerence is that our lowest value of CI corresponds to the highest concentration. Another probable cause for the diﬀerences, in addition to the fact that we have used a diﬀerent ﬂow 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 ﬁbre suspensions. They report simulated values 43 Chapter 4. Estimating Dr −3 5 −3 x 10 x 10 4 λ nL3/ln2r 4 C I 3 3 5 10 15 20 25 30 2 35 3 nL 3 nL Figure 4.9: Shown is a comparison of λ ln −) with CI (− −) vs nL3 . 2 r (− 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 ﬁbre aspect ratio, that is the ratio of length to diameter. The trend for λ is also interesting. We see that λ decreases monotonically with increasing ﬁbre concentration, something which is unexpected since the eﬀects of ﬁbre concentration should be accounted for in the interaction rate term in the diﬀusion coeﬃcient, that is, the nL3 part of Equation 4.2. This ﬁnding would suggest that the concentration dependence in Equation 4.2 is not quite correct, at least for a 3 nL linearly contracting channel. What is interesting, however, is to plot the λ ln 2 r term in Equation 4.2, alongside CI . This is shown in Figure 4.9. This plot shows several interesting points. The ﬁrst observation is that the diﬀusion coeﬃcient does scale with the velocity gradient α˙ for this one-dimensional ﬂow. 44 Chapter 4. Estimating Dr The second key observation is that the term in Equation 4.2 containing the fourth order orientation tensor eﬀectively acts as a ampliﬁcation factor to the diﬀusion coeﬃcient for each concentration. 4.3.1 Comparison of Rotary Diﬀusion Models The evolution of Ψ is now characterized for the experimental conditions tested. An example of the development of both the predicted and observed orientation distribution 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 different equations for the rotary diﬀusion coeﬃcient. The histograms of the orientation probability functions for each axial position were found to be smooth and clustered around a deﬁnite 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 ﬁt 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 nL3 = 8 3 nL = 8 7 7 6 6 5 5 4 4 Ψ Ψ 3 3 2 2 1 1 0 −1.5 −1 −0.5 0 φ 0.5 1 1.5 0 −1.5 −1 −0.5 3 6 6 5 5 4 1.5 0.5 1 1.5 0.5 1 1.5 0.5 1 1.5 4 Ψ Ψ 3 3 2 2 1 1 −1 −0.5 0 φ 0.5 1 1.5 0 −1.5 −1 −0.5 3 0 φ 3 nL = 24 nL = 24 7 7 6 6 5 5 4 4 Ψ Ψ 3 3 2 2 1 1 −1 −0.5 0 φ 0.5 1 1.5 0 −1.5 −1 −0.5 3 0 φ 3 nL = 32 nL = 32 7 7 6 6 5 5 4 4 Ψ Ψ 3 3 2 2 1 1 0 −1.5 1 3 7 0 −1.5 0.5 nL = 16 nL = 16 7 0 −1.5 0 φ −1 −0.5 0 φ 0.5 1 1.5 0 −1.5 −1 −0.5 0 φ Figure 4.10: Development of the orientation distribution along the contraction centerline for the φ projection. The experimental data are given at Lx = 0(◦), 46 x x x = 0.33( ), L = 0.66( ) and L = 1.0( ). The best ﬁt of Equation 4.4 is given L 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. 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 ﬁbre suspension ﬂowing through a linearly contracting channel. A simpliﬁed 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 diﬀusion coeﬃcient were given. Numerical estimates for the unknown closure parameters in the two models for the rotary diﬀusion coeﬃcient have been presented, and the utility of the two models evaluated by comparison with experiments. It was shown that the interaction coeﬃcient, that is, CI in Equation 2.11 ﬁrst increases linearly with concentration 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 diﬀusion coeﬃcient breaks down and that this breaking point is the result of mechanical interactions and ﬁbre ﬂocculation. 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 ﬁnding would suggest that the scaling term in 4.2 is not quite correct for this particular type of suspension ﬂow. A comparison of the two models for the rotary diﬀusion coeﬃcient 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 diﬀusion coeﬃcient were made in §4. These estimates were based on observations of ﬁbre orientation along the central streamline of the contraction. In this section, we use these estimates for Dr and CI to predict ﬁbre orientation over the entire xy-plane of the contraction and include the two-way coupling between the ﬂuid and ﬁbre 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 ﬁbre orientation model breaks down at such a high concentration due to mechanical interactions and ﬂocculation eﬀects. 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 diﬀerences. In §5.4 the eﬀects of the two way coupling on the resulting ﬂow ﬁeld 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 ˙ ∂ 2 Ψ ∂(φΨ) ∂Ψ ∂Ψ +v = Dr 2 + ∂x ∂y ∂φ ∂φ (5.1) ˙ was shown in §2 to be independent of θ and therefore The orientation velocity, φ, assumes that same form as that given by Equation 2.9, that is 1 φ˙ = 2 ∂v ∂u − ∂y ∂x sin(2φ) − ∂v ∂u 2 sin (φ) + cos2 (φ) ∂y ∂x (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 orientation distribution, that is Ψ(0, y, φ) = Ψexp (0, y, φ) (5.5) The ﬂuid ﬂow is described using Cauchy’s equations of motion for an incompressible ﬂuid, that is ∇·u= 0 ρ ∂u + u · ∇u ∂t = −∇P + ∇ · τ (5.6) (5.7) where ρ is the ﬂuid density, P is the pressure and τ is the stress tensor, that is the sum of both the Newtonian ﬂuid and ﬁbre contributions τ = μ(∇u + ∇uT ) + τ f ibre (5.8) where τ f ibre 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 ﬂow 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. , ∂v , ∂u , ∂v are used to directly com1. For each streamline, the values of u, v, ∂u ∂x ∂x ∂y ∂y ˙ y) according to Equation 2.9 at each streamline position, s(x, y). pute φ(x, 2. Equation 5.1 is then mapped onto the streamlines according to the following transformations ∂Ψ ∂s ∂Ψ = ∂x ∂s ∂x (5.9) ∂Ψ ∂Ψ ∂s = ∂y ∂s ∂y (5.10) This leads to the steady-state, quasi 1D form of 5.1 u ˙ ∂Ψ ∂s ∂ 2 Ψ ∂(φ(s(x, y))Ψ) ∂Ψ ∂s +v = Dr 2 + ∂s ∂x ∂s ∂y ∂φ ∂φ (5.11) 3. Equation 5.11 is discretized using forward diﬀerences with respect to the spatial streamline variable, s, and centered diﬀerences with respect to the orientation angle φ. 4. The ﬂow ﬁeld is computed in this device using a 2D commercial software package (FLUENT) with the gradient of the additional ﬁbre stress term treated as a momentum source. An iterative procedure is used whereby the ﬂow ﬁeld is initially determined for the pure ﬂuid alone (i.e. water with no ﬁbres present), 52 Chapter 5. Validating Dr after which Equations 5.1 and 5.2 are solved using the initial ﬂow ﬁeld data. The gradient of the ﬁbre stress is then computed based on the initial solution for Ψ, which is then fed back into the ﬂuid ﬂow equations, which are then solved again to produce a new ﬂow ﬁeld, 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 ﬂow ﬁeld 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 ﬂow ﬁeld is shown in Figure 5.1 below and is typical of all concentrations. nL3 = 8 0 10 −1 10 −2 ||Ui − Ui−1|| 10 −3 10 −4 10 −5 10 −6 10 −7 10 1 1.5 2 2.5 iteration 3 3.5 4 Figure 5.1: Plot of the L2 norm of the relative change in the ﬂow ﬁeld 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 diﬀerence over the contraction is given for solutions obtained with the 736 x 500 mesh and another solution obtained with a 1397 x 1000 mesh. The diﬀerence in mesh solutions over the channel is seen to be small with the greatest diﬀerence occurring near the inlet and away from the central streamline region. The greatest diﬀerence 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 diﬃcult to smooth with the numerical solver. This diﬀerence does considerably diminish further down the channel and is not believed to signiﬁcantly 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 numerical 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¯ and variance, σ 2 , of the orientation tire 2D contraction plane we use the mean, φ, 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 deﬁned in the traditional manner as follows π 2 φ¯ = 2 σ = φΨ(x, y, φ)dφ (5.12) ¯ 2 Ψ(x, y, φ)dφ (φ − φ) (5.13) −π 2 π 2 − π2 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 eﬀect of ﬁbre concentration on the mean orientation angle is somewhat expected since ﬁbre-ﬁbre interactions are assumed to randomize ﬁbre orientation but should not aﬀect 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 ﬁbre 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 ﬁbre 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 ﬁbre 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 reﬂects the fact that ﬁbres align in response to the large ﬂow gradients near the contraction exit. Near the walls however we see an increase in 56 Chapter 5. Validating Dr the isotropy of ﬁbre orientation, that is a high value of variance, suggesting that the walls eﬀectively randomize ﬁbre orientation. It is also interesting to note the large development region near the inlet where ﬂow gradients are much smaller. It can be interpreted that in this region, ﬁbre-ﬁbre interactions have a much more dominant eﬀect and help maintain the orientation isotropy of the suspension. The eﬀect of increasing ﬁbre concentration is reﬂected 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 ﬁbre-ﬁbre interactions, and ultimately a greater diﬀusion of orientation and a more isotropic orientation state. This eﬀect is most obvious in the earliest stages of the contraction where the ﬂow gradients are still relatively small. However, the eﬀects of the ﬁbre-ﬁbre 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 ﬁbre orientation, orientation variance, and orientation distribution at a set of 5 mm x 5 mm spatial grid cells ﬁlling the entire section of the channel. Figure 5.5 shows the mean ﬁbre orientation observed experimentally. The ﬁrst observation 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 ﬁbres 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 ﬁbre orientation in the wall layer region. That is, the experiments show that ﬁbres align with the wall in the wall layer whereas the model predicts a shift in ﬁbre 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 ﬁbre 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 predictions. That is, we see the greatest values of σ 2 near and around the inlet, and this eﬀect propagates a ﬁnite distance down the channel where the ﬂow gradients are lower. Near the exit, σ 2 becomes very small indicating a highly aligned suspension. The eﬀect of concentration is also similar to the model predictions, that is, increasing the ﬁbre concentration increases the size of the development region. Near the walls however, the experiments diﬀer 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 ﬁbres has received some attention over the years. For example, Hyensj¨o et al. (2007) modeled the inﬂuence of shear on ﬁbre orientation in turbulent headbox ﬂows. They showed that near the walls, the 60 Chapter 5. Validating Dr tendency is for the mean ﬁbre 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 ﬁbre 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, ﬁbre alignment, or orientation anisotropy as they called it, decreases signiﬁcantly in a small region near a solid boundary, suggesting the idea that solid boundaries have a randomizing eﬀect on ﬁbre orientation. In studies by H¨am¨al¨ainen & H¨am¨al¨ainen (2007), ﬁbre 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 justiﬁcation for this choice of boundary conditions were given. Studies of near wall ﬁbre orientation by Stover, & Cohen (1990) showed that ﬁbres may orient away from walls in laminar ﬂow, but only under certain conditions. They found that when ﬁbres with a low Jeﬀery orbit constant, that is highly aligned ﬁbres, came within a half ﬁbre length from the wall the ﬁbre remained close to, and aligned with the wall indeﬁnitely. But ﬁbres 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, ﬁbres might align in the wall direction or they might reorient away from the wall. It should be pointed out that 5 mm length ﬁbres were used in this study which is roughly the same length scale as the wall layer region. It is possible that for ﬁbre lengths comparable to the length scale of the wall layer, ﬁbres will align with the 61 Chapter 5. Validating Dr wall, where as smaller ﬁbres might not. This conjecture was not invetigated in this work. Other more subtle diﬀerences between the model predictions and the experimental observations, such as diﬀerences 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 ﬁbre orientations onto the xy-plane. It might be interesting to solve the full 3D ﬁbre 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 diﬃculty with this method is that it is quite computationally intensive, particularly near the solid boundaries where the ﬂow gradients become large in non-streamwise directions. 5.4 Eﬀect of 2 Way Coupling In this section we discuss the eﬀects of the two way coupling between the ﬂuid and ﬁbre phases. That is, we investigate the altered structure of the resulting ﬂow ﬁeld, 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 ﬂow ﬁeld for the diﬀerent concentrations under investigation. Shown also is a contour plot of the magnitude of the ﬂow ﬁeld for the pure ﬂuid, that is, water with no ﬁbres 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 eﬀect of the ﬁbre phase is to redirect mass near the walls towards the middle of the contraction, eﬀectively reducing the velocity magnitude along and around the centerline and creating something similar to a plugtype ﬂow condition. The predicted ﬂow ﬁeld with the ﬁbre phase is considerably diﬀerent than the ﬂow ﬁeld without ﬁbres, where the velocity is seen to be at a maximum along the centerline and decreases towards the walls, a phenomenon typical of Newtonian channel ﬂow. However when the ﬁbres 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 ﬂows of a yield stress ﬂuid, 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 ﬂow actually becomes slowest along the centerline. The eﬀect of increasing the ﬁbre concentration on this observation is a very subtle one, where the ﬂow becomes marginally slower as the concentration increases. This sort of ﬂow structure has been observed by other researchers when modelling ﬁbre suspension ﬂows. For example, Heath et al. (2007) experimentally studied the ﬂow 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 ﬂow 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 ﬂow of ﬁbre suspensions using non-Newtonian stress models have seen limited success. For example, Duﬀy (2003) reviewed many attempts at using non-Newtonian 64 Chapter 5. Validating Dr models such as shear thinning and thickening, Bingham plastic, and friction factorReynolds number method and showed that none of these models were suitable when describing ﬁbre suspension ﬂow. Duﬀy (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 conﬁrm the predictions made in this thesis with experimental measurements. 5.4.2 Orientation Distribution Here we look at the eﬀect 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 distribution 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 eﬀect of the two way coupling on the orientation distribution function is very small. The diﬀerences are summarized quantitatively in Table 5.1 below where we give the L2 norm of the diﬀerence between coupled and uncoupled solutions at the end of the contraction. Also given is the relative diﬀerence 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: Diﬀerence Between Coupled and Uncoupled Orientation Distributions nL3 ||ΨCoup − ΨN oCoup ||2 ||ΨCoup − ΨN oCoup ||2 /||ΨN oCoup ||2 8 0.3561 0.0135 16 0.3064 0.0119 24 0.1039 0.0045 66 Chapter 5. Validating Dr 7 6.5 6 6 5 5.5 Ψ Ψexit 4 5 exit 3 4.5 2 4 1 3.5 0 −2 Coupled Uncoupled −1.5 −1 −0.5 0 φ 0.5 1 1.5 3 −0.1 2 6 6 5 5.5 −0.05 0 φ 0.05 0.1 Coupled 4 5 Ψ Ψ exit exit3 Ψ 4.5 2 4 1 3.5 0 −2 −1.5 −1 −0.5 0 φ 0.5 1 1.5 3 −0.1 2 5 5 4.5 4.8 4 4.6 3.5 4.4 3 Ψ −0.05 0 φ 0.05 0.1 Coupled Uncoupled 4.2 exit exit2.5 4 2 3.8 1.5 3.6 1 3.4 0.5 0 −2 Uncoupled 3.2 −1.5 −1 −0.5 0 φ 0.5 1 1.5 2 3 −0.1 −0.05 0 φ 0.05 0.1 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 eﬀects of long range hydrodynamic ﬁbre-ﬁbre interactions on the orientation state of a semi-dilute, rigid ﬁbre suspension ﬂowing through a linear contracting channel under laminar conditions. The eﬀects of ﬁbreﬁbre 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 ﬁbre-ﬁbre interactions are assumed to be random in nature and eﬀectively result in a ﬂux which opposes the local mean orientation state of the ﬁbres thus mimicking diﬀusion-type process. In the ﬁrst part of this work, the diﬀusion coeﬃcient was measured along the central streamline of the contraction and the utility of two diﬀerent models for the rotary diﬀusion were compared. Each of the two models contains one unknown proportionality 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 distribution was solved and the solution ﬁt at the contraction exit for four diﬀerent ﬁbre 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 dimensions over the entire plane of the contraction using the measured values for the diﬀusion coeﬃcient from part one. A two way momentum coupling between the ﬂuid and ﬁbre phases was assumed and its eﬀect on both the suspension orientation state and structure of the ﬂow ﬁeld was investigated. 6.2 Diﬀusion Coeﬃcient The rotary diﬀusion coeﬃcient was measured along the contraction centerline as a function of ﬁbre concentration. In either case it was found to ﬁrst increase with increasing suspension concentration up to a maximum, and then decrease with concentrations above this point. This non-monotonic behavior was attributed to ﬁbre clumping or ﬂocculation, a mechanism not considered in the relationships for the rotary diﬀusion coeﬃcient. To be more speciﬁc, it is believed that at ﬁbre concentrations above a critical value, nL3 = 24 in this case, the assumption of purely hydrodynamic ﬁbre-ﬁbre interactions is no longer a valid one and that the eﬀects of mechanical contacts between ﬁbres play a more dominant role. The result is that ﬁbre mobility is greatly reduced, and that this reduction in mobility, combined with 69 Chapter 6. Summary and Conclusions the aligning eﬀects of the high ﬂow 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 interaction coeﬃcient, CI , is determined as a function of ﬁbre concentration. The results showed that CI increases linearly with concentration up to nL3 = 24 after which it suddenly decreases to a minimum value at nL3 = 32. This behavior indicates that CI is eﬀectively a measure of the diﬀusion coeﬃcient and it indicates a linearity in rotary diﬀusion 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 ﬁbre concentration. It was shown that λ decreases nonlinearly with increasing ﬁbre 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 diﬀusion coeﬃcient does not scale with nL3 , but perhaps with some other combination of n, L, and r. This ﬁnding is unfortunate as it would be highly desirable from a modelling standpoint to be able to determine a proportionality constant once for a particular ﬂow, then be able to reliably predict the eﬀect of varying other parameters without measuring a new proportionality constant. 70 Chapter 6. Summary and Conclusions 6.3 Comparison of Diﬀusion Models The utility of the two models for the rotary diﬀusion coeﬃcient were evaluated along the central streamline of the contraction by comparison with experiment. It was shown that both models give good predictions of ﬁbre orientation along the contraction once the diﬀusion coeﬃcient has been determined. It was shown quantitatively through computation of the χ2 goodness of ﬁt values that Koch’s model for Dr gives a marginal improvement in predicting ﬁbre 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 ﬂows as studied in this work. In general, it is concluded that the diﬀusion coeﬃcient must be determined empirically for a given set of suspension parameters and ﬂow ﬁeld, at least with the two models investigated in this work. This is a situation that continues to hinder the ﬁbre orientation problem. 6.4 2D Fibre Orientation The two dimensional ﬁbre 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 ﬁbre and ﬂuid phases. It was shown that ﬁ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 ﬁbre 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 ﬁbre orientation angle in a direction away from the wall. Experimental observations however showed no such behavior but instead showed that ﬁbres align with the wall indeﬁnitely. It is speculated that the discrepancy between the model predictions and experimental observations near the wall are based on the length of the ﬁbres used in the experiments. More speciﬁcally, it is believed that the ﬁbre is not aﬀected by the larger ﬂow gradients near the wall due to the fact that ﬁbre length is comparable to the length scale of the wall boundary layer. It is concluded that a more detailed study of wall-ﬁbre interactions is needed to determine sound relationships between ﬁbre length, aspect ratio and possibly concentration and near wall ﬁbre behavior. This is another situation that has in the past hindered, and will continue to hinder the ﬁbre orientation problem. 6.5 Two Way Coupling The eﬀect of the two way coupling on the orientation distribution and on the ﬂow ﬁeld has been investigated numerically. This was accomplished through an iterative procedure where the ﬂow ﬁeld was solved without the ﬁbre phase present and the orientation equations solved based on the initial ﬂow ﬁeld. The ﬁbre stress term was computed and used to recompute the ﬂow ﬁeld, this time with the ﬁbre phase present, after which the orientation equations were solved again with the updated ﬂow ﬁeld. The process was repeated until convergence was reached. The results showed that the eﬀect of the two way coupling on the 2D planar orientation distribution was 72 Chapter 6. Summary and Conclusions negligible. However the structure of the ﬂow ﬁeld was altered considerably. It was shown that the ﬂow ﬁeld changed from a parabolic-type proﬁle, typical of Newtonian channel ﬂows, to something which resembled a plug-type ﬂow. The change in ﬂow structure was attributed to a redirection of mass and momentum away from the walls towards the middle of the contraction, thus slowing the ﬂow away from the walls and creating a more uniform ﬂow distribution throughout the middle section of the channel. 73 Chapter 7 Future Work Several interesting ﬁndings have been made in this work. These include the observation that ﬁbres ﬂocculate at concentrations above some critical value, the shortcomings of the available models for the rotary diﬀusion coeﬃcient, the near wall behavior of ﬁbres, the boundary conditions on Ψ and the two-way interaction between the suspension orientation state and the resulting ﬂow ﬁeld. It is recommended that these issues all be further studied in greater detail. 7.1 Mechanical Fibre-Fibre Interactions In order to study the eﬀects of ﬂocculation 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 eﬀects of mechanical ﬁbre-ﬁbre 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 eﬀect of variable ﬁbre concentra- 74 Chapter 7. Future Work tions throughout the channel would also likely improve the predictions made in this work. 7.2 Rotary Diﬀusion Coeﬃcient As was shown in this work, neither of the two models for the rotary diﬀusion coeﬃcient were able to model diﬀerent 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 diﬀusion coeﬃcient could possibly eliminate this problem. Phan-Thien et al. (2002) addressed this issue in their work where they developed numerical correlations between the Jeﬀerey orbit constant and the concentration parameter c · r as opposed to nL3 , where c was the volume fraction of ﬁbres and r the ﬁbre aspect ratio. This may provide a starting point for such a study. 7.3 Near Wall Behavior of Fibres The behavior of ﬁbres near solid boundaries still remains a big question mark in the ﬁbre 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 ﬁbres 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 eﬀects of ﬁbre length, aspect ratio, suspension concentration and possibly ﬂow conditions on the near wall behavior of ﬁbres and ﬁbre suspensions. 7.4 Two way coupling The ﬁnal recommendation for future work would be to measure the ﬂow ﬁeld associated with the ﬂow of non-dilute suspensions in contraction type ﬂows 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 ﬂow, such as turbulent ﬂows for example, with improved conﬁdence. 76 Bibliography Advani, S.G., Tucker, C.L., 1987 The Use of Tensors to Describe and Predict Fibre Orientation in Short ﬁbre 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. ¨ ceri, S.I., Pipes, R.B. 1989 On the DeAltan, M.C., Advani, S.G., Gu 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 Proﬁle 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 CrossSection 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 SemiConcentrated Fibre Suspensions J. Rheol. 28(3), 207–227 Doi, M. & Edwards, S.F. 1984 The Theory of Polymer Dynamics Oxford University Press New York Duffy, G. G. 2003 The Signiﬁcance of Mechanistic-Based Models in Fibre Suspension Flow Nordic J. of Pulp Paper Sci. 18(1), 74–80. Folgar, F. 1983 Orientation behavior of ﬁbres in concentrated suspensions Ph.D. thesis, University of Illinois at Urbana-Champaign, 1983 Folgar, F. & Tucker, C.L. 1984 Fibre Orientation Distribution in Concentrated Suspensions: A Predictive Model J. Reinf. Plastics and Comp. 3 98–119 78 Bibliography ¨ma ¨la ¨inen, T. & Ha ¨ma ¨la ¨inen, J. 2007 Modelling of Fibre Orientation in the Ha Headbox Jet J. Pulp Paper Sci. 33(1), 49–53. Heath H.S., Olson, J.A., Buckley, K.R., Lapi, S., Ruth, T.J., & Martinez, D,M. 2007 Visualization of the Flow of a Fibre Suspension Through a Sudden Expansion Using PET AIChE Journal 53(2), 327–334. ¨ derberg, D. 2007 Modelling the Eﬀect of Shear Flow on Fibre Holm, R. & So Orientation Anisotropy in a Planar Contraction Overview of forming literature, 1990-2000 Trans. 12th Fundamental Research Symposium, Oxford 1 431–558 ¨ , M., Dahlkild, A., Krochak, P., Olson, J. 2007 Shear Inﬂuence Hyensjo on Fibre Orientation Nordic Pulp and Paper Research Journal 22, 376–382. ¨ m, A., Thurigaswamy, Martinez, D.M., Buckley, K., Jivan, S., Lindstro R., Ruth, T.J. & Kerekes, R.J. 2001 Characterizing the Mobility of Fibres During Sedimentation Pulp Pap. Fund. Res. Symp. (Oxford) 1 225–254 ¨ derberg, D. 2001 Overview of Forming Literature, 1990-2000 Norman, B. & So Trans. 12th Fundamental Research Symposium, Oxford 1 431–558 Jackson, W.C., Advani, S.G. & Tucker, C.L. 1985 Predicting the Orientation 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 Eﬀects on the Motion of Long Slender Bodies J.Fluid Mech. 209 435–462 79 Bibliography Koch, D.L. 1995 A Model for Orientational Diﬀusion in ﬁbre Suspensions Phys. Fluids 7(8) 2086–2088 Leal, L.G. & Hinch, E.J. 1971 The Eﬀect 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 Function 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 ¨ma ¨la ¨inen, J. P. 2004 Modeling a Olson, J.A., Frigaard, I., Chan, C., & Ha 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 Stiﬀ 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 Eﬀect of Hydrodynamic 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 Orientation in Suspensions Subject to Planar Extensional Flows Phys. Fluids 7(8) 1811–1817 Ranganathan, S. & Advani, S.G. 1991 Fibre-Fibre Interactions in Homogeneous 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 Headbox 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 Geometries: Flow/ﬁbre Coupling Canadian Journal of Chemical Engineering 80 1093– 1106 Yasuda, K., Kyuto, T. & Mori, N. 2004 An Experimental Study of FlowInduced Fibre Orientation and Concentration Distributions in a Concentrated Suspension 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 ﬂow loop is constructed entirely of plexiglass with fully smooth walls, 5 mm thick. Upstream of the experimental section is a short rectangular channel section where the ﬂow is allowed to become fully developed. A dimensioned 83 Appendix A. Experimental System sketch of the experimental section is shown below, see ﬁgure A.2. Figure A.2: A dimensioned sketch of the experimental section. The suspension is continously recirculated with a Masterﬂex B/T variable speed peristaltic pump (www.coleparmer.com). The speciﬁcations 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 ﬁbre phase suspended in a glycerine ﬂuid phase. The physical data for the suspending ﬂuid, 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 ﬁbres used in the experiments is given below: • 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 ﬁbre was exactly two pixels, and (c) the relative change in optical path due to the diﬀerences of index of refraction as light travels from the suspension to the camera lens is negligible. The speciﬁcations 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 ﬁbre 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 presented. This includes the programs used to track ﬁbres within the experimental images of the ﬂowing suspensions and the programs used to analyze the raw data generated by the ﬁbre tracking algorithms. Also included in this section is the code used to perform the Eularian numerical predictions of ﬁbre orientation. That is, the numerical solvers used to solve Equations 4.4 and 5.1 are given B.1 B.1.1 Fibre Tracking and Analysis Fibre Tracking The programs used to track ﬁbres 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 ﬁbres 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 ﬁbre in the black and white image. This function also ﬁlters images for unwanted debris and crossing ﬁbres. The values of the ﬁbre end points is determined by the function FindEndPoints.m called within idFibresE cor.m. These values are then saved in a seperate ﬁle. 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 %End vw loop....vw = view = [top,bottom] % 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 ﬁbre tracking algorithm is converted to position and orientation data with the code presented below. Included in this is the routine used to measure the ﬁbre 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 ﬁbre data from each image, that is the data representing the centre of area and two end points of each observed ﬁbre as measured with the code from §B.1.1, and computes the orientation of each ﬁbre. These orientations are binned into the correct spatial cell based on the centre of area of the ﬁbre 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 ye_b kt = kb = = y0 = ye ye_t ye_b + + + + 0.3*(y1-y0); 0.3*(y4-ye); (y0_t-ye_t)/(x2-x3)*x; (y0_b-ye_b)/(x1-x4)*x; %%% Top perimeter of Window %%% 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 ﬁbre 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 ﬁbre 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 solvePsiStr.m respecively. The code shown here for solving Equation 4.4 is for the case when Equation 4.2 is used to deﬁne the diﬀusion coeﬃcient. 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 ﬁbre stress term and it’s gradient is given below. This routine ﬁrst 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 ﬂuid momentum equations. The ﬁbre 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 ﬂow ﬁeld. The resulting ﬁle containing this data is later loaded into the Fluent environment to solve the ﬂuid ﬂow 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 |
File Format | 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 |
Graduation Date | 2008-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | 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:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0066399/manifest