MODELING OF FIBER MOTION IN PULP AND PAPER EQUIPMENT By Suqin Dong B. Sc., Zhengzhou University, China, 1989 M. Sc., Chinese Academy of Engineering Physics, China, 1992 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS OF THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MECHANICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA August 2002 © Suqin Dong, 2002 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mechanical Engineering The University of British Columbia 2324 Main Mall Vancouver, BC Canada V6T 1Z4 Date: August 2002 Abstract This thesis describes a computational tool developed to simulate the fiber motion in pulp and paper equipment with particular emphasis on pulp screens. It is focused on modeling the fiber concentration upstream of a screen slot and examining the effects that different slot shapes have on the passage of fibers through the narrow apertures. This provides insight into the screen performance and has been used in this research to predict the characteristics of the screen. This simulation tool includes three coupled models: the flow model to predict the flow field in the equipment, the fiber model to trace the fiber trajectory in the equipment and the wall model to deal with the case when a fiber touches the equipment wall. The fluid field in the screen was examined using two numerical methods: the Reynolds Averaged Navier-Stokes (RANS) method for the calculation of flow through various slot shapes and Large Eddy Simulation (LES) for the fiber concentration calculation upstream of the screen slot. A three-dimensional flexible fiber model was employed to track the fiber trajectory in the screen. A very general wall model was developed to deal with the fiber-wall interaction. This wall model could be used for the case where a fiber touches a wall of any geometry, while specifying the wall geometry only through the fluid flow calculation routine. The simulated fiber concentration distribution increases linearly near the wall and becomes approximately constant farther from the wall, in reasonable agreement with experimental observations reported elsewhere. Slot shape has the most important influence on fiber passage ratio. For the same kind of fibers, same flow conditions and same overall slot width, a slope-slope slot provides the best passage for the fiber among all the slot shapes tested. For these contour slots, the contour height is also very critical, and a contour depth of 0.5mm seems best for both the step-step contour and slope-slope contour for fiber lengths of 1 to 3 mm. For the contour slots, the upstream side of contour plays no significant role in the fiber behavior; only the downstream side is important. ii Tables of contents Abstract ii List of Tables vi List of Figures vii Acknowledgements xiv Chapter 1 Introduction 1 Chapter 2 Literature review 5 2.1 Screening Concepts 5 2.2 Screening Performance 5 2.3 Fiber Passage Ratio 8 2.4 Screen Plate Design 12 2.5 Particle Dispersion Two-Phase Flow 14 2.6 Pulp Fiber Models 16 Chapter 3 Fiber-dispersed turbulent flow simulation 18 3.1 Introduction 18 3.2 Fiber-dispersion turbulent flow simulation 18 3.2.1 Random Walk method 18 3.2.2 Large Eddy Simulation (LES) 23 Chapter 4 Fiber model, wall model and algorithm 25 4.1 Introduction 25 4.2 Fiber model 25 4.2.1 Configuration and position of fiber 26 4.2.2 Transformation matrix from the fiber body-fixed frame to the fixed reference frame 27 4.2.3 Dynamics 28 4.3 Wall model 32 4.3.1 The condition for the fiber to touch the wall 33 4.3.2 The force equation for the spheroid when touching the wall 34 4.4 Calculation procedure 37 4.5 Algorithms 40 ii i iv 4.5.1 Euler method 41 4.5.2 A Variable order Runge-Kutta Method (VRK) 41 4.6 Validation example 42 4.6.1 Fiber moving in a simple shear flow 42 4.6.2 A more general example 44 4.7 Discussion of the resistance force 47 48 4.8 Fiber flexibility Chapter 5 Numerical simulation of fiber concentration 49 5.1 Introduction 49 5.2 Flow simulation 50 5.2.1 Reynolds-Averaged Navier-Stokes (RANS) method 51 5.2.2 Large Eddy Simulation (LES) 56 5.3 Fiber motion and concentration in turbulent mean flow 60 5.4 Fiber motion in turbulent flow with diffusion of fibers 63 5.4.1 Concentration from Random Walk 64 5.4.2 Concentration from LES 65 5.5 Discussion 67 5.5.1 Effects of some factors on the fiber concentration 67 5.5.2 The orientation of fiber 70 5.6 Summary 73 Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 75 6.1 Introduction 75 6.2 Flow field simulation 75 6.2.1 Slot geometry and boundary conditions 76 6.2.2 Grid independence 78 6.2.3 Convergence 83 6.2.4 Results 84 6.3 Fiber motion in different slot shapes 100 6.3.1 Trajectories of 1mm rigid fiber (L/W = 2) 102 6.3.2 Trajectories of 3mm rigid fiber (L/W = 6) I l l V 6.3.3 Trajectories of 1mm flexible fiber (LIW = 2) 115 6.3.4 Trajectories of 3mm flexible fiber (L/W = 6 ) 120 6.4 Stapling 123 6.4.1 Stapling in the smooth slot 124 6.4.2 Stapling in the Cl slot 125 6.4.3 Stapling in the C2 slot 126 6.4.4 Stapling in L2 slot 128 6.5 Fiber passage ratio in different slot shapes 128 6.5.1 Passage ratio for 1mm rigid fiber 129 6.5.2 Passage ratio for 3mm rigid fiber 133 6.5.3 Passage ratios for flexible fibers 135 6.6 Summary 138 Chapter 7 Other applications 140 7.1 Introduction 140 7.2 Fiber orientation in the headbox 140 7.3 Fiber separation in hydrocyclone 142 Chapter 8 Summary and conclusions 145 8.1 Three models 145 8.2 Fiber concentration upstream of the slot 146 8.3 The flow and fiber passage in slots with different shapes 146 Chapter 9 Recommendations for future work 148 Nomenclature 149 Bibliography 155 List of Tables TABLE 4.1 THE RESISTANCE FUNCTIONS FOR A SINGLE SPHEROID (KIM AND KARRILA 1991) 30 TABLE 4.2 COMPARISON OF C P U TIME FOR EULER AND V R K METHODS 47 TABLE 5.1 DIFFERENT CASES FOR RANS CALCULATION 53 TABLE 5.2 RATIO OF FIBER LENGTH TO GRID SIZE IN X DIRECTION FOR RANS 53 TABLE 5.3 DIFFERENT CASES FOR LARGE EDDY SIMULATION CALCULATION 57 TABLE 6.1 GRIDS USED FOR TESTING GRID INDEPENDENCE FOR SMOOTH SLOT 79 TABLE 6.2 SYMBOLS USED IN ALL THE "FATE PLOTS" 101 TABLE 6.3 PASSAGE RATIO OF 1MM RIGID FIBER IN SLOTS WITH DIFFERENT SHAPES 132 TABLE 7.1 PROPERTIES OF NYLON FIBERS USED IN THE EXPERIMENTS (Ho 2001) 144 vi List of Figures FIGURE 1.1 A TYPICAL PRESSURE SCREEN 2 FIGURE 1.2 HOLES OR SLOTS IN A SCREEN PLATE (PROVIDED BY PAPRICAN) 3 FIGURE 2.1 SCREEN PERFORMANCE CURVES (2.1), (2.2) FOR DIFFERENT a AND /? 7 FIGURE 2.2 SCREEN PERFORMANCE CURVES (2.3) ; 8 FIGURE 2.3 DEFINITION OF PASSAGE RATIO 9 FIGURE 2.4 PASSAGE RATIO FOR 1 MM, 3MM RIGID FIBERS WITH SLOT WIDTH 0.5MM (KUMAR, 1991) 10 FIGURE 2.5 SOME COMMERCIAL CONTOUR DESIGNS (GOODING 1996) 13 FIGURE 3.1 8TH ORDER POLYNOMIALS CURVE-FITTING COMTE-BELLOT'S EXPERIMENT (1965).... 21 FIGURE 4.1 A FIBER CONSISTING OF SPHEROIDS CONNECTED THROUGH BALL AND SOCKET JOINTS (ROSS AND KLINGENBERG 1997) 27 FIGURE 4.2 CONNECTIVITY VECTOR AND TWO REFERENCE FRAMES 27 FIGURE 4.3 FREE-BODY DIAGRAM FOR SPHEROID / (Ross AND KLINGENBERG 1997) 28 FIGURE 4.4 SCHEMATIC DIAGRAM OF THE FIBER TOUCHING THE WALL 34 FIGURE 4.5 A N EXAMPLE OF THE FIBER TOUCHING THE WALL 37 FIGURE 4.6 FIBER'S INITIAL POSITION 39 FIGURE 4.7 FIBER ROTATES IN X Y PLANE 42 FIGURE 4.8 FIBER TRANSLATES AND ROTATES IN X Y PLANE 42 FIGURE 4.9 FIBER TRANSLATES AND ROTATES IN 3D 43 FIGURE 4.10 FLEXIBLE FIBER TRANSLATES 43 FIGURE 4.11 COMPARISON BETWEEN FLEXIBLE MODEL AND JEFFERY'S PREDICTION 44 FIGURE 4.12 GEOMETRY OF THE SYMMETRICAL CONVERGENT SECTION OF A HEADBOX 45 FIGURE 4.13 COMPARISON OF ANALYTICAL SOLUTION AND VRK SOLUTION 46 FIGURE 4.14 COMPARISON OF THE RELATIVE ERROR OF VRK AND THE FIRST EULER METHOD 47 FIGURE 5.1 BOUNDARY CONDITION FOR RANS OF THE CHANNEL FLOW 51 FIGURE 5.2 GRIDS OF RANCASE 1: NON-UNIFORM IN Z DIRECTION AND UNIFORM IN x DIRECTION. 54 FIGURE 5.3 COMPARISON OF MEAN VELOCITIES U FROM RANS FOR DIFFERENT CASES 55 FIGURE 5.4 u VELOCITY ACROSS THE CHANNEL FOR DIFFERENT CASES 56 vii Vlll FIGURE 5.5 COMPARISON OF U+ VS. Z + FROM L E S OF DIFFERENT CASES WITH SPALDING'S THEORETICAL SOLUTION AND THREE EXPERIMENTAL DATA : LIU ET AL. (Re r =184), COMTE-BELLOT ( ReT = 2340) AND WEI & WlLLMARTH (Re T =1655) 57 FIGURE 5.6 u'lut CALCULATED BY L E S OF THREE CASES, COMPARED WITH EXPERIMENTS OF: LIU ET AL. (Re R=184), COMTE-BELLOT (Re R=2340) AND WEI & WlLLMARTH (Re r =1655) 59 FIGURE 5.7 v'/uT CALCULATED BY L E S OF THREE CASES, COMPARED WITH EXPERIMENTS OF COMTE-BELLOT (ReT = 2340) (OTHER EXPERIMENTS ARE NOT AVAILABLE) 59 FIGURE 5.8 w 7 « r C A L C U L A T E D BY L E S OF THREE CASES, COMPARED WITH EXPERIMENTS OF: LIU ET AL. (Re R=184), COMTE-BELLOT (Re R =2340) AND WEI & WlLLMARTH (Re r =1655) 60 FIGURE 5.9 THE MEAN VELOCITY PROFILE IN xz PLANE 61 FIGURE 5.10 CENTERS OF FIBERS AT X Z PLANE WITH 62 FIGURE 5.11 CENTERS OF FIBERS AT THREE- DIFFERENT INITIAL HEIGHT. DIMENSIONAL SPACE WITH DIFFERENT INITIAL 62 FIGURE 5.12 CENTERS OF FIBERS AT THREE-DIMENSIONAL 62 FIGURE 5.13 PROBABILITY OF 1MM WITH DIFFERENT INITIAL HEIGHT AND DIFFERENT ANGLE 9. FIBER 2 D AND 3 D TURBULENT MAIN 62 FIGURE 5.14 CONCENTRATION FOR 1MM, 2MM AND 3MM FIBERS VS. CHANNEL HEIGHT, USING RANDOM WALK 64 FIGURE 5.15 CONCENTRATION FOR 1MM, 2MM , 3MM FIBERS VS. CHANNEL HEIGHT SCALED BY FIBER LENGTH, USING RANDOM WALK 65 FIGURE 5.16 CONCENTRATION FOR 1MM, 2MM , 3MM FIBER VS. CHANNEL HEIGHT SCALED BY FIBER LENGTH USING L E S 66 FIGURE 5.17 CONCENTRATION FOR 2MM FIBER WITH DIFFERENT NUMBER OF SPHEROIDS 67 FIGURE 5.18 CONCENTRATION OF 8,000,16,000 AND 24,000 3MM RIGID FIBERS 68 FIGURE 5.19 CONCENTRATION AT DIFFERENT STATIONS 69 FIGURE 5.20 CONCENTRATION FOR 3MM RIGID AND FLEXIBLE FIBER 70 FIGURE 5.21 PROBABILITY DENSITY CONTOURS FOR HEIGHT z / L AND ORIENTATION FOR 1MM FIBERS IN XZ PLANE 71 ix FIGURE 5.22 PROBABILITY DENSITY CONTOURS FOR HEIGHT Z/L AND ORIENTATION FOR 2MM FIBERS IN XZ PLANE 71 FIGURE 5.23 PROBABILITY DENSITY CONTOURS FOR HEIGHT z/L AND ORIENTATION FOR 3MM FIBERS IN XZ PLANE 72 FIGURE 5.24 THE COMPARISON OF AVERAGE ORIENTATION IN xz PLANE FOR 1MM, 2MM AND 3MM FIBERS 72 FIGURE 6.1 GEOMETRY OF A SMOOTH SLOT 77 FIGURE 6.2 SLOT GEOMETRY AND DIMENSIONS OF STEP-STEP AND SLOPE-SLOPE CONTOURS (ALL DIMENSIONS ARE IN MM) 78 FIGURE 6.3 Two POSITIONS FOR COMPARISON OF VELOCITIES WITH DIFFERENT MESHES 79 FIGURE 6.4 u VELOCITY COMPARISON ABOVE THE SLOT BETWEEN DIFFERENT GRIDS 80 FIGURE 6.5 w VELOCITY COMPARISON ABOVE THE SLOT BETWEEN DIFFERENT GRIDS 81 FIGURE 6.6 u VELOCITY COMPARISON ALONG THE SLOT CENTER BETWEEN DIFFERENT GRIDS 81 FIGURE 6.7 w VELOCITY COMPARISON ALONG THE SLOT CENTER BETWEEN DIFFERENT GRIDS 82 FIGURE 6.8 GRIDS FOR SMOOTH SLOT 82 FIGURE 6.9 RESIDUALS OF SIX VARIABLES FOR SEGMENT 1 OF SMOOTH SLOT 83 FIGURE 6.10 RESIDUALS OF SIX VARIABLES FOR SEGMENT 2 OF SMOOTH SLOT 84 FIGURE 6.11 STREAMLINE IN THE SMOOTH SLOT FOR DIFFERENT SLOT VELOCITIES 85 FIGURE 6.12 STREAMLINE IN Cl SLOT FOR DIFFERENT SLOT VELOCITIES 86 FIGURE 6.13 STREAMLINE IN L2 SLOT FOR DIFFERENT SLOT VELOCITIES 87 FIGURE 6.14 FLOW FIELD FOR STEP-STEP CONTOUR WITH DIFFERENT CONTOUR HEIGHT 88 FIGURE 6.15 FLOW FIELD FOR SLOPE-SLOPE CONTOUR WITH DIFFERENT CONTOUR HEIGHT 89 FIGURE 6.16 GEOMETRY FOR CONTOUR WITH SAME HEIGHT 0.5MM AND DIFFERENT SHAPE 90 FIGURE 6.17 FLOW FIELD FOR CONTOURS WITH SAME HEIGHT 0.5MM AND DIFFERENT SHAPES 91 FIGURE 6.18 Two CONTROL VOLUMES USED BY GOODING (1996) 92 FIGURE 6.19 PRESSURE LOSS FOR THE SMOOTH SLOT COMPARED WITH GOODING'S RESULTS 95 FIGURE 6.20 PRESSURE DROP COEFFICIENT FOR SMOOTH SLOT COMPARED WITH GOODING'S RESULTS 95 FIGURE 6.21 PRESSURE DROP FOR SMOOTH SLOT, Cl SLOT AND C2 SLOT 96 FIGURE 6.22 PRESSURE DROP COEFFICIENT FOR SMOOTH SLOT, Cl SLOT AND C2 SLOT 96 FIGURE 6.23 TURBULENCE INTENSITY IN SMOOTH SLOT FOR LOW SLOT VELOCITY (Vs =2.4 M/S).97 X FIGURE 6.24 TURBULENCE INTENSITY IN SMOOTH SLOT FOR HIGH SLOT VELOCITY (VS =5.4 M/S). 98 FIGURE 6.25 CONTOUR OF TURBULENCE INTENSITY FOR SMOOTH, CO, Cl, C2 SLOTS 99 FIGURE 6.26 CONTOUR OF TURBULENCE INTENSITY FOR CISMOOTH, C I S L O P E , C I T W O SLOPE AND L l SLOTS 100 FIGURE 6.27 SMOOTH SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION 6 ON PASSAGE OF 1MM RIGID FIBER FOR DIFFERENT SLOT VELOCITIES (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL" ": STAPLED) 104 FIGURE 6.28 SMOOTH SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION <p ON PASSAGE OF LMM RIGID FIBER FOR DIFFERENT SLOT VELOCITIES (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL " ": STAPLED) 105 FIGURE 6.29 Cl SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION 0 ON PASSAGE OF 1MM RIGID FIBER FOR DIFFERENT SLOT VELOCITIES (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL" ": STAPLED) 106 FIGURE 6.30 Cl SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION <p ON PASSAGE OF 1MM RIGID FIBER FOR DIFFERENT SLOT VELOCITIES (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL " ": STAPLED) 107 FIGURE 6.31 CO, Cl, C2 SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION 6, <p ON PASSAGE OF 1MM RIGID FIBER FOR SLOT VELOCITY VS = 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL " ": STAPLED) 108 FIGURE 6.32 L l , L2 SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION 6, cj) ON PASSAGE OF 1MM RIGID FIBER FOR SLOT VELOCITY VS = 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL" ": STAPLED) 109 FIGURE 6.33 C ISMOOTH, C I S L O P E , C1_TWO_SLOPE SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION 6, tp ON PASSAGE OF 1MM RIGID FIBER FOR SLOT VELOCITY VS = 5.4 xi M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL " ": STAPLED) 110 FIGURE 6.34 SMOOTH SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION 9 , <j> ON PASSAGE OF 3MM RIGID FIBER FOR SLOT VELOCITIES 2.4 M/S, 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL " ": STAPLED) 112 FIGURE 6.35 Cl SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION S, <f> ON PASSAGE OF 3MM RIGID FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL" ": STAPLED) 113 FIGURE6.36 C2 SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATIONS, <p ON PASSAGE OF 3MM RIGID FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (SYMBOLS SAME AS FIGURE 6.35) 114 FIGURE 6.37 SLOPE-SLOPE L2 SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATIONS , </> ON PASSAGE OF 3MM RIGID FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL " ": STAPLED) 115 FIGURE 6.38 SMOOTH SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION 9, <p ON PASSAGE OF 1MM FLEXIBLE FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL " ": STAPLED) 116 FIGURE6.39 Cl SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATIONS, <j> ON PASSAGE OF LMM FLEXIBLE FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL " " : STAPLED) 117 FIGURE6.40 C2 SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATIONS, $ ON PASSAGE OF LMM FLEXIBLE FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL " ": STAPLED) 118 FIGURE6.41 L2 SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATIONS, <j> ON PASSAGE OF 1MM FLEXIBLE FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (GREEN SYMBOL "+": X l l PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL" ": STAPLED) 119 FIGURE 6.42 SMOOTH SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATIONS , </> ON PASSAGE OF 3MM FLEXIBLE FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL " ": STAPLED) 120 FIGURE6.43 Cl SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATIONS, <j> ON PASSAGE OF 3MM FLEXIBLE FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "O": GOES AWAY FROM THE CHANNEL; RED SYMBOL" ": STAPLED) 121 FIGURE 6.44 C2 SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATION e, </> ON PASSAGE OF 3MM FLEXIBLE FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (SYMBOLS SAME AS FIGURE 6.43) 122 FIGURE6.45 L2 SLOT: THE EFFECTS OF INITIAL HEIGHT AND ORIENTATIONS, ON PASSAGE OF 3MM FLEXIBLE FIBER FOR SLOT VELOCITIES 2.4 M/S AND 5.4 M/S (GREEN SYMBOL "+": PASSES THROUGH THE SLOT; BLUE SYMBOL "o " : GOES AWAY FROM THE CHANNEL; RED SYMBOL" ": STAPLED) 123 FIGURE 6.46 VERTICAL STAPLE AND HORIZONTAL STAPLE 124 FIGURE 6.47 FIBER STAPLING IN SMOOTH SLOT 125 FIGURE 6.48 FIBER STAPLING IN Cl SLOT 126 FIGURE 6.49 FIBER STAPLING IN C2 SLOT 127 FIGURE 6.50 FIBER STAPLING IN L2 SLOT 128 FIGURE 6.51 PASSAGE RATIO OF 1MM RIGID FIBERS IN SMOOTH SLOT 130 FIGURE 6.52 PASSAGE RATIO OF 1MM RIGID FIBERS IN C l , C2 SLOT 131 FIGURE 6.53 PASSAGE RATIO OF 1MM RIGID FIBERS IN DIFFERENT SLOT SHAPES 131 FIGURE 6.54 PASSAGE RATIO OF 3MM RIGID FIBERS IN SMOOTH SLOT 134 FIGURE 6.55 PASSAGE RATIO OF 3MM RIGID FIBERS IN C l , C2 SLOT 134 FIGURE 6.56 PASSAGE RATIO OF 3MM RIGID FIBERS IN DIFFERENT SLOT SHAPES 135 FIGURE 6.57 PASSAGE RATIO OF 1MM FLEXIBLE FIBERS IN SMOOTH SLOT 136 FIGURE 6.58 PASSAGE RATIO OF 1MM FLEXIBLE FIBERS IN C l , C2 SLOT 137 FIGURE 6.59 PASSAGE RATIO OF 3MM RIGID AND FLEXIBLE FIBERS IN SMOOTH SLOT 137 FIGURE 6.60 PASSAGE RATIO OF 3MM RIGID AND FLEXIBLE FIBERS IN C l , C2 SLOT 138 FIGURE 7.1 GRIDS AND THE MEAN FLOW FIELD IN A SYMMETRY HEADBOX 141 FIGURE 7.2 ORIENTATION IN TWO PLANES NEAR THE HEADBOX EXIT 141 FIGURE 7.3 DIAGRAM OF A HYDROCYCLONE (WANG 2002) 142 FIGURE 7.4 BAUER'S HYDROCYCLONE: COMPARISON OF MEAN COARSENESS IN INLET AND OUTLET STREAMS (WANG 2002) 144 Acknowledgements I would like to thank all those who have given me support to make this work possible. First, I would like to thank my research supervisors, Dr. Martha Salcudean and Dr. Ian Gartshore for their guidance, kind support, invaluable suggestions and encouragement. They have made my Ph.D. experience a wonderful period of my life, which I will cherish forever. I would also like to thank my research committee members, Dr. Carl F. Ollivier-Gooch, Dr. James A. Olson and Dr. Brian Wetton for their advice and helpful suggestions. Many thanks to Dr. Eric Bibeau and Dr. Paul Nowak for their useful discussions and help, and to my colleagues, Zheqiong Wang, Mohammad R. Shariati and Mike Georgallis, for their friendship and help. Special thanks to Dr. Matthew W. Choptuik for his generosity and technical support in the use of the PIII/Linux computer cluster in the Physics and Astronomy Department at the University of British Columbia. Thanks also to Alan Steeves for his help and technical support in using the facilities in the Mechanical Engineering Department at the University of British Columbia. I also wish to acknowledge the financial assistance from UBC (University Graduate Fellowship), Forest Research of BC and the Natural Sciences and Engineering Research Council. Finally, I would like to thank my parents, my husband Xiaosi and my daughter Manning for their love, understanding, unconditional support and help throughout this work. xiv Chapter 1 Introduction A pulp screen is a piece of equipment in the pulp and paper industry used either to remove contaminants from the pulp suspension or to separate fibers having different properties. Contaminants including fiber bundles, bark and plastic specks are introduced when fibers are separated from the wood by mechanical or chemical pulping processes. Contaminants significantly affect the strength and smoothness of the paper and must be removed before the final paper is produced. Therefore pulp screens are very important in the pulping and papermaking processes. With the increasing demand for more uniform and higher quality paper products, pulp screens are more and more used for "fractionation" which separates fibers according to their properties such as length, coarseness and flexibility. The pressure screen is the most common type of pulp screen and is the subject of this study. Figure 1.1 is a typical pressure screen. The feed (pulp suspension) enters tangentially at the top of the screen and travels through the gap between the screen plate and the rotor. "Accept" pulp passes through the screen plate to the accept outlet and the "Reject" pulp goes into the reject outlet. The rock trap stops heavy debris at the beginning from destroying the screen plate and rotor. Dilution water is introduced in some pressure screens to avoid thickening problems. The most important parts of the screen are the screen plate and rotor. The screen plate has holes or slots (Figure 1.2) which let the accept pulp go through. The holes or slots may have different shapes as discussed in a later section of this thesis. The rotor introduces pressure pulses to help prevent the screen plate from plugging and it also tangentially accelerates the pulp at the surface of the screen plate. 1 Chapter 1 Introduction 2 Figure 1.1 A typical pressure screen. Contaminant removal efficiency and fiber separation efficiency are the main criteria o f screen performance, with the fiber passage ratio determining the removal (separation) efficiency. The relationship between the fiber passage ratio and efficiency w i l l be discussed in chapter 2. Objective of the thesis: The overall objective o f this thesis is to develop a practical numerical tool which can be used to study fiber motion in pulp and paper equipment, and especially to examine the upstream fiber concentration and the effects o f the slot shapes on fiber passage. Chapter 1 Introduction 3 Figure 1.2 Holes or slots in a screen plate (provided by PAPRICAN) In order to achieve the objective, several steps are necessary: • The flow field in the equipment must be modeled; • The fibers that move in the flow field must be modeled; • Because the fibers often touch the equipment wall, a general wall model needs to be developed. "General" means that the wall model can be used for equipment of any kind of geometry without specifying the wall for each geometry, other than through the flow field calculation; • The flow, the fiber and the wall models need to be coupled together to simulate the complete motion of the fiber in the equipment. Following this introductory chapter, Chapter 2 is the literature review, covering several aspects of the screen: the general definition of the screen, screen performance, fiber passage ratio and screen plate design. Particle dispersion in turbulent flow and fiber model are also reviewed in chapter 2. Chapter 3 describes the two flow simulation methods used in this thesis for fiber dispersion in turbulent flow: Random Walk and Large Eddy Simulation (LES). Chapter 4 presents the fiber model, Chapter 1 Introduction 4 wall model and the corresponding algorithms. Chapter 5 describes the numerical simulation of fiber concentration upstream of the slot. In Chapter 6, a complete study of the effect of the screen plate shape on the flow field and passage ratio is described. Chapter 7 introduces other applications of the current numerical tool. Chapter 8 presents the summary and conclusions. Recommendations for future work are given in Chapter 9. Chapter 2 Literature review 2.1 Screening Concepts Screening is widely used to separate particles. There are two kinds of screening: solid-solid screening, which separates one type of solid particles from another, and solid-liquid screening, which separates solid particles from a suspending liquid. General screening principles are described in Tappi's Technical Information Sheets (1989). Depending on the dimension of the particles and the screen apertures, screening can be one of two types: barrier screening and probability screening. For barrier screening, only the particles whose dimensions are smaller than the apertures can go through. For probability screening, a particle which has one dimension larger than the aperture size can still go through the aperture provided it has the right approach conditions. Pressure screens (Figure 1.1) used in the pulp and paper industry are solid-solid screens. They are widely used to remove contaminants from the pulp and separate fibers of different types. The latter is called fiber fractionation. The simplest pressure screens have one inlet (feed) and two outlets (accept and reject). Possible contaminants have a broad range of sizes, whereas wood fibers are usually l-3mm long and 20-45 jum wide. The screen aperture is about 0.1-0.5 mm for a slot width and 0.8-2.5 mm for a hole diameter. Contaminant removal is accomplished by both barrier screening and probability screening and fiber fractionation is probability screening. 2.2 Screening Performance The contaminant removal efficiency (E) and reject rate (R) are the two most common parameters which describe screening performance. E is the ratio of the mass of contaminants in the reject flow to the mass of contaminants in the feed flow. R is the overall fraction of the feed pulp 5 Chapter 2 Literature review 6 that leaves in the reject stream (Tappi's Technical Information Sheets 1989). The relationship between these two parameters is used to judge the screen performance. Good performance is obtained with high contaminant removal efficiency at a low reject rate. An E-R curve is obtained experimentally in industry. There are two models commonly used to relate E and R : E = (2.1) 1-a+aR or E = R" (2.2) where a, ft are constant screen indices, found experimentally. A higher value of a and lower value of p denote better screen performance as shown in Figure 2.1. Equation (2.1) was first presented by Nelson (1981) with no derivation or description of assumptions. Gooding (1986) derived (2.1) for P solid-solid screening with a = 1 — - , where Pc Pp are the passage ratios of contaminants and pulp fibers as defined in section 2.3. He also referred to equation (2.1) as a "mixing flow" model under P the assumption of constant — and perfect mixing in the screen annular zone, which means that the concentration of pulp fibers in the reject stream is the same as that in the feed stream. Equation (2.2) is based on Kubat and Steenberg's flat screen analysis (1955) and was derived by Gooding (1986) P for a pressure screen with P = — and was called a "plug flow" model. The assumptions for (2.2) are that there is no mixing in the screen annular zone but that the flow in the radial direction is perfectly mixed. The importance of Gooding's derivations of equation (2.1) and (2.2) is that they directly relate the screen performance to the passage ratio. Equation (2.2) can also be expressed to describe both barrier screening and probability screening (Olson 2000): E = RP' r j _ ^ f c b a r C f < ) c Chapter 2 Literature review 7 where C / c is the total contaminant concentration upstream of the slot, Cf.c.prob is the contaminant concentration following probability screening and Cf ,c.har is the contaminant concentration following barrier screening. Figure 2.2 compares the efficiency of pure barrier screening, pure probability screening and both of these together. It appears that some percentage of barrier screening would provide the highest efficiency. Decreasing aperture dimensions can enhance barrier screening and attain high efficiency. On the other hand, smaller aperture dimensions can definitely decrease the capacity and passage. The choice of a screen plate is therefore a trade off between the efficiency and the capacity. The challenge of the current screen plate design is to keep the aperture size small and at the same time to increase the capacity. Recently designed contour slots make this possible. Figure 2.1 Screen performance curves (2.1), (2.2) for different a and /? . Chapter 2 Literature review 8 M a s s R e je ct ra tio R Figure 2.2 Screen performance curves (2.3). 2.3 Fiber Passage Ratio The "Passage ratio" was first defined and rated as an important performance parameter by Gooding (1986), and Gooding and Kerekes (1989) to relate the contaminant removal efficiency and reject rate, as explained in section 2.2. It was defined as the ratio o f the concentration o f pulp libers or contaminants in the flow through the aperture (C a ) and the corresponding concentration upstream o f the aperture ( C u ) . Gooding and Kerekes also experimentally studied the fiber passage through a single aperture using a single slot screen. A s illustrated in Figure 2.3, the average upstream flow velocity Vu comes from the rotation o f the rotor and is tangential to the screen plate, and Vs is the average flow velocity through the slot. They filmed the trajectories o f individual nylon fibers and deduced that screening takes place by two mechanisms: a "wal l effect" and a "turning effect". The wal l effect occurs because o f the fiber-wall interaction. There is a thin "exit layer" near the wal l as shown in Figure 2.3, where Hex is the exit layer height. The concentration in this layer is lower than Chapter 2 Literature review 9 that in the mean flow. The longer and shorter fibers have different concentrations in this layer. Thus the wal l effect alone can result in some screening effect. The turning effect describes the way in which the aperture entry flow affects the passage o f different types o f fibers. For example, long flexible fibers tend to bend and enter the slot and short rigid fibers tend to rotate and to be swept back into the channel. Figure 2.3 Definition o f passage ratio. Kumar (1991) experimentally demonstrated the factors that affect the passage ratio. He found that the passage ratio varied significantly with the ratio of fiber length L to slot width W . The exponential curve in Figure 2.4 (the A-curve) occurs when LIW < 2 and the S-shape curve (the B-curve in Figure 2.4) occurs when LIW > 2. Kumar also showed that under ideal conditions, Chapter 2 Literature review 10 the passage ratio was wel l characterized by a dimensionless penetration number Pe, given by P JJL e~VuL 1.2 Vs/Vu Figure 2.4 Passage ratio for 1mm, 3mm rigid fibers with slot width 0.5mm (Kumar 1991). The above studies about the passage ratio are based on global variables. Olson (1996), and Olson and Kerekes (1997) also used some local variables to express the passage ratio as: 1 j~ v(z) C ( z ) dz pf =T I —— P { Z )T (2-4) where v(z) is fluid velocity at height z above the wal l , C ( z ) is the upstream concentration o f fibers at height z reflecting the "wal l effect" and p(z) is the probability o f fiber passage at height z reflecting the "turning effect". Olson (1996) theoretically estimated the probability o f passage function and experimentally determined C ( z ) as a function of fiber length: Chapter 2 Literature review 11 C{z)_\mzlL if zlL<\lm Cu [l if zlL>\lm where m is experimentally determined as 3.22. The probability of passage function was explored by Olson (1996) using a two-dimensional fiber trajectory simulation with random initial orientations and was approximated by: (l if z < Hex *(H if z .Hex ( 2 6 ) where Hex is the height of the exit layer. Equation (2.6) shows that if the fiber centers are initially within the exit layer, they will pass through the aperture; otherwise they will not pass through the aperture. This expression neglects the effects of the entry flow of different slot shapes. As will be shown later in this thesis, this effect can be important. Based on (2.5) and (2.6), (2.4) is also expressed as (Olson and Wherrett 1998): — if Hex < ml P = \ 2 C c (2-7) 1 otherwise IP where C is the product of two experimentally determined constants and Pe is the penetration number based on the rotor tip velocity Vt, i.e. VW P (2.8) VtL Equation (2.7) shows that for a constant C, P is a function of Pe alone. Lawryshyn (1997) expressed the probability of passage function in the following formula, based on the assumption that the all the fibers are initially along the flow: fl if z < Hcri P ( Z ) = \0 if z>Hcri ( 2 9 ) where Hcri is below the exit layer height. Chapter 2 Literature review 12 The performance equations (2.1), (2.2) were also extended to describe the effects of reject ratio, accept flow rate, contour type, slot width and slot velocities on the fiber fractionation efficiencies (Olson et al. 1998, Olson et al. 2000, Allison and Olson 2000, Olson 2001). 2.4 Screen Plate Design The screen plate is the most important part of the screen and therefore its design is critical for screen performance. The apertures of the screen plates are either holes or slots (Figure 1.2). The size of the aperture is chosen to be small enough to increase the contaminant removal efficiency letting more contaminant go to the reject outlets, but large enough to increase screen capacities allowing more pulp to pass through the aperture. Screens with holes, the first pressure screens used, have typical hole diameters of 1.8mm. They are cheaper, more durable and do better fractionation by fiber length than slotted screens (Olson 2000a). However they tend to increase the consistency in the reject flow, which may cause screen failure, and are much less effective for the removal of cubical debris contaminants. Slotted screen plates have typical slot widths of less than 0.5mm, and are used to remove cubical debris contaminant, but with lower capacity than the screens with circular holes. Contoured screen plates were introduced in the early 1980s (Boettcher 1987, Frejborg 1987, Hopper 1989) since they were found to increase greatly the screen capacity. The feed-side surface of the screen plate has various configurations of ridges and grooves. Various manufacturers have their own different contour patterns (Figure 2.5). The traditional screen plate having no contours near the hole or slot is called the "smooth screen plate". The appearance of the contoured screen plates allows the slot width to be reduced further, increasing the barrier screening. The increase in capacity produced by the contoured screen plate is attributed to the increase of turbulence generated by the slot flow, which breaks up the flocks near the apertures (Halonen et al. 1989). Yu and Defoe (1994a, 1994b) experimentally investigated the flow pattern and the single Chapter 2 Literature review 13 fiber motion near a contoured screen plate. Gooding (1996) has assessed the flow resistance of screen plate apertures using computational fluid dynamics (CFD), flow channel tests and pilot plant trials. Little work has been done on the effects of different contour configurations on the fiber behavior. Kumar (1991) experimentally compared the passage ratio for two heights of the Lehman Contour. Tangsaghasaksri (1994) analyzed the flow characteristics for smooth and contour slots and experimentally examined the effect of some factors on fiber passage. There has been no thorough theoretical (CFD) examination of the fiber motion near different contoured screen slots nor have the effects of the different contoured shapes on the fiber motions been studied numerically. The specific objective of this thesis is to study numerically the fiber motion in contoured slots of different kinds and the effects of the contoured slots on fiber passage, and provide guidance of better slot shape and the importance of contour elements, e.g. contour height, contour upstream and downstream sides. (a) L A M O R T M I C R O V O R T E X C O N T O U R (b) M O D I F I E D F I N C K H C O N T O U R (c) L E H M A N C O N T O U R (d) C A E P R O F I L E C O N T O U R » Flow Z l -•Flow Flow fe, Flo Figure 2.5 Some commercial contour designs (Gooding 1996). Chapter 2 Literature review 14 2.5 Particle Dispersion Two-Phase Flow The fiber-dispersed flow in the screen is a two-phase turbulent flow. In general, there are two approaches to simulate two-phase flow. One is an Eulerian approach, in which the solid phase (particles) is regarded as a continuum and the governing equations are solved for both phases. The other is a Lagrangian approach, which treats the particles as distinct entities within the fluid phase. The Random Walk simulation uses the Lagrangian frame to solve the particle trajectory equations as the particle interacts with a succession of discrete turbulent eddies. Randomly drawn turbulent fluctuations are used to determine the instantaneous velocity of the flow field. Kallio and Reeks (1989) used a Random Walk model to calculate the particle concentration in two-dimensional turbulent flow using a linear Stokes drag and Saffirian lift forces for small spheres. The velocity field was obtained from the log-law and the near wall velocity fluctuation velocities normal to the wall were obtained by assuming a Gaussian distribution with rms values that were curve-fitted from experimental results. Kroger and Drossinos (2000) followed this approach to simulate particle distributions in isothermal and heated turbulent boundary layers. Both the Eulerian and Lagrangian approaches need good predictions of the underlying turbulent flow. Three numerical methods can be used for modeling the turbulence: Direct Numerical Simulation (DNS), Reynolds-Averaged Navier-Stokes (RANS) and Large Eddy Simulation (LES). DNS provides a complete time dependent solution of the Navier-Stokes and continuity equations. DNS assumes that the size of the smallest eddies in any turbulent flow is of the order of the Kolmogorov microscale and specifies the grids so that these eddies can be resolved. Unfortunately, DNS can only be used for very small Reynolds number problems using current computer resources. i Simulations of particles in wall bounded flows using DNS with low Reynolds number have been performed by McLaughlin (1989), Pedinotti et al. (1992), Brooke et al. (1992, 1994), Chen et al. (1995) and Haarlem et al. (1998). Chapter 2 Literature review 15 For RANS methods, the mean velocity and mean pressure are solutions of the Reynolds Averaged Navier-Stokes and continuity equations, the entire spectrum of velocity fluctuations are modeled by a single series of assumptions such as the well-known K-E model. RANS can be used to calculate the particle trajectory in turbulent flow when it is coupled with other tools to get the effect of the individual velocity fluctuation components. The third method to simulate the flow field is LES. In LES, the large scales of motion are calculated directly from the Navier-Stokes equations while the small scales of motion are modeled. LES employs a mesh resolution coarser than that used by DNS but fine enough to resolve the large-scale eddies and uses a subgrid-scale (SGS) model only for those eddies whose size is smaller than the mesh spacing. LES is much more accurate than RANS and it can also be applied to relatively high Reynolds number problems. Wang and Squires (1996) used LES to investigate the particle transport in fully developed turbulent channel flow with Reynolds numbers, Re r , of 180 and 644. Most of the particles considered in previous papers about particle-laden turbulent flow are assumed to be rigid, smooth spheres with a radius that is small compared to the Kolmogorov length scale of the flow. The density of the particles is usually larger than that of the fluid. The particles tend to accumulate near the wall. Brooke et al. (1994) claimed that the near wall concentration buildup comes from particles that arrive in free flight from distances far from the wall and are trapped in the viscous sublayer where normal velocity fluctuations of the fluid are very small. There are only a few theoretical studies of highly elongated particles or fibers in turbulent flow. Kagermann and Kohler (1982) considered non-interacting rigid spheroids suspended in a homogeneous turbulent flow. A renormalized expansion is used to derive a kinetic equation whose coefficients are determined by the Lagrangian correlation function of the turbulent velocity field and by the translational diffusion coefficients and the orientational relaxation frequency of the spheroid. The assumption is made that the particle is smaller than the smallest scales of turbulence and the velocity field is statistically homogeneous, stationary and Gaussian with zero mean. Krushkal and Chapter 2 Literature review 1«3 Gallily (1988) calculated the orientation distribution function of small non-spherical aerosol particles using the Fokker-Planck equation. 2.6 Pulp Fiber Models Jeffery (1922) first investigated the motion of a rigid, neutrally buoyant, ellipsoidal particle in a Stokes flow. The motion equation will be displayed in chapter 4 to compare with the fiber model used in this thesis. Pulp fibers have high aspect ratios and can have considerable flexibility and could not be modeled by Jeffery's equation (Mason 1954). Wherrett (1996) modeled a fiber as a series of spheres based on the work of Yamamoto and Matsuoka (1993, 1994) for 2-D flow. The translation and rotation equations of each sphere determined the fiber motion. The results for a rigid fiber matched Jeffery's theory (1922). Wherrett's computation time was somewhat long, his model was two-dimensional and had no wall interaction model, all of which limit its usefulness for real industrial applications. Olson (1996) developed a theoretical model of fiber passage and verified the predictions of his model experimentally. His model was restricted to rigid fibers. Lawryshyn (1996) introduced a method to calculate the dynamic movement of a flexible fiber in a two-dimensional plane based on Euler-Bernoulli beam bending theory, i.e. small deflection beam theory. Stockie and Green (1998) simulated the motion of flexible pulp fibers using the immersed boundary method, which is able to capture the influence of the fiber on the fluid. Their work was restricted to two-dimensional simulations. In order to simulate industrial applications, a three-dimensional flexible fiber model is necessary, together with a corresponding wall model, since the influence of walls is central to the operation of some of the pulp and paper equipment. In order to account for these characteristics, we have employed a flexible fiber model proposed by Ross and Klingenberg (1997) and developed a corresponding general wall model. The fibers were modeled as linked spheroids Chapter 2 Literature review 17 (not spheres) to save computational time. The details of this fiber model will be described in chapter 4. Chapter 3 Fiber-dispersed turbulent flow simulation 3.1 Introduction In the current simulation, we assume that fibers are affected by the surrounding flow but have no influence on the flow. The flow field could be first simulated and the fibers' motion could then be described using the calculated flow information. A Random Walk method and LES are coupled with the fiber model and wall model to model the fiber motion in the channel upstream of the screen. The fiber model and wall model will be described in chapter 4. For the fiber motion in different shaped slots, only the mean velocities from the RANS are used, because there is no available experimental data for the turbulence intensity needed in the Random Walk simulation, and the current LES solver could not be used for a complicated geometry like a screen slot. Section 3.2 displays the Random Walk method (section 3.2.1) and the Large Eddy Simulation (LES) method (section 3.2.2). The RANS method is included in the introduction of the Random Walk method (section 3.2.1) and will not be displayed separately. 3.2 Fiber-dispersion turbulent flow simulation 3.2.1 Random Walk method In the Random Walk method, the turbulent flow velocity is decomposed into mean and fluctuating parts. 18 Chapter 3 Fiber-dispersed turbulent flow simulation 19 (a) Mean velocity The mean velocities are given by numerical solutions of the Reynolds Averaged Navier-Stokes and continuity equations, which are given by dU, dx, L = 0 (3.1) ac/,. du, I dp i d , du, , „ ^ + u ^ = - - P ^ ; - P ^ - ^ + ^ ( 3- 2 ) where U, is the mean velocity vector, P is the mean pressure, [i is the dynamic viscosity and p is density. Ty = —pujuj is known as the Reynolds-stress tensor, which introduces 6 unknown quantities as a result of Reynolds averaging. /' = 1,2,3 refers to the x (streamwise), y (spanwise) and z (normal) directions, respectively. The usual summation notation applies here. The Reynolds stress tensor needs a model or another equation for closure, and often takes the form of a Boussinesq approximation: Tv=2MTSStt-lPK6v (3.3) where |1T is eddy viscosity, SSy is the mean strain-rate tensor which is given by sso=i y dx j dXj (3.4) and K =-j«,w, is the turbulence kinetic energy. Stj is often called the Kronecker delta, and equals one or zero according to /' equals / or not. The standard K — £ model was used to close the system (3.1) and (3.2). Requirements for the use of this model include that the effective viscosity must be isotropic and that the Reynolds number must be high enough to produce a fully turbulent flow. The eddy viscosity is expressed as | l x = p C u K 2 / e (3.5) Chapter 3 Fiber-dispersed turbulent flow simulation 20 The turbulence kinetic energy K is dK TT die dU, d , . die at dXj dXj oXj dXj (3.6) The dissipation rate e is given by (3.7) Closure coefficients are given as their standard values (Launder and Spalding 1974): C, =1.44,C2 =1.92, CM = 0.09, aK =1.0, aE =1.3 . (3.8) The K — £ turbulent model is not valid in the near-wall region where the influence of laminar viscosity is important. The RANS solver uses the wall-function treatment described by Launder and Spalding (1974), assuming different velocity profiles (log-law when y+ > 11.36, linear when y+ < 11.36) in the near-wall region. The RANS solver used in this thesis was developed by Nowak (1992) at the Department of Mechanical Engineering, University of British Columbia. A finite volume method in conjunction with general curvilinear grids and domain segmentation method is used in the solver. The physical tangential velocity components are used as the dependent variables for the momentum equation and discretization is done by directly using the coordinate-invariant governing equations. The boundary conditions for the channel and for the slot are different and will be discussed in chapter 5 and chapter 6, respectively. (b) Fluctuation velocity Gosman and Ioannides (1981), Shuen et al. (1983, 1985), Adeniji-Fashola et al. (1988) assumed that the turbulence is isotropic and calculated the fluctuations by a Gaussian distribution whose variance is given by Chapter 3 Fiber-dispersed turbulent flow simulation 21 (3.9) where K is the turbulent kinetic energy o f the flow. Because the turbulence field in the channel is anisotropic near the wal l , equation (3.9) may not be accurate. Instead, we use the variances measured by Comte-Bellot 's experiment (1965) for fully developed turbulent channel flow with R e r =2,340 based on the friction velocity (or Rb =57,000 based on the bulk velocity), which is close to the condition simulated here. Unfortunately, there is no data in the region where the channel height z is less than 0.1 times o f the half channel height. In order to fill this gap, we add one point (0,0) at the origin and curve-fit the resulting experimental data using 8 t h order polynomials:M/w T = £io+a1(v//7)+a2(v///)2+...+i%Cy//j)8, where uT is the friction velocity. The fitting coefficients for the three velocities u,v,w are (-1.927e-4, 5.856e+l, -5.609e+2, 2.615e+3, -6.871e+3, 1.064e+4, -9.620e+3,4.698e+3, -9.564e+2), (1.695e-4, 3.403e+l, -3.360e+2,1.678e+3, -4.806e+3, 8.194e+3, -8.213e+3, 4.463e+3, -1.013e+3) and (1.863e-4, 2.270e+l, -2.031e+2, 9.522e+2, -2.595e+3,4.242e+3, -4.097e+3, 2.155e+3, -4758e+2) respectively. Figure 3.1 shows the fitting curve o f the experimental data. 3 u 7 u t 2 .5 2 0 .5 0 0 0 .1 0 .2 0 .3 0 .4 0 .5 0 .6 0 .7 0 .8 0 .9 Z 1(0 .5 H ) Figure 3.1 8 t h order polynomials curve-fitting Comte-Bellot 's experiment (1965). Chapter 3 Fiber-dispersed turbulent flow simulation 22 The vorticity fluctuations are also required for the motion of the fibers. We examined two different approximations: Type I assumes that the total vorticity is given by the formulae: a,; =o.5(— - — ) , co =0.5(— - — ) , « ' = 0 . 5 ( — -— ) > w h e r e w'.v'.w'are the total velocities dz dy ' y dx dz ' 2 3y dx including the mean velocities and the fluctuation velocities. Type II assumes that the vorticity fluctuation components are randomly drawn from a Gaussian probability density of zero mean and £ standard deviation. Wilcox (1988) suggested g for his K - CO model and identified it as 0.09 K" 0.09 K-"the RMS fluctuating vorticity" (Wilcox and Alber 1972). These two types of approximation give similar results. The instantaneous velocity and angular velocity are assumed to influence the fiber motion during a given time period, the residence time, before a new fluctuation component is sampled. The residence time of the fiber in the present eddy is determined by r = min(r,,L, IV), 3 i C 4IC 2 where Tt is turbulence time scale j , L, is the approximate size of a local eddy — , and V = ^u2r + v2r + w2r , where (ur ,vr,wr) is the relative velocity of the fiber with respect to the fluid. The fiber will stay in the same eddy until 7 expires. Although there are many time scales which could be defined for the turbulence, only the turbulence parameters K and e are calculated in the RANS method, so that only f is available for the current calculation as a turbulent time scale. Chapter 3 Fiber-dispersed turbulent flow simulation 23 3 . 2 . 2 L a r g e E d d y S i m u l a t i o n ( L E S ) The equations governing the transport of the large eddies, obtained by filtering the Navier-Stokes equations, are aT° <310> oyuj ^_ du~i _ 1 dp + 1 3 3«, ^ ^ dt J dxj p dx( p dxj dxj 'j C 3- 1 1) where ut is the filtered velocity, p is the filtered pressure , T,- is the subgrid-scale (SGS) stress tensor. Equations (3.10), (3.11) have similar forms to the Reynolds Averaged Navier-Stokes equations (3.1), (3.2). The first subgrid scale stress model proposed by Smagorinsky (1963) has been used in the present work. His model assumes that the SGS stress follows a gradient-diffusion process, similar to molecular motion. Consequently, T{j is given by the eddy viscosity assumption *V - j V " * = -2v^ij (3-12) wherevv =(C J/ JAJ 2 |s | , = ^2SiJSiJ ,StJ = - ( | ^ L + T ^ " ) is the filtered rate of strain tensor and Cs is the Smagorinsky constant (0.1 < Cs < 0.24, Rogallo and Moin 1984). Here = (Ax, * Ax2 * Ax3)3 is filter width, Ax,(/= 1,2,3) is /-direction mesh size. / , =l-exp(-^+ /25) is the damping function needed in the vicinity of the wall. In the current calculation, we set Cs = 0.1. The LES solver was developed by Feng (2001) in the Department of Mechanical Engineering at the University of British Columbia. The governing equations are discretized on a staggered grid using the finite volume method in curvilinear coordinates. The staggered grids avoid oscillations in the computed pressure. A fractional step method (Ferziger and Peric 1996, Kim and Chapter 3 Fiber-dispersed turbulent flow simulation 24 Moin 1985) is employed. The pressure gradient and the incompressibility constraint are integrated implicitly in time. The space derivative is approximated by the second order central scheme. The convective and diffusive terms are treated explicitly by a two-stage Adams-Bashforth scheme. The Poisson-like pressure equation is solved to satisfy the mass continuity. The fast Fourier transform is applied in the direction of periodicity, and the cyclic reduction method is then used to solve the 2D Poisson equation. Parallel computing is used for the LES. The domain is decomposed in one direction and the data is distributed into multiprocessors. The communication occurs only at block boundaries. MPI (Message Passing Interface) is used to exchange the data among the processors. The parallel solver runs well on the PIII/Linux cluster (University of British Columbia) and the MACI Compaq Alpha Cluster (University of Calgary). Chapter 4 Fiber model, wall model and algorithm 4.1 Introduction As described in chapter 2, we have employed a three-dimensional flexible fiber model of Ross and Klingenberg (1997) to account for the flexibility of pulp fibers and to decrease the computation time. Section 4.2 describes this fiber model. Section 4.3 develops a corresponding wall model to deal with the situation when the fiber touches the wall. Section 4.4 shows the calculation procedure of the fiber trajectory. Section 4.5 introduces the algorithms of the fiber model and wall model. Section 4.6 gives some examples validating the fiber model. Section 4.7 is a discussion of the hydrodynamic force in the fiber model. Section 4.8 defines the fiber flexibility and introduces the four fiber types studied in this thesis. 4.2 Fiber model The main assumptions for the fiber model are: • The fiber is infinitely thin • no hydrodynamic interactions between spheroids • Stokes flow assumptions for the hydrodynamic force and torque acting on the spheroids • no wall penetration, ensured by an appropriate force on the fiber normal to the wall • a tangential friction force equal to 0.2 times the normal force was applied to fibers touching the wall In the model, each fiber is modeled as a chain of spheroids linked together by ball and socket joints. The motion of the fiber is determined by solving each spheroid's translation and rotation equations which are derived from Newton's second law and the law of moment of momentum. 25 Chapter 4 Fiber model, wall model and algorithm 26 In order to apply the model to pulp and paper equipment, a wall force is added to the model to deal with the fiber-wall interaction. 4.2.1 Configuration and position of fiber Each fiber consists of N rigid spheroids connected through N — I joints (see Figure 4.1). The lengths a,b(b < a) are the major and minor axes of each spheroid. The fiber can bend and twist much like a real fiber because of the rotational freedom in each joint. The position of the fiber is defined by the center of each spheroid and the orientation of each spheroid. The center of spheroid / (z = 1,2," • ./V) is defined by the Cartesian vector F in a fixed world reference frame £ ( > v ) (see Figure 4.2). The orientation of the spheroid is defined using Euler parameters qu,q2.,q3j, and q4j. Euler parameters are derived from Euler's theorem (Wittenburg 1977). Euler's theorem states that any orientation of body-fixed base e can be obtained by a rotation £ from the reference frame em (which has same origin as the body-fixed frame) about some unit vector T: c . c c c tfi = C O S 2 ' ^ = r * s i n 2 ' ^ = I ^ s i n 2 ' * 4 = r * s m 2 ( 4 1 ) In the model, 9n,92i>93i> ^ #41 orient spheroid 1 relative to the fixed reference frame e , and qu>l2i>V3i> a n d #4, (« = 2,---,AO orient spheroid I relative to the body-fixed frame of spheroid z'-l. The x-axis in a body-fixed frame of spheroid i is along the major axis of this spheroid pointing to joint i+l. Chapter 4 Fiber model, wall model and algorithm 27 4 . 2 . 2 T r a n s f o r m a t i o n m a t r i x f r o m t h e f i b e r b o d y - f i x e d f r a m e t o t h e f i x e d r e f e r e n c e f r a m e The direction cosine matrix G(. which transforms base £ to base is expressed in term o f Euler parameters: q \ + qf2 - 0.5 q l 2 q l 3 - q n q „ q i 2 q i t + q n q i 3 G = Qui* ~ Wn 1t#« + M a 9a + 9M - 0 5 (4.2) The transformation matrix 77V which transforms the fiber frame to the fixed reference frame can be obtained by: Chapter 4 Fiber model, wall model and algorithm 28 TNX=GX, (4.3) TN=TN*G i = 2,---,N 4 . 2 . 3 D y n a m i c s Free-body diagrams are used to apply Newton's law and the law of moment of momentum in order to get the dynamical equations of each spheroid. Figure 4.3 is a free-body diagram for spheroid i. Figure 4.3 Free-body diagram for spheroid i (Ross and Klingenberg 1997). Using the coordinate system of Figure 4.2, Newton's second law is m'ii=%+±S±j (4.4) 7=1 and the law of conservation of momentum takes the form: H , = M . + t ^ . ( c , x X . + Y . ) (4.5) 7=1 where F . is the resultant external force acting through the center of mass, including the direct hydrodynamic force F 1 . ( * ) , indirect hydrodynamic interparticle force F . ( p ) and body force F . ( g ) . Only the hydrodynamic force is considered in the current simulation, since hydrodynamic interactions and particle inertia are neglected. M ; is the corresponding resultant external torque. Chapter 4 Fiber model, wall model and algorithm 29 w,-, I*/, H ; are the mass, translational acceleration, and time rate of change of the angular momentum of ellipsoid / , respectively. X A , X c are the internal constraint forces in joint b and c respectively andY6, Y , are the corresponding resultant internal torques in joint b and c respectively. Bending torque is considered in the current calculation, and is defined by: %*>=kBlfl-<ta>)n (4.6) where k B is a bending constant which takes the form of ^ - , where E is Young's modulus and I is the area moment of inertia. Both rigid and flexible fibers can be modeled by changing this constant. 0; is the angle between the adjacent spheroids, 6 ^ is its equilibrium angle ( the angle for which the bending torque is zero), n is the unit vector normal to the plane of bending. Sy is the element of 0 the connectivity matrix S = -1 1 0 0 -1 1 0 '•. 0 -1 1 0 -1 C.. ( / , / '= ! , • • • , iV ) are a set of body-fixed connectivity vectors which extend from the center of spheroid i to joint / (see Figure 4.2). The hydrodynamic force and torque on the spheroid i are (Kim and Karrila 1991): F? =6nm\xAdidj + YA(SIJ-dJji-V1?-Vj), ~MT = 8/T/AJ3[xAdtclj + YA(STJ -dtdj)]• 0J°] -Q)j)-%K^Yheijkd{dk~E)k , (4.7) where U(p is fluid velocity, and £ w and Q ( o o ) are the rate of strain tensor and vorticity respectively, whose general forms are: du c dv ax ~dy 1 2 aw 2ydy dx 1 (du + dw dz OX; \>E23 ~ T dv aw dz dy dw du dx dz >"3 =" \(du dv 2\dy dx (4.8) Chapter 4 Fiber model, wall model and algorithm 30 where ju is the dynamic viscosity, di is the unit directional vector along the axis. eijk is the component of the alternating tensor and it can be defined by the cross-product vector relationship axb = £ykajbkbeing valid for arbitrary vectors a and b. XA,YA,XC,YC,YH are resistance functions, which are dependent on the prolate spheroid eccentricity e = ~b } where a and b are the major and minor axes of the spheroid. Near-spheres (e<\) and needles (e —>1) are two asymptotic cases of prolate spheroids. The appropriate resistance functions are in listed in Table 4.1, where Lie) = l n f . [l-e) Table 4.1 The resistance functions for a single spheroid (Kim and Karrila 1991). Resistance Functions Prolate-spheroids (0 < e < 1) Near-spheres(e < 1) Needles (e —> 1) £ = yl\-e2 < l , £ ' = (ln(2/£))"' XA ^e3[-2e + (\ + e2)LY 1 2 g 2 1 7 g 4 + 5 175 AE (S-6E)Ee2 6 - 3 £ 12-12TS' + 3£: 2 yA ye 3 [2e + (> 2-l)/:]- , , 3e2 57e4 1 + ••• 10 700 BE 4E2£2 6 + 3E U + 12E + 3E2 xc ^e3(l-e2)[2e-(\-e2 , 6e2 27e4 1 + + ••• 5 175 2e2 (2-2E)eA 3 3£ yC | e3 (2 -e2)[-2e + (\ + e2)/,]"' , 9e2 18e4 1 + + ••• 10 175 2E E2e2 + -+••• 6-3E \2-\2E + 3E2 yH |e5[-2e + (l+e2)Z,|"' £ l _ £ l + ... 2 5 2E {%-5E)Ee2 6 - 3 / i 12-UE + 3E2 Chapter 4 Fiber model, wall model and algorithm 31 The corresponding rotational equation is (4.9) where and N N Qu =sitcf -Af rbl. S £ £ d r A f -C^Af -bu j=l 7=1 k=\ 7=1 A f +F(f +TTT-Af -(A)- -f( Af -dV +F? k=l (4.10) r l —(°°) — (~) N -+^[xcdidj+Yc(Sij-didj)\aJ -%niiJYn'Eijkd,dkEjk 7=1 du = 0 - d y 0 dyjj -dx;; ij — ij dx ; ; 0 > bii = 0 b Z y - b y - b z ; i byij -bX:; IJ IJ bx ; i 0 (4.11) The translation equation of motion for each spheroid is — N _ Vi=Ya>jXbji-(A)-1-7=1 where = -j7tab2Apg is the body force that is ignored in the current calculation, because the density difference between the fiber and the water Ap is very small. Fjw) is the force due to the wall. If the fiber does not touch the wall, it is zero. If the fiber touches the wall, this force will be calculated. The details of this calculation are shown in section 4.3. The matrices dy and bjj are given by Zj=-&c)ijXj=-(CTcV)ij (4-13) N N / _ \ N =1 k=l 7=1 Ah) -(• •Uj ') - (g) >(w) + Fj +FTj (4.12) Chapter 4 Fiber model, wall model and algorithm 32 where C is an NxN matrix with components Cy=S0Cy ,vy =StJ -UN, and the connectivity matrix r. is defined as: T = -1 -1 -1 0 -1 -1 -1 ... _ i ... _ i ... _ i 0 -1 (4.14) After angular velocities C0i (i = !,•••,N) sue calculated, the relative angular velocity gt (i = !,•••,N) can be obtained through: N (4.15) Euler parameters can be updated according to the following relationship: <i\i The Euler parameters are normalized via 0 In 0 Iii 0 (4.16) 0 3*. 9 a = V 2 2 2 2 $1/+02, +03/+ 04, , 7=1,2,3,4, i = l,...,N (4.17) 4.3 Wall model Fibers frequently touch the wall in pulp and paper equipment, so that a wall model that can efficiently deal with the fiber-wall interaction needs to be developed. For the smooth slot geometry, Olson (1996) developed a two-dimensional wall model based on the assumption that a reaction force normal to the wall is exerted on the fiber to stop it passing through the solid wall, and the friction Chapter 4 Fiber model, wall model and algorithm 33 force tangential to the wall is proportional to the normal force on the fiber. Lawryshyn (1997) used a similar approach. Because pulp and paper equipment can have many different wall geometries, a three-dimensional universal wall model has been developed, which can deal with the fiber interaction in any wall geometry. Because the fiber in the current simulation consists of several spheroids, the fiber-wall interaction depends on the spheroid-wall interaction, which makes the current wall model more complicated. Section 4.3.1 introduces the way we judge how the fiber touches the wall. Section 4.3.2 gives the equations of the fiber when it touches the wall. 4 .3 .1 T h e c o n d i t i o n f o r t h e f i b e r t o t o u c h t h e w a l l The wall surface is composed of grid cells. If the grid used for the flow calculation is coarse, it can be refined into equal intervals in any direction. Figure 4.4 is a schematic diagram showing the fiber touching the wall. The fiber touches the wall if the system of the single spheroid and the refined grid line segment has a solution. The equation of the refined grid line segment is X — XQ ~h (X\ (X\ XQ ), y = y0+a1(y1-y0), 0<a] ^ z = z 0 + « , ( z 1 - z 0 ) (4.18) where (x0, y0, z0), (JC, , y,, z,) are the position of the two vertices of the line segment. The equation of the spheroid in fiber coordinates is 2 2 2 a2 b2 b2 (4.19) If the transformation matrix which transforms the fixed world coordinate to the fiber coordinate is TR, then 77? = 77V~', where 77V is defined in (4.3). (xc,yc,zc) is the fiber center in the fixed reference coordinate. The position in the fiber coordinate becomes xf f Vf = tr(1,2) tr(l,3)Yx-x ^ tr(2,\) tr(2,2) tr(2,?>) y-yc (4.20) Chapter 4 Fiber model, wall model and algorithm 34 So the spheroid equation becomes [tr(U)(x-xc)+tr(l2)(y-yc)+trimz-zc)]2 [trimx-xc)+tr(2,2)(y-yc)+H2,3)(z-zc)]2 b2 | [/K3,l)(x-x c)+?K3,2)(y-v c)+/K3,3)(2-z c)] 2 ft2 If the system of (4.18) and (4.21) has a solution, then we say the fiber touches the wall. F 2 Figure 4.4 Schematic diagram of the fiber touching the wall. 4.3.2 T h e f o r c e e q u a t i o n f o r t h e s p h e r o i d w h e n t o u c h i n g t h e w a l l When the fiber touches the wall, the assumption for the forces added to the fiber is that the normal force is always perpendicular to the wall, and the tangential friction force is proportional to the normal force. If the i th-spheroid touches the wall, then the normal force for this spheroid is: Chapter 4 Fiber model, wall model and algorithm 35 F,=fc,-f,= fy where fc, = fy is the unit normal vector that is generated from the normal to the contacted wall surface and f is the magnitude of the normal force. The unit tangential vector ty is: ty \ t z J t tv Ivel (4.22) where V (fc) tv = V unori • fy Jwy UJ (fa) unori = V • fy i UJ (4.23) and u, v, w are the ith-spheroid velocities. The coefficient of the total force added to the / th-spheroid is (ft*) (fa) fty = fy + TJ- ty J t z , i [fa) f z , (4-24) where n is the coefficient of wall friction that is taken as 0.2 in all the simulations. The added total force is FT, = fty yftzj f, (4.25) Chapter 4 Fiber model, wall model and algorithm 36 where / , needs to be calculated. Introducing (4.25) into (4.10) and (4.12), the new velocities consist of two parts: Vt=VnlHl+V, old i wall i 7 (4.26) where VM is exactly the form of (4.12) without wall contact. Vwall comes from the fiber-wall contact, which takes the form of 7=1 fYAf\o)wallkxbkJ 7=1 yt=l 7=1 (4.27) o)wall . is obtained from the similar form of (4.10) with the same Qu and different Z), 3!=-£d,ixTf7'-A?>-(Ar±(^ 7=1 L *=' ^ After applying (4.25) to (4.27), (4.28), (4.26) can also be written as (4.28) (4.29) where W, takes / , from (4.27). Taking the product of (4.29) and the unit normal vector fc- = fy \fZJ , we get: fci*Vi=fci*Voldi+fci*Wi-fi (4.30) The LHS of (4.30) must be zero according to the assumption that the fiber can not go through the wall. fct • VM . and fc, • VVt could be calculated from the above analysis. If /' = 1, (4.30) becomes a single equation with the unknown / , ; if / > 2, (4.30) becomes a system of equations with the unknowns fni>2. After /) is solved, the rotational and translational equations (4.9) and (4.12) Chapter 4 Fiber model, wall model and algorithm 37 are calculated to predict the motion o f the fiber when it touches the wal l . A n example is shown in Figure 4.5. The fiber passes through the slot after contacting the slot wal l . Fiber passes through the slot after contacting slot wall -0.0025 0.0025 Figure 4.5 A n example o f the fiber touching the wall . 4.4 Calculation procedure The fiber model together with the wall model needs to be coupled with the flow model in order to predict the fiber trajectories. Each spheroid's dynamic equations are solved. The flow model imposes forces and torques on individual spheroids so that the fiber may experience different forces in its different parts. . The whole procedure to track a fiber is: (1). Initialization o f each spheroid's position, velocities and Euler parameters ( la) . The fiber's position is given by three variables: fiber center Pc, polar angle 6, which is the angle between the fiber main axis and Z-axis, and the azimuthal angle </>, which is the angle Chapter 4 Fiber model, wall model and algorithm 38 between the Y-axis and the projection of the fiber main axis on XY plane (see Figure 4.6). Each spheroid's position is calculated through the following procedure: h = P c - N - a - ^ n k P\ =PQ +<* u^nit (4.31) P^P^+la-d^, i = 2,--,N where P0 is the position of fiber front end and dmit = sinc9sin0 sinS c o s </> c o s t ? is the unit direction of the fiber main axis. (lb). Initial velocities for each spheroid are initially set to the fluid velocities at the center of each spheroid. The fluid velocities at the center of each spheroid are obtained by interpolating the velocities at the corners of the cell in which the center of the spheroid lies. (lc). Initial Euler parameters can be calculated from (4.2) if transformation G. (/ = 1, • • •, N) is known. G, transforms the first spheroid frame x-y-z to fixed reference frame X-Y-Z. The elements of G, are direction cosines of the three axes of x-y-z. The direction cosine of the x-axis drx is just dmit in (la). The equation of the plane perpendicular to the *-axis through P0 is: dfx»(R-OP0) = 0 (4.32) dr can be chosen from the above plane as long as (4.32) is satisfied. d?z can be found from: y dr =df XdP (4.33) G. (Z = 2,- • N ) transforms the frames of spheroids / to frame of spheroid i-l. Because initially the frames of all the spheroids are equal, so G ; (i = 2, • • •, N ) is a unit matrix. Chapter 4 Fiber model, wall model and algorithm 39 X Figure 4.6 Fiber's initial position. (2) . Obtain the fluid velocities U(°°], the rate o f strain tensor E °° and vorticity £2*°°' at the center o f each spheroid, by interpolating the corresponding values at the corners o f the cell in which the spheroid centers lie. (3) . Calculate d.. and b.. using (4.13). (4) . Transform d. and b. to a fixed reference frame using the transformation matrices (4.3). Chapter 4 Fiber model, wall model and algorithm 40 (5). Using table 4.1, find the resistance tensors: 67U/la[xAd.dj + YA(Sjj,- dtdj)\ Snpa3[xAdldJ + YA(StJ - d,dj)\ and - Snjua3YH£i (6) . Calculate the bending torque using (4.6). (7) . Compute angular velocities from (4.9) and translation velocities from (4.12). (8) . Calculate the time derivatives of the Euler parameters from (4.16) using the relative angular velocities in (4.15). (9) . According to the translation velocities from (7), judge if the spheroid touches the wall in the next time step using the methods described in section 4.3. If the spheroid touches the wall, calculate the wall force and go to (10a) to get the next step spheroid position and Euler parameters; if several spheroids touch the wall at the same time, a system equation is solved to find the wall force for each wall-touched spheroid. The method of wall force calculation is described in section 4.3. If the spheroid does not touch the wall, go to (10b) to get the new spheroid position and Euler parameters. (10) . Different update procedure (10a). Update spheroid position and Euler parameters using first Euler method as shown in section 4.5. (10b). Update fiber position and Euler parameters using a Variable order Runge-Kutta Method as shown in section 4.5. (11) . Repeat (2) to (10) to get the fiber trajectory in the flow field. 4.5 Algorithms The translation equation (4.12) is a first-order ordinary differential equation in spite of its complicated form. The general form of the first-order ordinary differential equation is: Chapter 4 Fiber model, wall model and algorithm 41 y =f(t,y), y(0) = v (4.34) Two different algorithms have been used for discretizing (4.34) during the course of the research. Section 4.5.1, 4,5.2 introduce the Euler method and the Variable Order Runge-Kutta Method, separately. 4.5.1 Euler method The Euler method for solving (4.34) is: where h is the constant time step. This Euler method is simple and easy to implement, but it is only a first-order scheme and the time step is also fixed. For the fiber motion in the present pulp and paper equipment, a more accurate solution is needed. The computation time is also expected to be large for a large numbers of fibers. So higher order methods with a variable time step are needed. 4.5.2 A Variable order Runge-Kutta Method (VRK) The fiber frequently touches the wall in the present equipment. The right hand side of (4.34) changes rapidly in the near wall region. In this region, a higher order algorithm may have poor performance (see comments in Cash and Karp 1990). It is therefore necessary to use variable order methods for different regions. The variable order Runge-Kutta Method is a family of explicit Runge-Kutta formulas. Each member of the family consists of a fifth-order formula that includes imbedded formulas of orders from 1 to 4. A proper order formula is chosen by calculating the solution at several different orders before the full Runge-Kutta step is computed. An error tolerance parameter is specified as 10"9 in the code. Value of 10"8 and 10"10 were also tried but no significant change was found. The detailed algorithm and the comparison with other algorithms are included in the work of Cash and Karp (1990). Chapter 4 Fiber model, wall model and algorithm 42 4.6 Validation example 4 .6 .1 F i b e r m o v i n g i n a s i m p l e s h e a r flow According to Jeffery's analysis (1922), for a prolate spheroid with aspect ratio re in an ambient simple shear flow U = (Gy, 0, 0) where G is a constant, the angular motion o f the spheroid is described by: t a n f l = i , & — — and tan = tan(2>T •£-), (4.36) •^re cos 0+sin S 'p where 9 is the angle between the fiber's major axis and the vorticity axis (Z axis), 9 is the angle between the Y-axis and the XY-projection o f the fiber axis, Tp is the orbit period and C is the orbit constant, determined by the initial orientation, C = tan<90Jcos2 d0 +-^-sin2 <p0. (4.37) Figure 4.7 Fiber rotates in X Y plane. Figure 4.8 Fiber translates and rotates in X Y plane. Chapter 4 Fiber model, wall model and algorithm 43 This analysis predicts that a prolate spheroid w i l l repeatedly rotate through the same orbit, and the orbit period is independent o f the initial orientation or orbit constant C . Figure 4.7 shows the fiber motion in the X Y plane for a non-zero velocity field u = y, where v is the second coordinate o f each spheroid's position. The fiber rotates in the X Y plane because the fiber is initially set in the X Y plane and the fiber center is set at the origin. If the fiber center is not set at the origin, then the model predicts the translation and rotation in the X Y plane, as demonstrated in Figure 4.8. Furthermore, i f the fiber is injected randomly in physical space, the fiber w i l l translate and rotate in 3-D, as shown in Figure 4.9 for a rigid fiber, and in Figure 4.10 for a flexible fiber. The calculated orbit compares closely to Jeffery's analysis, as shown for a rigid fiber in Figure 4.11. Figure 4.9 Fiber translates and rotates in 3D. Figure 4.10 Flexible fiber translates and rotates in 3D. Chapter 4 Fiber model, wall model and algorithm 44 5 sphero ids , a spec t ratio = 2 eo=7t/5o Figure 4.11 Comparison between flexible model and Jeffery's prediction. In a simple shear flow, Skjetne et al. (1997) showed some other examples o f a single fiber trajectory calculated by the current fiber model and comparison to Jeffery's analysis (1922). 4 . 6 . 2 A m o r e g e n e r a l e x a m p l e Here, the fiber model and algorithm are verified by a more general example. Olson's (1996) analytical solution is given for the fiber moving in extensional flow that is similar to the flow in the paper machine headbox. The analytical solution and the fiber motion calculated by the fiber model and the algorithm are compared. The C P U time for the fiber motion with the two different algorithms is also shown. Chapter 4 Fiber model, wall model and algorithm 45 4.6.2.1 Analytical solution The flow field of the symmetric 2D convergent section of the headbox (Figure 4.12) (Olson 1996) is approximated by U(x) = - U0 V(x) = o -y (4.38) Figure 4.12 Geometry of the symmetrical convergent section of a headbox. By assuming that the center of the fiber moves along the centerline, Olson (1996) determined the analytical alignment of the fiber to be : 6{x) = arctanJ sin(20o) r cos(20o) + l ^ 1 - 2 — + (—) 2 V x J r/+l (4.39) where 6 = 00 at t = 0. Chapter 4 Fiber model, wall model and algorithm 46 4.6.2.2 Numerical solution In the example, the contraction angle ft = 16.7°, the inlet channel height is 150mm and the outlet channel height is 15mm. So the contraction length is 250 mm. Taking the initial alignment to be 60 = 60° , Figure 4.13 shows a comparison o f the analytical solution with the simulation angle from the Variable order Runge-Kutta Method ( V R K ) . Figure 4.14 is the comparison o f the relative error o f the alignment with V R K and first order Euler for different time step dt. For V R K , the maximum time step is 3e-2 seconds, and the minimum time step is 7e-4 seconds. From Figure 4.14, for the Euler method, the bigger time step causes bigger relative errors. If the Euler algorithm takes the minimum time step of V R K (7e-4), the relative error is much greater than that o f V R K . I f we want a solution by Euler to be o f the same accuracy as that o f V R K , a minimum time step le-6 seconds must be taken. Table 4.2 is a comparison o f C P U time for V R K and Euler with different time steps. It shows that i f the 1 s t order Euler algorithm is used to get the solution that is as accurate as that o f V R K , a much smaller time step must be taken and the C P U time is increased by about 4000 times. For these reasons, V R K is used in all following examples and calculations. 50 4 0 30 20 1 0 0 0 _l I I I I I I I I I I I I I I I I I L 0.05 0.1 0.1 5 x position along the co n ve rg e n t se ctio n (m e te r) 0.2 Figure 4.13 Comparison o f analytical solution and V R K solution. Chapter 4 Fiber model, wall model and algorithm 47 0.03 0 .025 0.02 § Q) 9> 0 .015 (0 0.01 0 .005 ~i 1 1 r 0.05 1 1 r 0.1 - 1 — i r 0.15 ' 1 0.2 - 1 — r VRK — — — — E u l e r dt=1 e-3 E u l e r dt= 7e-4 E u l e r dt= 1 e-4 E u le r dt= 1 e-5 E u l e r d t = 1 e - 6 / / / / / / • / / / / s / s s 1 h i: I-I I I / -i s / _l I I I I I I I I I ' I ' ' I I , , 0.03 0 .025 0.02 0.01 5 0.01 0 .005 0 0.05 0.1 0.1 5 0.2 x position along the convergent section (meter) Figure 4.14 Comparison o f the relative error o f V R K and the first Euler method. Table 4.2 Comparison of C P U time for Euler and V R K methods. Algorithm Dt(s) CPU time (s) V R K 7e-4< dt <3e-2 0.19 le-3 0.66 7e-4 1.37 E U L E R I e-4 7.65 I e-5 80.24 le-6 808.61 4.7 Discussion of the hydrodynamic viscous force The resistance force and torque equations (4.7) were based on the Stokes flow assumption ( K i m and Karri la 1991). The velocity fields were expressed as a weighted integral o f the Green's Chapter 4 Fiber model, wall model and algorithm 48 function. A multi-pole expansion was taken for the Green's function. The simplified solution was obtained by applying the multi-pole expansion for the particle having a high degree of symmetry, e.g. the sphere and the spheroid. If the spheroid becomes a sphere, a = b, and e = 0 . From table 4.1, XA = 1, YA = 1, Y" = 0 and (4.7) is reduced to — — ( - ) F i =67tna(Uj - Vj), —•(*) , .•—•(-) — . Mi =$7Ujua -(Glj -CDj) (4.40) which is the general form of force and torque acting on a sphere, displayed in many general fluid mechanics textbooks. For the case of the linear flow field, Kim and Karrila (1991) showed that " as the spheroid rotates in the linear field, its axis traces a periodic trajectory known as a Jeffery orbit". 4.8 Fiber flexibility Fiber flexibility Ff is an important fiber property affecting the pulp and paper process. It is defined as SlUr EI ''iff where Sliff is fiber stiffness, defined as the product of elastic modulus Em and the area moment of inertia / . Four types of fiber with different length and flexibility are studied in this thesis: 1mm rigid fiber: EJ = 329xl0~12 Nm2, lmm flexible fiber: £'m/ = 3.29xl0"12 Nm2, 3mm rigid fiber: Eml = 329xlu"12 Nm2, 3mm flexible fiber: EJ = 3.29xl(T12 Nm2. They are similar to those examined by Kumar (1991). Chapter 5 Numerical simulation of fiber concentration 5.1 Introduction Understanding the fiber motion and fiber concentration in the flow approaching the screen plate is a key part of predicting passage ratio and efficiency of the screen. From the literature of numerical simulations in chapter 2, most of the particles previously studied in particle-laden turbulent flow are assumed to be rigid spheres, small in size compared to the Kolmogorov length scale of the flow and with density larger than that of the fluid. This is not the case for the pulp fibers. For softwoods, the average fiber length is 2-5 mm with width 30-45 ju m; for hardwoods, the fibers can be 1-2 mm long and 20 fi m wide. So the aspect ratio is as high as 50-100 (Kumar 1991). Pulp fibers also have considerable flexibility. The objective of this chapter is to study numerically the fiber concentration in a channel flow similar to the flow approaching the screen slot. Several experimental measurements have been carried out of the fiber concentration in turbulent flow similar to that approaching a screen. Gooding (1986), Gooding and Kerekes (1989) investigated the motion of 3mm fibers at five heights near a wall and found that there are fewer fibers in the "exit layer" than in the mainstream flow. Olson (1996) experimentally determined the detailed distribution of fiber concentration for a range of fiber lengths. The channel velocity in his experiment was 7.1 m/s and the channel Reynolds number Rb was about 71,000, based on the channel bulk velocity and half channel height. Dyed nylon fibers of length 3.1, 2.0 and 1.1 mm with diameter 40 microns were used. The fiber concentration was less than 5000 fibers/liter so that the fiber-fiber interaction was negligible. In this chapter, the fiber concentration is simulated by both Random Walk and LES methods, and is compared to Olson's (1996) experimental measurements. 49 Chapter 5 Numerical simulation of fiber concentration 50 We assume that the fiber suspension is dilute and we ignore the fiber-fiber interaction and the effects of the fiber on the turbulent flow. "Dilute" means Nc < 1, where Nc is called "crowding factor" representing the number of fibers in the volume swept out by the length of a single fiber (Kerekes and Schell 1992). If a spherical particle of density p^and radius a moves in a fluid of molecular viscosity ju and friction velocity uT, the Stokes number (Crowe et. al. 1985) is defined as St = Y~ : the ratio of the particle response time rp = 2p9p° and some fluid time scale rf = -— (Rouson and Eaton 1994, H is the channel width). If the Stokes number is very small, the particle follows the flow. If it is very large, the particle will not be affected by turbulent fluctuations. If the Stokes number is in the order of one, the particle will be affected by the turbulence but cannot be the flow tracer (Fessler, et. al., 1994). In the current simulation, if we take a as the typical long radius of the spheroid, the response time is about 0.05 s and St is about 1.9. So the fiber is expected to be affected by the turbulence but it does not exactly follow the flow. Section 5.2 describes the boundary conditions and grid independence of the two flow simulation methods that were introduced in chapter 3. Section 5.3 examines the fiber motion and concentration in turbulent mean flow. Section 5.4 explores the fiber concentration in turbulent flow using Random Walk and LES methods. Section 5.5 is a discussion of the effects of some factors on the fiber concentration. A brief summary in section 5.6 concludes this chapter. 5.2 Flow simulation This section describes the boundary conditions and the degree of grid independence for the two flow simulation methods that were used for the concentration calculation. The first method is the Reynolds Averaged Navier-Stokes (RANS) method (section 5.2.1) and the second is the Large Eddy Simulation (LES) method (section 5.2.2). Chapter 5 Numerical simulation of fiber concentration 51 The height of the channel is chosen to be 0.02m (z direction) in both RANS and LES calculations, which is the about the same as that of Olson (1996). In the RANS calculation, the length of the channel is 0.4m (x direction) and the spanwise size (y direction) is 0.01m. In the LES calculation, the length of the channel is 0.048m (or 0.060m) and the span-wise size is 0.02857m. The Reynolds number Rb is 65,000 for RANS, and for LES, Rb is 61,000 and RT is 2,820 based on the friction velocity and half channel width. The Reynolds number for LES can not be specified exactly at the start of the calculation and it is found only when the calculation is complete. Those flow characters are similar to the upstream flow condition of an industrial screen (Kumar 1991, Olson 1996, Gooding 1996). 5 . 2 . 1 R e y n o l d s A v e r a g e d N a v i e r - S t o k e s ( R A N S ) m e t h o d Boundary conditions: The boundary conditions are specified as following (see Figure 5.1): inlet Vu - 6.5 m/s K , £ given top bottom X outlet = o Figure 5.1 Boundary condition for RANS of the channel flow. Chapter 5 Numerical simulation of fiber concentration 52 Inlet: uniform velocity Vs = 6.5 m/s which is similar to the upstream average velocity of industrial screen (Kumar 1991, Olson 1996, Gooding 1996). The turbulent kinetic energy and the dissipation rate of turbulent kinetic energy at entry are specified as: k = l.5*(I*Vu)2 c>> ( 5 1 ) £ = — where 1 is turbulence intensity taken as 0.05, Ls is turbulence length scale taken as GCtH 12, where Ctt is a constant equal to 0.03. The boundary values for k and e have no significant effect on the downstream flow field. Outlet: "gradient condition", which means the normal gradients of all the variables (velocities, pressure and turbulent scalars) are zero. Top and bottom: no-slip wall conditions, all the velocities are set to zero. Spanwise direction: symmetry boundary conditions are used, which means the mass, momentum and scalar fluxes are zero. Grid independence: A stretching grid generation technique (Moin and Kim 1982) was used to create 100 non-uniform grids in the z direction. The transformation that gives the location of grid points in this direction is: Zj = 0.5H (tanh(^ artanh A)l A + \) (5.2) where = - l + 2(y-l)/(N2-l), j = \,2,---,N2, N2 is the total number of grid points in the z direction and A is the adjustable parameter (0 < A < 1). Larger A results in finer grids near the wall. Three kinds of grids in the z direction are used: ztypel: A = 0.95 and the first grid size is 0.039mm; ztype2: A = 0.98846 and the first grid size is 0.013mm; ztype3: uniform grids 0.2mm. Chapter 5 Numerical simulation of fiber concentration 53 For grid ztypel, two different kinds of grids in the x direction are used (0.75mm uniform for Rancasel and 0.5mm uniform for Rancase2). For grid ztype2, 0.75mm uniform grid and even larger non-uniform grids are used in Rancase3 and Rancase4, separately. For grid z_type3, 0.75mm uniform grids are used for Rancase5. Table 5.1 shows the cases with different grids in the x and z directions. Table 5.2 includes the ratio of the fiber length to the grid size in the x direction. Figure 5.3 is the grids of Rancasel. Table 5.1 Different cases for RANS. CASE Zgrid Xgrid Rancasel Ztypel (non-uniform) Uniform, 0.75mm Rancase2 Ztypel Uniform, 0.50mm Rancase3 Ztype2 (finer, non-uniform) Uniform, 0.75mm Rancase4 Ztype2 Non-Uniform, average 2mm Rancase5 Ztype3 (uniform) Uniform, 0.75mm Table 5.2 Ratio of fiber length to grid size in x direction for RANS. 1mm 2mm 3mm Rancasel,3,5 (0.75mm) 1.333 2.667 4.0 Rancase2 (0.5mm) 2.0 4.0 6.0 Chapter 5 Numerical simulation of fiber concentration 54 0.02 0.015 z o . o i B 0.005 -0.01 J L 1 ' 1 , 1 1 0.01 Figure 5.2 Grids o f Rancasel: non-uniform in z direction and uniform in x direction. Figure 5.3 is a comparison o f u+ from R A N S for different cases with Spalding's single composite formula (1961) that covered the entire wall-region with kx = 0.40, B = 5.0: z =u +e (5.3) The results for Rancasel and Rancase2 are about the same, and the results for Rancase3 and Rancase4 are about the same. These show that the grids in x direction do not matter for the mean velocity. It is easy to understand that i f the flow is fully-developed, the mean flow field should not change much for different grid sizes. Obvious differences appear for cases with different grids in the z direction. Uniform grids in Rancase5 give the results closest to the analytical solution. The finest grids near the wal l in Rancase3 and Rancase4 give the least accurate results. This could be explained by the following: the standard K-e model produces inaccurate uT very close to the wal l , e.g. for Chapter 5 Numerical simulation of fiber concentration 55 Rancasel and Rancase2 ( z + for the first node is about 13), uT =0.33 ; for Rancase3 and Rancase4 (for the first node z + =5), uT =0.35. But Rancase5 produces uT =0.29 which is probably most accurate because K - E is very powerful for relatively high z + ( z + =58 for Rancase5). So u+ is closest to the theoretical solution for Rancase5 but it is farthest from the theoretical solution for Rancase3 and Rancase4. But Figure 5.4 shows the mean velocities are very close to each other for those five different cases. The maximum difference occurs between Rancase3 (4) and Rancase5 where the relative error is under 6%. So we assume the five types o f grids w i l l not significantly affect the fiber behavior and only the grids in Rancasel are used for the Random Walk calculation o f fiber concentration. zplus Figure 5.3 Comparison of mean velocities u from R A N S for different cases. Chapter 5 Numerical simulation of fiber concentration 56 12 lo E XT 10 £ .O) '33 .c •3 8 c c to (A (/> 2 o re o o o > Rancasel Rancase2 RancaseS Rancase4 Rancase5 0.005 0.01 0.015 Channe l height (m) 0.02 Figure 5.4 u velocity across the channel for different cases. 5 . 2 . 2 L a r g e E d d y S i m u l a t i o n ( L E S ) Boundary conditions: The boundary conditions for L E S are as follows: periodic conditions are used for the inlet and outlet (x direction), south and north (y direction). W a l l conditions are set at both the top and bottom walls. Grid independence The same grids in the z (wall) direction (100 non-uniform grids of ztypel and ztype2) as that for the R A N S study are used in the z direction to test the grid independence. Two kinds o f uniform grids in the x direction (0.6mm and 0.75mm) are used. The uniform grid in the y direction is 0.57 mm. Table 5.3 shows the three test cases. Lescasel used ztypel in the z direction and 0.75mm Chapter 5 Numerical simulation of fiber concentration 57 uniform grids in the x direction. Lescase2 and Lescase3 use the same kind o f grids in the z direction (ztype2), but different grids in the x direction (0.75mm for Lescase2 and 0.60mm for Lescase3). Table 5.3 Different cases for Large E d d y S imula t ion calculat ion. C A S E Zgrid Xgrid Lescasel Ztypel Uniform, 0.75mm Lescase2 Ztype2 Uniform, 0.75mm Lescase3 Ztype2 Uniform, 0.60mm — Liu et al.'s experiment zplus Figure 5.5 Comparison o f w + vs. z + for L E S of different cases with Spalding's theoretical solution and three experimental data sets: L i u et al. ( R e T = 184), Comte-Bellot ( R e r = 2340) and W e i & W i l l m a r t h ( R e r =1655). Chapter 5 Numerical simulation of fiber concentration 58 Figure 5.5 displays the comparison of u+ vs. z+ from LES of three cases in Tables 5.3 with theoretical and experimental solutions. The theoretical solution is Spalding's formula (5.3). Three experimental data sets are from Liu et al.'s experiment (Rer=184, 1991), Comte-Bellot's experiment (ReT =2340, 1965) and Wei & Willmarth's experiment (ReT =1655, 1989). The results for Lescase2 and Lescase3 are nearly the same, which shows the similar phenomena as RANS: the grids in x direction do not matter too much. The result for Lescasel is a little lower than that for Lescase2 and Lescase3 for z+ > 20 . The region where the LES results and other results differ most is the buffer region. This is an area which is hard to predict. From Figure 5.5 we can observe a big difference between the simulations of RANS and LES: finer grids in the z direction produce worse results for RANS (Figure 5.3), but results in better results for LES (Figure 5.5). This demonstrates the difficulties of the K - e model and the power of LES near the wall, because the K - e model requires the first grid must not be too small (z+ >30) but there is no grid limit for the LES. Figure 5.6, 5.7 and 5.8 show the three fluctuation velocities calculated by LES for the same three cases in Table 5.3 compared with Liu et al., Comte-Bellot and Wei & Willmarth's experiments. The three cases of grids in Table 5.3 produce almost the same results. Considering the performances in these figures, the grid of Lescase3 was chosen to calculate the fiber concentration. We also see that there are some differences between the simulated fluctuations and the experimental data. The Reynolds numbers for the experiments and the simulations are not exactly the same, and also the LES data includes only the contributions from the resolved eddy motions. Chapter 5 Numerical simulation of fiber concentration 59 0.5 — <3- — L iu et a l . 's e x p e r i m e n t L e s c a s e 2 : z-z ty p e 2 , x-0 .7 5 m m L e s c a s e 3 : z -z ty pe 2 , x - 0 . 6 0 m m L e s c a s e l : z - z t y p e l , x - 0 . 7 5 m m ^ C o m te-B e Mot's e x p e rim e n t £> W e i a n d W illm a r t h ' s e x p e rim e n t 0 ' — 1 — 1 — 1 L 1 1 i i i i i i_ J I I I L. 0.1 0.2 0.3 z/H 0.4 0.5 Figure 5.6 u'/ur calculated by L E S of three cases, compared with the experiments of: L i u et al . ( R e r =184), Comte -Be l lo t (Re r = 2340) and W e i & Willmarth ( R e r =1655). 1 .5 ro 0.5 3 = L e s c a s e 2 : z-z ty p e 2 , x-0 .7 5 m m L e s c a s e 3 : z-z type 2 , x - 0 . 6 0 m m L e s c a s e l : z - z t y p e l , x - 0 . 7 5 m m ^ C o m te-B e Dot's e x p e rim e nt Figure 5.7 v'/uT calculated by L E S of three cases, compared with the experiments o f Comte-Bellot ( R e T = 2340) (other experiments are not available). Chapter 5 Numerical simulation of fiber concentration 60 z/H Figure 5.8 Wlut calculated by LES of three cases, compared with experiments of: Liu et al. ( R e r = 184), Comte-Bellot ( R e r = 2340) and Wei & Willmarth ( R e r = 1655). 5.3 Fiber motion and concentration in turbulent mean flow In this section, we w i l l examine the fiber motion and concentration in turbulent mean flow without fluctuation. Figure 5.9 shows the mean velocity profile in main xz plane calculated by R A N S method. The flow is fully developed for about x > -0.1m. Chapter 5 Numerical simulation of fiber concentration 61 0 . 0 1 5 N 0 .0 1 0 . 0 0 5 m m "V I t " \ a t a r . is? L» Sr, ,.A a t ig . 2 - 0 . 1 5 0 . 1 -0 .0 5 0 X 0 . 0 5 0 . 1 0 . 1 5 0 . Figure 5.9 The mean velocity profile in xz plane. If a 1 mm fiber is initially set in the xz plane (6 is 70 0 and <p is 90° in Figure 4.6) and its initial height is 1mm, then the fiber w i l l rotate and translate along its center with nearly the same height as its initial height. If the same fiber is initially put at 0.5mm, it also rotates and translates in xz plane keeping nearly the same height 0.5mm. But i f the fiber is initially set below the half fiber length (for example 0.2mm), it w i l l rotate and translate and move up to 0.5mm and then keep this height. Figure 5.10 shows the centers o f fibers with these three initial heights. So there are no fibers below half fiber length 0.5mm. If the concentration is based on the fibers only in the xz plane, it must become zero below the height o f the half fiber length. But i f the fiber is initially set in three-dimensional space, the results are different: Figure 5.11 shows the motion o f fibers with the same angle 0 and different angles (j). The fibers move to different heights. The bigger the <p, the higher the fiber moves. The effect o f angle 6 is smaller than that o f <p: Figure 5.12 shows the centers o f the fibers with the same angle (f) and different angle 6. Chapter 5 Numerical simulation of fiber concentration 62 50 100 153 I i i I | i i i I | i i i I | I I I Lo • • • • • • • • • • • • • • • • • • • d 08r-E E i c s 06r-0 n S a t y a t t O m n A intialiyatC.5nTTi 0 HSaByatOiTm 0 50 100 150 x position along channel (rrm) • I • i • ilp 50 100 150 Q8 0.8 I t 06 £ 0.6 h Q4 J 0.4 f-.0), (70,60), (70.80) (70,80) 02 02faoooooooori)flfB956rJt • • • • • • • • • • • • • • • • • (70,60) (70,0) n ' 1 1 1 1 ' 1 1 1 I . • : . I • 0 50 100 150 x position along channel (mm) 0.6 0.4 0.2 Figure 5.10 Centers of fibers at XZ plane with different initial height. Figure 5.11 Centers of fibers at three-dimensional space with different initial height and different angle 0. 50 100 150 t n a m«aigiairffiyrgtg|^rii'rffr1i'ry'iyrS'rvy|^"ti (20,60), (50,60), (70,60) 08r-E E I 0.6 .g04 0.2 I I I I I I i I 100000 O00O0O000OO1 (70,60) (80,60), (90,60) 0.8 0.6 0.4 0.2 50 100 150 x position along channel (mm) Figure 5.12 Centers of fibers at 3D with different initial height and different angle 0. 5h > re o 0 4 O i 1 3 c 8 2 c O o 1 _g 1mm fiber in 2D - 9 — 1mm fiber in 3D 1 2 Channel height z (mm) Figure 5.13 Probability of 1mm fibers moving in turbulent mean flow with 2D and 3D initial orientations. Chapter 5 Numerical simulation of fiber concentration 63 Figure 5.13 is the probability of fibers at each height interval in two-dimensional and three-dimensional space. 8000 fibers are randomly chosen between the wall and 2mm height with random orientation. We can see that in 2D, the probability is zero below the half fiber length; it is one above that height and it tends to be large (about 5.2) at that height. In 3D, the concentration is about 1 across the whole channel height (the isolations of the curve results from the initial distribution which is not exactly uniform). But the real flow is turbulent and both mean and fluctuation parts of the turbulence should be considered for the calculation of fiber concentration. 5.4 Fiber motion in turbulent flow with diffusion of fibers In order to record the fiber positions, several stations are set along the channel length. The channel height is divided into equal height intervals. The concentration is calculated by N , , where N M N, is the number of fibers at / th height interval and Nme is the average number of the fibers at each height interval. 8000 rigid fibers of 1mm, 2mm, 3mm are initially chosen across the channel with random height and random orientation angles 0 and $ such that the number of fibers is about the same at each height interval. It should be noted that the angle 6 could not be randomly chosen near the wall due to the fact that the fiber can not go through the wall. If a fiber initially goes through the wall with a random 9, we abandon this 6 and chose another 6 with the same height, so that it can not go in the wall. Note also that by choosing a random orientation and equal number of fibers at each height, we are assuming a perfectly mixed fiber/fluid mixture at the entry. This is probably incorrect for any real application, which will be strongly affected by upstream conditions. The perfectly mixed assumption is used here for simplicity and in the absence of better information. For a fixed length and diameter fiber, the number of the spheroids employed in the fiber can be relatively small to reduce the computation time, which is one of the benefits of modeling fibers as linked spheroids. The results for using different number of spheroids are slightly different due to the Chapter 5 Numerical simulation of fiber concentration 64 treatment of the wal l , but the global results remain unaffected (Examples are shown in section 5.5). For the current simulation, 2-4 spheroids are used for 1mm fiber simulation and 2-6 spheroids are used for 2 and 3 mm fiber simulation. 5 . 4 . 1 C o n c e n t r a t i o n f r o m R a n d o m W a l k The fibers are injected into the mean flow field after the flow field is converged by R A N S . The flow field and the fiber motion are calculated separately. Figure 5.14 presents the comparison o f concentrations for different fiber lengths. Except for the area very near to the wal l , the concentration increases, has a peak at a distance o f about one-half fiber length and then tends to be I, The distributions for different fiber lengths are different. The smaller the fiber length, the steeper the slope in the concentration curve. If the height is scaled by fiber length, the slopes collapse (Figure 5.15). This agrees reasonably well with Olson's fitted curve for his experiment (1996). Olson's fitted concentration curve is C/c J 3 ' 2 2 z , i ' z / i < " 3 ' 2 2 ' (5.16) m l l, z/Z->l/3.22. v ' > O O .o 2 E o 2 c c o o H ft 1 m m f i b e r w i t h 2 s p h e r i o d s — — — — 2 m m f i b e r w i t h 2 s p h e r o i d s 3 m m f i b e r w i t h 2 s p h e r o i d s 1 H-,!l »UVj h ' i ' U i.K 1 1 1 1 J i l t 1 AT m 5 10 15 C h a n n e l h e i g h t z ( m m ) 2 0 Figure 5.14 Concentration for 1mm, 2mm and 3mm fibers vs. channel height, using Random Walk. Chapter 5 Numerical simulation of fiber concentration 65 1mm fiberwith 2 spheriods 2mm fiberwith 2 spheroids 3mm fiberwith 2 spheroids Olson's fitting curve for his experiment data | i i i i i i i i i i i 2 4 6 C hannel height z IL Figure 5.15 Concentration for 1mm, 2mm, 3mm fibers vs. channel height scaled by fiber length, using Random Walk. 5 . 4 . 2 C o n c e n t r a t i o n f r o m L E S After the flow field from L E S is statistically steady, the fibers are put in. Then each L E S time step calculation is followed by one step of the fiber motion. One large difference between the flow field o f the Random Walk and the L E S is this: as long as the mean flow field is calculated by the Random Walk, it w i l l keep unchanged for the whole procedure o f the fiber concentration calculation. But for L E S , the flow field is updated each time step for the fiber concentration calculation. Figure 5.16 shows a comparison of concentration for 1mm, 2mm, 3mm fibers. W e see that the concentration curve o f 2mm and 3mm fibers are very close to Olson's fitted curve. The 1mm fiber does not have the same curve. The grid in the x direction (0.6mm) is probably too coarse for the 1mm fiber simulation. Chapter 5 Numerical simulation of fiber concentration 66 For Random Walk calculation, the flow field was calculated using the standard k — £ model for isotropic effective viscosity, which is not exact for the near wall region. There is no experimental data for the velocity fluctuations very close to the wall. This limits the application of the Random Walk for fiber motion in the channel and other equipment. LES will constitute a better technique to simulate the fiber-fluid interaction in the future, particularly as better and faster computer resources are developed. L E S : 2mm fiber — — — - L E S : 3mm fiber — LES: 1mm fiber Olson's fitting curve for his experiment data TO VW I I 2 3 4 Channel height z/L Figure 5.16 Concentration for 1mm, 2mm, 3mm fiber vs. channel height scaled by fiber length using LES. Chapter 5 Numerical simulation of fiber concentration 67 5.5 Discussion 5.5.1 E f f e c t s o f s o m e f a c t o r s o n t h e f i b e r c o n c e n t r a t i o n This section describes the effects of several factors on the concentration, calculated by the Random Walk method. Number of spheroids constituting the fiber: For the same length fiber, the number of spheroids that form the fiber does not affect the global result. Figure 5.17 is the concentration for 2mm fiber with different numbers of spheroids. More spheroids produce more accurate solutions, because a shorter radius of the spheroid makes the wall force in the wall model more accurate. The wall force is added to the center of the spheroid no matter which part of the spheroid touches the wall. But the computational time is greatly increased with more spheroids. For practical calculations, 3 or 4 spheroids are enough for 2-3 mm fibers. 2mm fiber with 4 spheroids -~~ — — _ 2mm fiber with 3 spheroids 2mm fiber with 2 spheroids Olson's fitting curve for his experiment data 2 O o 0 2 4 Channel height z/L 6 Figure 5.17 Concentration for 2mm fiber with different number of spheroids. Chapter 5 Numerical simulation of fiber concentration Number of fibers: In the above calculations, 8,000 fibers are used for each fiber length calculation. In order to verify that 8000 fibers are enough to get the needed information, 16,000 and 24,000 fibers are also tested for 3mm rigid fibers composed o f 2 spheroids. Figure 5.18 shows a comparison o f the concentration with 8000, 16,000 and 24,000 fibers with the error bars that stands for the standard deviation from the mean o f these three sets o f data. The result shows that the three sets o f data produce no significant difference along the whole channel width. More data results in smoother curve away from the wal l (2 mm < z < 18 mm), but no obvious effects on the slope near the two walls. So we use 8,000 fibers for the calculation. Channel height z(rn m) Figure 5.18 Concentration o f 8,000, 16,000 and 24,000 3mm rigid fibers. Chapter 5 Numerical simulation of fiber concentration 69 Stations: The concentration in the previous calculations can quickly attain its steady state. Figure 5.19 is the concentration for 3mm fibers at 100mm, 150mm and 198mm along the channel. It shows that the concentration curves at these three stations are very similar. — — — — 3mm fiber with 4 spheroids at station x=1 00mm _ 3mm fiber with 4 spheroids at station x=1 50mm > 3mm fiber with 4 spheroids at station x=198mm to O U 2 -Channel height Figure 5.19 Concentration at different stations. Fiber flexibility: Figure 5.20 is a comparison of concentration for 3mm rigid fibers and flexible fibers. It shows that the flexibility has not much influence on the fiber concentration in this simple channel flow although some fibers were observed to bend significantly. Chapter 5 Numerical simulation of fiber concentration 70 — 3mm rigid fiberwith 4 spheroids _ 3mm flexible fiberwith 4 spheroids to O O 2 -0 H. i i i I i i i i I i i i i I i i i i I i i i i I i i i i I i i i i I i i i i I i i i i I i i i ill 0 2 4 6 8 10 12 14 16 18 20 Channe l height (mm) Figure 5.20 Concentration for 3mm rigid and flexible fiber. 5.5.2 The orientation of fiber The angle o f fiber orientation relative to the channel wall and the height in xz plane and xy plane could easily be obtained during the concentration calculation. The probability density contours were generated by binning the fibers into 15 angle bins between -K12 and 7t 12 and 20 height bins between 0 and 2.0L. Figure 5.21, 5.22 and 5.23 show the contour plots o f the probability for fiber orientation and length for 1mm, 2mm and 3mm fibers in xz plane, respectively. Figure 5.24 is the comparison o f the average orientation in xz plane. Chapter 5 Numerical simulation of fiber concentration 7 1 P r o b a b i l i t y d e n s i t y c o n t o u r s in xz p l a n e with t u r b u l e n c e , L=1 -1.5 -1 -0.5 0 0.5 J L -0.5 0 0.5 Orientation [radians] m m \ 1.8 1.6 1 .4 1 .2 1 0.8 0.6 0.4 5 ° 2 Figure 5.21 Probability density contours for height z / L and orientation for 1mm fibers in xz plane. P r o b a b i l i t y d e n s i t y c o n t o u r s in xz p l a n e with t u r b u l e n c e , L=2mm 0.2 1 .5 -0.5 0 0.5 Orientation [radians] 1 .5 Figure 5.22 Probability density contours for height z / L and orientation for 2mm fibers in xz plane. Chapter 5 Numerical simulation of fiber concentration 72 r o b a b i l i t y d e n s i t y c o n t o u r s in xz p l a n e with t u r b u l e n c e , L = 3 m m •1.5 -1 -0.5 0 0.5 1 1.5 ^-.0.11 7 ; -0.5 0 0.5 Orientation [radians] Figure 5.23 Probability density contours for height z / L and orientation for 3mm fibers in xz plane. 2.5 -0.5 —i 1 .5 JZ ro 1 1 0.5 r -• 1 m m f i b e r O 2 m m f i b e r _ 3 m m f i b e r -0.5 0.5 • o • c» v V 5 v I V -57 r V V V V V -I V V V J I I 1 —1 I I l_ Orientation [radians] 2.5 H1 .5 H i H 0 . 5 0.5 Figure 5.24 The comparison o f average orientation in xz plane for 1mm, 2mm and 3mm fibers. Chapter 5 Numerical simulation of fiber concentration 73 Those figures show that the fibers are preferentially positively oriented in xz plane and as fiber length increases, the preferred orientation becomes much more distinct, which is similar to Olson's observation (1996). When the fiber at a height below half of the fiber length, it will rotate and move to a higher position with a slightly positive orientation because of the impact with the wall. At this new height, the fiber rotates slowly before being transported away by turbulent fluctuations, which results in the preferred positive fiber orientation (Olson 1996). 5.6 Summary Fiber motion and concentration in the turbulent channel flow have been studied. The flow distribution was computed first, in one of two different ways: either the Reynolds Averaged Navier-Stokes equations were solved with the k - e model of turbulence or a Large Eddy Simulation (LES) was used. The distribution of fiber concentration was then obtained using the computed flow in two corresponding ways: a Random Walk method was developed to use the k - e predictions, and a direct calculation of fiber motion was made, based on the detailed flow field computed by the LES method. The Random Walk method uses experimental RMS distributions for the three components of velocity close to the wall, as these components are not predicted separately by the k - e turbulence model. From the Random Walk method, with fibers of 1mm, 2mm, 3mm lengths, the concentration results scaled by the fiber length show a linear region near the wall and an almost constant region above a height of about half fiber length. Similar results were observed in the experiments of Olson (1996). A peak between the linear and constant concentration regions is observed in the simulation, a feature which also appears (with somewhat lowered magnitude) in Olson's (1996) data. Using the LES method with different length fibers, the concentration curves of the 2mm and 3mm fibers are similar to the experiments. The 1mm fiber results do not fit well with other computed or measured results, probably because the grid size for the LES calculation in the x direction Chapter 5 Numerical simulation of fiber concentration 74 (0.6mm) is too coarse to model the turbulence adequately at the scale of the smallest fibers. A finer grid could not be used for the current LES calculation because of computer memory limitations. The LES method is more general, as it does not rely on experimental data and gives a better representation of the turbulence characteristics, but it is computationally much more expensive than RANS method. As the LES method is made more general and as computer speed and memory increase, the LES method will be preferable. Currently the RANS method is still the practical computational procedure for turbulent flow in complex geometries. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 6.1 Introduction Several studies have been conducted on fiber passage through an aperture. Gooding (1986), Gooding and Kerekes (1989) filmed the trajectories of 3mm long rigid fibers through a smooth slot. Yu and Defoe (1994b) observed the paths of a small number of fibers approaching a few different screen slots. Tangsaghasaksri (1994) numerically calculated the fiber trajectories before the fibers contact the wall and thus did not give a whole picture of fiber trajectories. Olson (1996) developed a two-dimensional theoretical model to simulate the probability of passage through a smooth slot. His work will be discussed in some detail later in this thesis. Lawryshyn (1997) used CFD to examine the passage into a smooth slot of fibers with only one initial orientation (assuming the fibers are all along the flow). This thesis used the three-dimensional simulation tool described in chapters 3 and 4 to examine the passage through several slot shapes, using fibers of different length, with random initial orientation and placed initially at a range of initial heights. Section 6.2 examines the flow field through different slot shapes, and provides results including the velocity vector, hydraulic resistance and turbulence intensity. Section 6.3 shows the "fate plot" of the fiber trajectories as a function of initial height and orientation. Section 6.4 analyzes the stapling phenomena. Section 6.5 explores fiber passage for several different conditions and slot shapes. Section 6.6 is the summary of this chapter. 6.2 Flow field simulation This section examines the flow field for different shapes of screen slots. As explained in chapter 3, the flow in all the slots is calculated using the RANS model. 75 Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 7 6 6.2.1 Slot geometry and boundary conditions The flow field consists of three parts as shown in Figure 6.1 for a smooth slot. In this figure, Part 1 is the channel. For all the calculations, the channel height H is 20mm; this is based on the distance between the rotor and the screen plate in an industrial screen as shown in Figure 1.1. The channel length / is about 400mm. The upstream length /„ , which is the distance between the inlet of the calculation domain and the middle of the slot, is 300mm and it is set long enough so that the entry flow attains a fully developed profile. Part 2 is the slot with width W = 0.5 mm for all the calculations; Slot height is d = 1 mm. Part 3 is a small chamber with width 3mm and length 60mm to ensure the flow is fully developed at the outlet 1. The boundary conditions are as follows: Inlet: uniform velocity 6.5 m/s. As already noted, the entry length is long enough so that this flow becomes essentially fully developed. The turbulent kinetic energy and the dissipation rate of turbulent kinetic energy at entry are specified same as (5.1) for the channel flow calculation Outlet 1: "gradient condition", which means the normal derivatives of all the variables are set equal to 0, and a fixed mass flow rate is specified. Outlet2: "gradient condition". Walls: no-slip wall conditions. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 77 inlet outlet 2 Part 1 W Part 2 Part 3 outlet 1 Figure 6.1 Geometry o f a smooth slot. The contoured slots we have examined include the step-step contoured slots and the slope-slope contoured slots as shown in Figure 6.2. The step-step contour is a Lehman contour where the upstream and downstream sides are normal to the screen plate. In the slope-slope contour, the upstream and downstream sides are replaced with two sloped sides. For the step-step contour, three heights are examined: 0.1mm (CO), 0.5mm ( C l ) , 1.2mm (C2). For the slope-slope contour, two heights are examined: 0.5mm ( L l ) , 1.3mm (L2), which are as about the same height as C l , C2 respectively. L 2 is a C A E profile which was provided by Wherrett in C A E (1999). For all the calculations, the upstream velocity is 6.5 m/s and, to examine the effects o f slot velocities on fiber passage, six slot velocities are used: 1.5 m/s, 2.4 m/s, 3.3 m/s, A Amis, 5.4 m/s, 6.5 m/s, as used by Lawryshyn (1997). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 7 8 Wc=3.0 dc--> = 0 . 1 , 0 . 5 , 1 . 2 e => w Step-step contour slot geometry: CO, C l , C2 <= n k V W=3.5 ^ ^ = 0 . 5 , 1 . 3 n e r> w Slope-slope contour slot geometry: L l , L2 Figure 6.2 Slot geometry and dimensions of step-step and slope-slope contours (all dimensions are in mm). 6 . 2 . 2 G r i d i n d e p e n d e n c e Stretched grids are created to most suitably calculate the details of the flow patterns in every part of the geometry. Fine grids are used for the area where there are large flow variations. The simulations of the flow field should be relatively independent of the grid size. For example, for the smooth slot calculation, three sets of gradually finer grids (Table 6.1) are tested. The grids in set "Grid 2" is 1.2 times finer in the x and z directions than that in "Grid 1". Grids in set "Grid 3" is 1.2 times coarser in x and z directions than that in "Grid 1". Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 79 Table 6.1 Grids used for testing grid independence for smooth slot. Grids Parti cells (xxz) Part2 cells (xXz) Part3 cells (xxz) Total cells 3 168 x96 6 xlO 22 x42 136896 1 208 x 116 8x12 28 X 5 0 204992 2 250 x140 10 X 1 4 34 x60 297440 Four locations are selected on which to compare the calculated velocities (Figure 6.3). Position 1: -2<x<2 , z = 0.1; Position2: -2 < x <2, z = 0.3. Position 3: -2 < x < 2, z = 0.5; Position 4: x = 0, - 0.5 < z < 0.5 . T O . ' l " i— 0.3 0.5 p o s i t i o n 3 p o s i t i o n 2 p o s i t i o n 1 fc.25 50.25 0.5 Figure 6.3 Two positions for comparison o f velocities with different meshes (all dimensions are in mm). Figure 6.4 is a comparison o f the mean u velocities at position 1, 2 and 3. It is shown that the velocities from the finer " G r i d 2" is closer to the results o f " G r i d 1" than that from the coarser " G r i d 3". Figure 6.5 is a comparison o f the mean w velocities at position 1, 2 and 3. N o obvious Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 80 differences are observed. Figure 6.6, 6.7 are a comparison of velocity u and w at position 4 along the slot centerline. The above plots show that the results o f " G r i d 1" do not change significantly from that o f " G r i d 2". The maximum error in the results from " G r i d 1" and " G r i d 2" is about 5%. " G r i d 1" is therefore used for the simulation o f the smooth slot. Figure 6.8 is the grid near the smooth slot. Similar grids are used for the contour slot simulations. -0.002 -0.001 0 0.001 x position (m m ) 0.002 Figure 6.4 u velocity comparison above the slot (position 1, 2 and 3) between different grids. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes Figure 6.5 w velocity comparison above the slot (position 1,2 and 3) between different grids z position (mm) Figure 6.6 u velocity comparison along the slot center (position 4) between different grids. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 82 o 1 1 II — — g r i d l -gr id2 (f iner) ^ — _ g rid 3 ( c o a rse r) ^ -^^ "^ --_i i 1 i i i i 1_ I I I I I I I J I I I I I I I L_ -0.0004 -0.0002 0 0.0002 0.0004 z position (m m ) Figure 6.7 w velocity comparison along the slot center (position 4) between different grids. Figure 6.8 Grids for smooth slot. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 83 6.2.3 Convergence The convergence criterion is that the residuals of all the variables are enough small. Figure 6.9 and 6.10 are the convergence histories for the six variables in Part 1 and 2, separately. These six variables are mass, u velocity, v velocity, w velocity, turbulent kinetic energy K and the dissipation rate of turbulent kinetic energy £. It can be seen that the residuals of all the variables drop over 12 orders of magnitude over the first 10,000 iterations and then keep constant values. As usual, £ is the variable that is hardest to converge. Figure 6.9 Residuals of six variables for Part 1 of smooth slot. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 84 Ite rations Figure 6.10 Residuals of six variables for Part 2 of smooth slot. 6 . 2 . 4 R e s u l t s 6.2.4.1 Flow field The smooth slot is the simplest case of all the screen slots. Figure 6.11 shows the streamlines in the flow near the smooth slot with different slot velocities. It is clear that a vortex appears on the upstream side of the slot no matter how big the slot velocity is. The larger the slot velocity, the smaller the size of the vortex. Figure 6.12 shows the streamlines for the Cl slot for the same six slot velocities as used in Figure 6.11. For all the slot velocities, there is no vortex in the slot and the flow in it is very smooth. Upstream of the contoured slot, there is an obvious vortex whose size changes with the slot velocity, larger for the lower slot velocity and smaller for the higher slot velocity. There is also a small vortex in the corner downstream of the contour whose size does not change much with the variation of the slot velocities. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes Figure 6.11 Streamlines in the smooth slot for different slot velocities (X, Z dimensions: mm) Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 86 -0.0015 -0.001 -0.0005 0 0 0005 0.001 0.0015 -0.0015 -0.001 -0.0005 0 0.0005 0.001 0.0015 X X Figure 6.12 Streamlines in Cl slot for different slot velocities (X, Z dimensions: mm). Figure 6.13 shows the streamlines in slot L2. The size of the vortex in the upstream side of the contour decreases with an increase of the slot velocity. The size of the relatively small vortex in Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 8 7 the downstream region of the slot also decreases with the higher slot velocity and the vortex disappears from the slot for a slot velocity of 5.4 m/s. Figure 6.13 Streamlines in L2 slot for different slot velocities (X, Z dimensions: mm). In order to observe the effects of changing the height of the step-step contour and slope-slope contour, CO with the shallower height (0.1mm) and C2 with deeper height (1.2mm) are investigated to compare with Cl (0.5mm). The shallower height Ll (0.5mm) is compared to L2 (1.3mm). In Figures 6.14, (a), (b), (c) show streamlines in CO, Cl, C2 for the low slot velocity (2.4 m/s) and (d), (e), (f) are streamlines in CO, Cl, C2 for the high slot velocity (5.4 m/s). It is shown that the flow field differs considerably with changing height for the step-step contour. The flow features in CO slot are closer to that in the smooth slot. There is a small vortex in the upstream corner and a relatively bigger vortex in the slot. The vortex in Cl is much smaller than that in C2 and the Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 88 flow approaching the slot is smoother in C 1 . A small vortex appears in the downstream side o f the slot in C 2 for the low slot velocity and disappears for high slot velocities. -0-0015 -0.001 -0.0005 0 0.0005 0.001 0 0015 -0 0015 -0 001 -0 0005 0 0 0005 0.001 0.0015 X X Figure 6.14 F low field for step-step contour with different contour height ( X , Z : mm). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 89 Figure 6.15 shows the streamlines in slope-slope slot LI, L2 with low slot velocity (2.4 m/s) (a, b) and high slot velocity (5.4 m/s) (c, d). The streamlines, the size and the shape of the vortex are also quite different for the slope-slope contour with different height. The above results show that contour height is one critical factor in the design of the contour slot. -0 001 0 0 001 0 002 0 003 -0.001 0 0.001 0.002 0.003 -0.001 0 0 001 0 002 0 003 -0.001 0 0.001 0.002 0.003 Figure 6.15 Flow field for slope-slope contour with different contour height. For contours of the same height 0.5mm, Figure 6.16 shows four different shapes. Clsmooth has the same upstream side as that of Cl slot and the same downstream side as that of the smooth slot. Cl slope also has the same upstream side as that of Cl but a similar downstream side as that of LI. Cltwoslope is similar to LI except that the slope angle is different. Cl, Cl_slope and Cltwoslope have the contour size 3.0mm, being the distance between the upstream edge and Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 90 downstream edge o f the contour. Figure 6.17 shows a comparison o f streamlines o f C l smooth, C l_s lope and C l t w o s l o p e for two slot velocities: 2.4 m/s and 5.4 m/s. It appears that i f the contour has the same upstream side (e.g. a, b, d and e), no matter what the downstream side looks like, the size o f the upstream vortex is about the same for the same slot velocity and decreases with an increase in the slot velocity. If the contour has a different upstream size, even i f the contour has the same downstream side (e.g. b, c, e and f), the size o f the vortex is much different. This difference is more obvious for the high slot velocity (e.g. the vortex disappear in (f)). H o w these flow fields affect the fiber behavior w i l l be discussed in section 6.3. 3.0 m m C l slot C l smooth C l slope C l _ t w o _ s l o p e Figure 6.16 Geometries for contours with the same height (0.5mm) and different shapes. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 91 Figure 6.17 F low field for contours with the same height (0.5mm) and different shapes ( X , Z dimensions: mm). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 92 6.2.4.2 Hydraulic resistance The hydraulic resistance of a screen plate has an important effect on the screen performance. Thomas and Cornelius (1982) experimentally studied the laminar flow resistance in pulp screens. The flow resistance was assessed using a non-dimensional pressure drop coefficient: K = - ^ (6.1) where the pressure drop Ap was the difference between the pressure in the mainstream flow above the slot and that in the plenum below the slot. They observed the existence of a vortex on the upstream side of the slot wall, which is similar to the vortex in Figure 6.11. They also noted that, with an increase of upstream velocity or a decrease of slot velocity, the size of the vortex and the pressure drop coefficient increased. Gooding (1996) studied the hydraulic resistance using three methods: computational fluid dynamics (CFD), flow channel experiments with a highly simplified screen configuration, and pilot plant trials with an industrial pulp screen. A one-dimensional energy equation was applied to the two control volumes shown in Figure 6.18. Control volume 1 was used for his CFD calculation and control volume 2 was used for the flow channel experiment and was also used for the current CFD calculation. C o n t r o l v o l u m e I C o n t r o l v o l u m e 2 Figure 6.18 Two control volumes used by Gooding (1996). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 93 The one-dimensional energy equation yields: pl + a,(\pV2) + pgZx = p2+a2(\pV22) + pgZ2 +Ap (6.2) where /?,, V, andZ, represent the average pressure, velocity and elevation at plane 1 and p2, V2 and Z2 are the average pressure, velocity and elevation at plane 2. Planes 1 and 2 were chosen at places where the flow is fully-developed. or, and GC2 are kinetic energy correction factors to account for the cross-stream velocity gradients in the flow coming into and out of the control volume (Roberson and Crowe 1985). For the flow going into the slot, Gooding (1996) assumed the power-law turbulent velocity profile near the screen plate and derived the kinetic energy correction factor or,: or, (\pV2) = 0.520 r pV}'AWu*Vj'4>\ 1/4 (6.3) where yc is half of the channel height. Assuming the flow at plane 2 is turbulent and fully-developed, the kinetic energy correction factor a2 was expressed as (Roberson and Crowe 1985): V a2(^pV2) = 1.05 \/2p V A2 J (6.4) The total pressure loss is thus expressed through: Ap = P x - p2 + or, (1 /2pV2) -a2(1 /2pV2) + pg(Zl - Z2) (6.5) and the total pressure drop coefficient is then obtained by (6.1). The total pressure loss may include the upstream pressure drop, pressure drop at the slot entry, pressure drop along the slot depth and through the plenum. Figure 6.19 compares the total pressure drop obtained from the current CFD calculations (filled circle symbols) for an upstream velocity of 6.5 m/s with the slot entry pressure drop from Gooding's CFD results (square, gradient symbols stand for pressure drop for upstream Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 94 velocity 7.5 m/s, 5 m/s separately). The current pressure drop for a smooth slot lies between the two lines of Gooding's results, which shows the approximate agreement between the two different simulation methods. The entry pressure drop is the major part of the total pressure drop and other losses are relatively small. Figure 6.20 includes the pressure drop coefficients associated with the pressure drop in Figure 6.19. The current pressure drop coefficients match reasonably well with Gooding's results. The trend in both methods is that the hydraulic resistance decreases with increasing slot velocities, which is similar to Thomas and Cornelius (1982)'s observations. Figure 6.21 compares the pressure loss with the smooth slot and the contour slots CO, Cl , C2. Figure 6.22 shows the pressure drop coefficient for the smooth slot, CO slot, Cl slot and C2 slot. The pressure drop and hydraulic resistance for the CO slot is about the same as that of the smooth slot. The pressure drop for the Cl slot is greatly reduced compared to that of the smooth slot. The hydraulic resistance of the C2 slot is almost the same as that of the smooth slot and CO slot at slot velocities below 3 m/s, and is reduced considerably for high slot velocities. The drop is smaller than that for Cl slot. The above results show that the contour slots generally decrease the hydraulic resistance, but the contour height may be the most important variable in the design of the contour slots. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes Figure 6.19 Pressure loss for the smooth slot compared with Gooding's results. 20 *X 15 £ 9 O 5 a. 10 2 -a £ 3 (fl (fl e 5 • smooth slot: current calculation, Vu = 6.5 m/s • smooth slot: Gooding's calculation, Vu = 7.5 m/s \7 smooth slot: Gooding's calculation, Vu = 5 m/s V ••7'n, 0 I ' 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 , 1 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Normalized slot velocity, Vs/Vu Figure 6.20 Pressure drop coefficient for smooth slot compared with Gooding's results. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 70 60 (0 0_ T 40 w o 8? 5 30 n </> 22 Q- 20 10 r --e-- B -Current calculation: Smooth slot C 1 sot C 2 slot _ _<*>_ _ CO slot 0 *~ ' ' 1 1 1 1 1 1 1 1 1 ' 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 2 3 4 5 6 7 8 Slot velocity Vs (m/s) Figure 6.21 Pressure drop for smooth slot, C l slot and C2 slot. 1 5 c o 10 £ o u C L 2 •D w in CD 5 h Current calculation: Smooth slot C 1 slot C 2 slot CO X X 0.2 0.4 0.6 0.8 1 Normalized slot velocity Vs/Vu 1.2 Figure 6.22 Pressure drop coefficient for smooth slot, C1 slot and C2 slot. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 97 6.2.4.3 Turbulence intensity Turbulence intensity is defined as (Gooding 1996): 1/2*7 (6.6) where K is the turbulence kinetic energy. Gooding (1996) used / , to identify the contributions to pressure drop coefficient K for a smooth slot and a step-step contour slot. This section first shows the turbulence intensity in the smooth slot with low and high slot velocity. The present results are similar to those o f Gooding (1996). Most importantly, this section examines the turbulence intensity for contours having the same open area and different height (CO, C l and C2) and contours having the same height and different open area ( C l s m o o t h , C l s l o p e , Cl_two_slope and L I ) . This could identify how a single factor like the contour height or open area affects the flow structure. Figure 6.23 and Figure 6.24 shows the turbulence intensity for the smooth slot with slot velocity 2.4 m/s and 5.4 m/s separately. It is clear that the turbulence intensity for the low slot velocity is much bigger than that for high slot velocity so that K increases with I . The trend is also observed by Gooding (1996). 0.001 0 . 0 0 0 5 0 . 0 0 0 5 I— -0 .0 01 I— - 0 .00 1 5 - 0 . 00 1 5 -0 .00 1 0 . 0 0 0 5 0.00 1 0.00 1 5 Figure 6.23 Turbulence intensity in smooth slot for low slot velocity (Vs = 2.4 m/s). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 9 8 - 0 . 00 1 5 -0 .00 1 - 0 . 0 0 0 5 0 0 . 0 0 0 5 0 .00 1 0.00 1 5 X Figure 6.24 Turbulence intensity in smooth slot for high slot velocity (Vs =5 .4 m/s). Figure 6.25 includes the turbulence intensity for smooth, CO, C l and C2 slot with slot velocity 5.4 m/s. For the CO slot, the flow structure is similar to that o f smooth slot, except that the turbulence intensity is bigger in downstream corner o f the slot and the plane just below the vortex. There is no obvious effect o f the downstream corner o f the contour in CO slot. The larger the contour height, the bigger the effects o f the downstream corner ( C l , C2). It is interesting that the total turbulence intensity in the slot C l is less than that in CO and smooth slot. But for the C 2 slot, the turbulence intensity is much bigger both in the slot and in the downstream part o f the contour. Figure 6.26 displays the turbulence intensity for the slot contours with same height and different open area, al l for the same slot velocity o f 5.4 m/s. The turbulence intensity in the C l_smooth downstream slot corner is less that that in other contours and is even less than that in the smooth slot with same slot velocity. C I s l o p e and C l t w o s l o p e have the same shape o f downstream contour. Turbulence intensity in the C l slope case is about the same as that in the C l j w o s l o p e . It seems that the upstream contour for the same open area does not produce Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 99 Figure 6.25 Contour o f turbulence intensity for smooth, CO, C l , C 2 slots. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 100 0.0011 0.00051 N -0.0005 -0.001 -0.0015 -0.0015 -0.001 -0.0005 I Level I 6 0.50 5 0.40 4 0.30 3 0.20 2 0.10 1 0.00 0 00005 0.001 00015 X 0.001 0.0005 N '-0.0005 h -0.001 -0.0015 r--0.0015 -0.001 -0.0005 0 X 1 1 1 1 1 1 ' 1 1 1 1 1 ' 0.0005 0.001 00015 0.001 0.0005 N '-O.Cttfir--cool h -0.0015 h 0.001 0.00051 N '-0.0005 -0.001 h -0.0015 h Level I • 6 0.50 L1 M 5 0.40 U 4 0.30 1 3 0.20 1 2 0.10 | 1 0.00 -0.0015 -0.001 -0.0005 0 0.0005 0.001 O0015 X 0.001 0.002 0.003 Figure 6.26 Contour o f turbulence intensity for C I s m o o t h , Cl_s lope , C l t w o slope and L l slots. 6.3 Fiber motion in different slot shapes This section examines how the different flow fields affect the trajectories o f fibers with different length and flexibility. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 101 As explained in Chapter 4, the fiber position in three-dimensional space is characterized by three factors: the center position, angle 8 (which is the angle between the fiber main axis and the Z-axis), and angle <f> (which is the angle between the projection of the fiber main axis in XY pane and Y-axis). For each different flow field, different number of fibers were tested, and 4000 fibers are chosen with fixed center positions in the Y direction (5mm) and X direction (-2mm for 1mm fiber and -3mm for 3mm fiber, chosen enough close to the slot to examine only the effects of the entry flow on fiber trajectories), random height from the wall to 1mm above the wall, and random 6 and <p from 0°to 180°. For the case when the height above the wall is below the half fiber length, 6 cannot be chosen for the whole range. For fixed height, if one choice of 6 is impossible (the fiber goes into the wall), another 6 is randomly chosen to make sure that it is above the wall initially. Each fiber has three possible fates: it passes through the aperture, continues down the channel or staples at the contour corner or the slot corner. It must be pointed out that the stapling of fibers is only crudely simulated by this numerical work since stapled fibers are least likely to be well represented by the Stokes flow assumption and are also likely to be affected by neighboring fibers in actual conditions. All the figures in this section are "fate plots" (as introduced by Olson (1996) for the smooth slot) showing the fiber trajectory as a function of initial height and orientations 6 and <P. The symbols used in the plots are in Table 6.2. The pink solid line is the exit layer height for each slot velocity. Table 6.2 Symbols used in all the "fate plots". Symbol Meaning Green symbol"+" Passes through the slot Blue symbol "o" Goes away from the channel Red symbol "n" Stapled Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 102 6.3.1 Trajectories of 1mm rigid fiber (L/W = 2) Figure 6.27 and Figure 6.28 are fate plots for 1mm rigid fibers in the smooth slot with six different slot velocities as a function of initial height and orientation 0, (j) respectively. It is shown that more fibers go into the slot having bigger slot velocities. Fibers with initial orientation 0 > 7t 12 have a greater tendency to go into the aperture than for other 0 orientations. If the fiber is initially above the exit layer, it goes into the aperture only if (/> is close to 7112 , which means that the fiber is initially in or close to the XZ plane. In two-dimensional calculations of a smooth slot, Olson (1996) observed that the fibers with initial negative orientation (0Q < 0) could be drawn from a region beyond the exit layer, where 0O is the angle between the fiber and X-axis. Olson's observation is a special case of the current three-dimension observations. If the fiber is initially within the exit layer, it is seen that the fiber tends to pass through the aperture when its initial angle (j> is close to 0 or it/2 , which means that the fiber is approximately perpendicular to the flow. This three-dimensional phenomenon could only be observed in three-dimensional calculations such as those presented here. Some fibers staple in the region between where they pass into the aperture and where they exit from the channel. This will be discussed in the next section. Figure 6.29 and Figure 6.30 are fate plots for the 1mm rigid fiber in the Cl slot with six different slot velocities. The difference of the fiber trajectory between the smooth slot and the Cl slot cases is that for Cl slot, almost all fibers initially within the exit layer will go into the aperture no matter what their initial orientation is. For fibers above the exit layer, the fibers with 0 > n 12 also have more chance to pass through the aperture. From Figure 6.30, the fibers are stapled near the exit layer and the three-dimensional phenomenon is not obvious. Figure 6.31 gives comparisons of the fate plot of 1mm rigid fibers in CO, Cl and C2 slots with the same slot velocity 5.4 m/s. It is seen that the fate plot in CO is very close to that in the smooth slot. The obvious effect of initial 0 > n 12 inClis not observed in C2. Only the exit layer Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 103 height that separates the fibers which go through from those which do not go through and staples some fibers is important. Figure 6.32 is the similar set of fate plots for the slope-slope contour LI and L2. It is obvious that the number of stapled fibers is greatly reduced. There is no fiber stapled in the L2 contour. More fibers above the exit layer in LI go through the aperture than that in L2. Figure 6.33 are the similar fate plots for Cl_smoofh, Clslope and Cl_two_slope contours. Almost no fibers staple near the exit layer. The trajectory in Cl slope and Cl_two_slope is almost the same, which can be predicted from the similar downstream flow field and turbulence intensity shown in Figure 6.17 and Figure 6.26. The Clslope and Cl_two_slope allow more fibers to go into the aperture than the Cl smooth slot, which means that the downstream contour is important in design. The above fate plots show that the initial height and orientation have different effects on the fiber trajectories for different slot shapes. Overall, the assumption that fibers in the exit layer enter the slot is a reasonable approximation, but slot shape is also important. Generally, contour slots let more fibers go into the aperture than the smooth slot. For contour slots, more fibers staple in the step-step contour than in the slope-slope contour. For the same shape contour, contour height is an important factor in determining the fiber trajectory. For the contour of the same height, the angle of downstream side of the contour is also important. These findings could provide useful guidance for future contour designs. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 104 cfd! Vt-1.5 m/s Cfd2 V» 2 .4 m/s 1.5 2 2.5 Initial angle e (radians) 1.5 2 2.5 Initial angle 6 (radians) cfd3 Vs-3.3 m/s Cfd4 Vs -4 .4 m/s 1.5 2 2.5 Initial angle e (radians) 1.5 2 2.5 Initial angle 8 (radians) cfdS Vs-5.4 m/s 1.5 2 2.5 Initial angle e (radians) 1.5 2 2.5 Initial angle e (radians) Figure 6.27 Smooth slot: The effects o f initial height and orientation 0 on passage o f 1mm rigid fiber for different slot velocities (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "a": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 105 Cfd1 V s = 1 .5 m/s 1.5 2 2.5 Initial angle 41 (radians) Initial angle $ (radians) ctd3 Vs-3.3 m/s cfd4 Vs-4.4 m/s 1.5 2 2.5 Initial angle 0 (radians) 1.5 2 2.5 Initial angle 4 (radians) cfd5 Vs-5 .4 m/s cfd6 Vs-6.5 m/s 1.5 2 2.5 Initial angle <!• (radians) 1.5 2 2.5 Initial angle 0 (radians) Figure 6.28 Smooth slot: The effects of initial height and orientation <p on passage of 1mm rigid fiber for different slot velocities (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "o": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 106 cfd2 V s - 2 . 4 m/s 1.5 2 2.5 3 Initial angle 8 (radians) 1.5 2 2.5 Initial angle 0 (radians) Cfd3 V s - 3 . 3 m/s cfd4 V s - 4 . 4 m/s 1 1.5 2 2.5 3 Initial angle 6 (radians) 1.5 2 2.5 Initial angle 6 (radians) cfdS V s » 5 . 4 m/s ctd6 V s - 6 . 5 m/s £ E — 0.6 £ • f 0.5 •C « 0.4 ±s = 0.3 0.2 0.1 1.5 2 2.5 Initial angle e (radians) 1.5 2 2.5 Initial angle e (radians) Figure 6.29 C l slot: The effects of initial height and orientation 0 on passage of 1mm rigid fiber for different slot velocities (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "a": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 107 c f d l Vs=1.5 m/s cfd2 V s = 2 4 m/s 1.5 2 2.5 initial angle o (radians) 1 1.5 2 2.5 Initial angle o (radians) cfd4 V s - 4 . 4 m/s 1 1.5 2 2.5 Initial angle « (radians) 1.5 2 2.5 Initial angle $ (radians) cfd5 V s - 5 . 4 m/s Cfd6 V s - 6 . 5 m/s 1.5 2 2.5 Initial angle $ (radians) 1.5 2 2.5 Initial angle 0 (radians) Figure 6.30 C l slot: The effects of initial height and orientation (j) on passage of 1mm rigid fiber for different slot velocities (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "o": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 108 CO: cfd5 Vs-5.4 m/s CO: C f d 5 Vs-5 .4 m/s 1.5 2 2.5 3 Initial angle e (radians) 2 2.5 Initial angle $ (radians) C1: cfdS Vs-5.4 m/s C1: cfd5 Vs-5 .4 m/s Initial angle e (radians) C2: CfdS Vs-5 .4 m/s C2: cfdS Vs>5.4 m/s * 0.4 as e ~ 0.3 0.2 r 0.1 -1.5 2 2.5 Initial angle e (radians) Figure 6.31 CO, C l , C2 slot: The effects o f initial height and orientation 6, (j) on passage o f l m m rigid fiber for slot velocity Vs = 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "o": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 109 L1: cfcfi Vs=5.4 rrfs .c e 0.5 0.4 03J 0.2 -0.1 -1 1.5 2 25 3 Initial ancje 8 (radans) L1: cfcfi Vs=5.4 rrfs 15 2 25 3 Initial antje oi (radans) LZ- cfcfi Vs=5.4 rrfs L2- cfcfi Vs=5.4 rrfs Injrjal andje 6 (radans) Initial andje 0) (radans) Figure 6.32 L l , L2 slot: The effects of initial height and orientation 0, <p on passage of 1mm rigid fiber for slot velocity Vs = 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "a": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 110 C1-smooth: cfdS Vs-5.4 m/> C1-smooth: cfd5 Vs=5.4 m/s 1 1.5 2 2.5 Initial angle e (radians) C1-slops: cfd5 Vs-5.4 m/s Cl-slopa: cfd5 Vs-5.4 m/s £ 0 . 5 a 0 .4 as c — 0.3 0 . 2 0.1 1.5 2 2.5 3 Initial angle e (radians) C1-two-slopa: cfd5 Vs-5.4 m/s C1-two-sk>pa: Cfd5 Vs-5.4 m/s t o s t A 1.5 2 2.5 3 Initial angle 9 (radians) 1.5 2 2.5 Initial angle o (radians) Figure 6.33 Clsmooth, Clslope, Cltwoslope slot: The effects of initial height and orientation 0, </> on passage of 1mm rigid fiber for slot velocity Vs = 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "a": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 111 6 . 3 . 2 T r a j e c t o r i e s o f 3 m m r i g i d f i b e r (L/W = 6) Figure 6.34 is a series of fate plots for 3mm rigid fibers in the smooth slot with the same two slot velocities as that for 1mm rigid fiber. Again, more fibers go into the slots with higher slot velocities. The tendency that fibers with initial 9>KI2 go into the slots is even greater than that for 1mm rigid fibers in the smooth slot. Fibers with initial <j> close to 0 or n are more likely to go into the slot, which means that the fibers in or close to the vertical plane go into the slots more easily. It appears that long rigid fibers have a stronger "three dimensional phenomenon" than short rigid fibers in the smooth slot. Figure 6.35 shows similar fate plots for the Cl slot. More fibers go into the slots compared to the smooth slot. The fibers in or close to the main plane (XZ plane) initially below the exit layer tend to be stapled. The comment about the "three dimensional phenomenon" still holds. This phenomenon can also be seen in Figure 6.36 (fate plots for C2 slot) and Figure 6.37 (fate plots for L2 slot). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 112 Smooth siat: dd2 Vs=Z4 rrfs Smooth sk± dd2 VSF2.4 rrfs 1.3 1.5 1.6 1.7 1.8 hitE. age 0 (radans) 1.9 Smooth siot cfd5 Vs=5.4 rrfs 1.1 1.2 13 1.4 1.5 1.6 1.7 1.8 1.9 intaande 9 (radans) 1.5 2 25 3 Htiaarge c[) (radans) Smoothskt drJ5 VSF=5.4 mfs Q3 Q2 Q1 15 2 25 3 Irirjaande o (radarE) Figure 6.34 Smooth Slot: The effects of initial height and orientation 9 , <f> on passage of 3mm rigid fiber for slot velocities 2.4 m/s, 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "o": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 113 Cl slot drJ2 Vs=2.4 mfs ci slot cfcE Vs=24 nYs mBaiangk} 0 (radais) Hbiangje | (ratfans) Figure 6.35 C l Slot: The effects of initial height and orientation 6, 0 on passage of 3mm rigid fiber for slot velocities 2.4 m/s and 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "o": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 114 C2statcfd2 Vs=24 rrfs C2sfctdd2 Vs=24 mfs 1 2 1.3 1.4 1.5 1.6 1.7 13 1.9 Irifeiande 9 (radars) C2sk±cfd5 VSr=5.4 mfs 0.9 08 Fi 07 ? 5. 06 £ •5 05 2 04 0.3 021-0.1 . » > • V i 1.1 12 13 1.4 1.5 1.6 1.7 1.8 1.9 2 21 Inttaarge 9 (radans) 15 2 25 3 Mbaanrje o (radans) 15 2 25 3 Mbaande <p (radans) Figure 6.36 C2 Slot: The effects o f initial height and orientation9, <j> on passage o f 3mm rigid fiber for slot velocities 2.4 m/s and 5.4 m/s (symbols same as Figure 6.35). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 115 L2 s lo t cfcE Vs=2.4 mfs I 2 s k ± c f d 2 V s=Z4 mfs 12 1.3 1.4 1.5 1.6 1.7 1.8 1.9 H U anQje 9 (radars) L2sSctcfrJ5 Vs=5.4 rrfs 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 irifeiangk) 9 (radars) 15 2 25 3 Hfcaange (!) (radars) 15 2 25 3 Irtiaian^e (j) (radars) Figure 6.37 Slope-slope L 2 Slot: The effects o f initial height and orientation 0 , <p on passage o f 3mm rigid fiber for slot velocities 2.4 m/s and 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "o": Stapled). 6.3.3 T r a j e c t o r i e s o f 1 m m flexible fiber (L/W = 2 ) Figure 6.38 shows similar fate plots of 1mm flexible fiber for smooth slot with slot velocities 2.4 m/s and 5.4 m/s. Compared to Figure 6.27 and 6.28 for fate plots o f 1mm rigid fiber in the Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 116 smooth slots, slightly more fibers go into the slot and fewer fibers are stapled. This result is also true for 1mm flexible fiber in C l slot as shown in Figure 6.39. It is interesting that 1mm flexible fibers are stapled more easily in slot C2 with low slot velocities than for high slot velocities (see Figure 6.40); but 1mm flexible fibers are stapled more easily in the L 2 slot with high slot velocities than for low slot velocities (see Figure 6.41). cM2 Vs=2.4mfs 1.5 2 2.5 3 Initial angle 8 (radians) 1 1.5 2 2.5 3 initial angle c> (radians) dri5 Vs=5.4 rrVs cfd5 Vs=5.4 m/s 1 1.5 2 Z 5 3 Initial angle e (radians) 1.5 2 25 3 Initial angle $ (radians) Figure 6.38 Smooth Slot: The effects of initial height and orientation6, <j> on passage o f 1mm flexible fiber for slot velocities 2.4 m/s and 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "o": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 117 Cl: cftC Vs=24 rrfs II t S 2 25 3 Initial angja 6 (radans) 0.5 1 __15 2 25 3 Initial anc^ e e (radans) Cl: cfcfi Vs=*4 rrfs Cl: cfcfi Vs=54 rrfs Initial arge e (radans) Initial angje co (radans) Figure 6.39 Cl Slot: The effects of initial height and orientation0, <p on passage of 1mm flexible fiber for slot velocities 2.4 m/s and 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "o": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 118 C2: dc2 Vs=24 rrYs C2: cfc2 Vs=2.4 tris Initial angle 6 (radars) Initial angle 6 (radars) C2: c f c S Vs=5.4 r r f s (2: d c E Vs=5.4 m s Initial angle 6 (radars) Initial angle o (radars) Figure 6.40 C2 Slot: The effects of initial height and orientation 0, <p on passage of 1mm flexible fiber for slot velocities 2.4 m/s and 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "p": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 119 L 2 c f c E V s = 2 4 iris 0.5 1 1.5 2 25 3 Initial angje e (radans) 12: c f cE V s = 2 4 r r f s 05 1 1.5 2 25 3 Initial arcje <|> (radans) L2 : c f c S V s = 5 4 r r f s L 2 c f c E V s = & 4 r r f s Initial antje e (radans) Inibal ancje $ (radans) Figure 6.41 L 2 Slot: The effects o f initial height and orientation 0, <p on passage o f 1mm flexible fiber for slot velocities 2.4 m/s and 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "a": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 120 6.3.4 Tra jec tor ies o f 3 m m flexible fiber (L/W = 6) Figure 6.42 shows similar fate plots of 3mm flexible fibers in the smooth slot for slot velocities 2.4 m/s and 5.4 m/s. More fibers are stapled in the main X Z plane with 9or/2 . Figure 6.43 and 6.44 are fate plots o f 3mm flexible fibers in C l and C 2 slots with slot velocity 5.4 m/s. For the C l slot, fewer fibers in the X Z plane become stapled and more fibers go into the slots, compared to the 3mm rigid fibers. For the C 2 slot, there are fewer fibers stapled in the X Z plane below the exit layer, but more fibers become stapled near the exit layer. The fate plots o f 3mm flexible fibers in L 2 slot is very similar to that o f the 3mm rigid fibers as shown in Figure 6.37. Smooth slot: cfd2 Vs=2.4 m/s Smooth slot: cfd2 Vs=2.4 m/s • " i •' fc': Initial angle Q (radians) Smooth slot: cfd5 Vs=5.4 m/s 1 i 1.3 1.4 1.5 1 15 2 25 3 Initial angle (|) (radians) Initial angle 6 (radians) 1.5 2 2.5 3 Initial angle (J) (radians) Figure 6.42 Smooth Slot: The effects o f initial height and orientation 9, <j> on passage o f 3mm flexible fiber for slot velocities 2.4 m/s and 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "o": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 121 ClsktcfrJ2 Vs=24 mfs Cl stat cfcE Vs=2.4 rrfs 1.3 1.4 1.5 1.6 1.7 1.8 1.9 Irtfaanrie 0 (radars) Cl slot cfd5 Vsr=5.4 mfs 15 2 25 3 Initial ande (j) (radars) a stat cfd5 Vs=5.4 mfs 1.1 12 13 1.4 1.5 1.6 1.7 13 1.9 2 21 inHaarge 6 (radans) 15 2 25 3 marge (j) (radans) Figure 6.43 C l Slot: The effects o f initial height and orientation9, (j> on passage o f 3mm flexible fiber for slot velocities 2.4 m/s and 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol "a": Stapled). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 122 (2 stat cfcE Vs=24 rrfs C2 stat cfcE Vs=24 rrfs intra age 9 (radans) 1 15 2 25 3 Htfaage (j) (radars) (2 stat cfcfi Vs=5.4 mfs C2 stat cfcfi Vs=3.4 mfs 1.1 12 1.3 1.4 1.5 1.6 1.7 18 13 2 21 iritiaarge 0 (radars) 2 2 5 3 iritiaage (I) (radars) Figure 6.44 C 2 Slot: The effects o f initial height and orientation 0, $ on passage o f 3mm flexible fiber for slot velocities 2.4 m/s and 5.4 m/s (symbols same as Figure 6.43). Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 123 L2sJotcfc2 Vs=Z4 mfs L2sictcfrJ2 Vs=24 mfs 1.5 1.6 1.7 18 inittaiande 9 (radians) L2 slot cfdS Vs=5.4 mfs 1.1 12 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 Intiaarrje 9 (radars) 15 2 25 3 iritiaJande ()> (radars) 15 2 25 3 Irirjal ande (j) (radars) Figure 6.45 L2 Slot: The effects of initial height and orientation 6, (j> on passage of 3mm flexible fiber for slot velocities 2.4 m/s and 5.4 m/s (Green symbol "+": Passes through the slot; Blue symbol "o": Goes away from the channel; Red symbol " Q " : Stapled). 6.4 Stapling From the above section, some fibers become immobilized by force balance and the forces may include hydrodynamic forces and the wall forces. In the calculation, it was assumed that the Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 124 fiber is stapled i f it does not leave through either the slot or the channel exit within a large number (e.g. 200,000) o f time steps. Gooding and Kerekes (1989), Olson (1996) and Kumar (1991) also observed the stapling phenomena in a single smooth slot. Kumar (1991) described stapling o f three types: vertical staple, horizontal staple and the mix o f both (see Figure 6.46) and observed that the slot to channel velocity ratio Vs I Vu, is the key variable in determining the type o f stapling. Vertical staples are formed when Vu is great than Vs and horizontal staples o f the 1mm rigid fibers only appear when Vs >VU. The 3mm rigid fibers stapled less often than 1mm rigid fibers, and short flexible fibers were observed not to staple under high slot velocities. A l l Kumar 's observations were based on the geometry o f a single smooth slot. Stapling in other slot geometries, based on the present numerical methods, w i l l be considered in this section. We w i l l discuss for each slot geometry, where and at what height and orientation stapling appears. (a) vertical staple (b) horizontal staple Figure 6.46 Vertical staple and horizontal staple. 6 . 4 . 1 S t a p l i n g i n t h e s m o o t h s l o t Both vertical and horizontal stapling can be observed in the smooth slot (see Figure 6.47). The 3mm rigid fibers tend to be horizontally stapled and 3mm flexible fibers tend to be stapled vertically at the downstream corner o f the slot. The 1mm flexible fibers are seldom stapled. Most o f the 1mm stapled rigid fibers are vertically stapled. But we also observed a few 1mm rigid fibers Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 125 vertically stapled at the open area o f the slot, which was not observed by Kumar (1991) but was also seen by Olson (1996). The turbulence intensity in the downstream corner changes considerably especially for the high slot velocity (see Figure 6.24). This may be the reason why stapling appears mostly at the downstream corner o f the slot. From the fate plots Figure 6.27 and Figure 6.28, most stapled 1mm rigid fibers appear at 8>KI2 for low slot velocities and at both 6>KI2 and 0 < n 12 for high slot velocities. There are a few fibers stapled at the orientation close to n i l . F rom the fate plot o f the 3mm rigid fibers in Figure 6.34, most 3mm rigid fibers stapled close to the wal l at 0 < n 12 . The 3mm flexible fibers tend to be stapled in the main plane with 0 < 1/2. smooth slot: cfd2 1mm rigid fiber • . . . • • - 3 - 2 - 1 0 1 2 3 X ( m m ) smooth slot: cfd5 3mm rigid fiber . -, t—, , , i , , • , 11 , , i H i i — — i . i . . . . smooth slot: cfd5 [ 1 1 1 1mm rigid fiber L • . . . i , • J • 1 1 . . J . 3 - 2 - 1 0 1 2 3 X (mm) : smooth slot: cfd5 3mm flexible fiber i i . i Figure 6.47 Fiber stapling in smooth slot. 6 . 4 . 2 S t a p l i n g i n t h e C l s l o t Figure 6.48 shows some examples o f fibers stapled in the C l slot. The 1mm rigid fibers staple at the downstream corner o f the slot and they can staple at any orientation and at a height very Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 126 close to the exit layer (see Figure 6.29). The stapled fibers separate the fibers that go into the slot from those that pass down the channel. Fewer 1mm flexible fibers become stapled (see Figure 6.39). A s shown in Figure 6.48, 3mm fibers staple at the downstream corner o f the slot or lie between the slot and the downstream corner o f the contour. From Figure 6.35, the 3mm rigid fibers tend to staple at 6<n12and cj> close to nil, which means the fibers are in or close to the main plane ( X Z ) . The 3mm flexible fibers still have the same stapling tendency as the rigid fibers although this similarity is less obvious (see Figure 6.43). C1 slot: cfd5 1mm rigid fiber _ i i _ X (mm) Figure 6.48 Fiber stapling in C l slot. 6 .4 .3 S t a p l i n g i n t h e C 2 s l o t Most o f the fibers stapled in the C2 slot at the downstream corner o f the contour (Figure 6.49), which is very different from the stapling in the C l slot where most fibers stapled at the Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 127 downstream corner of the slot. The two slots C l and C2 have different flow fields (see Figure 6.14). The stagnation streamline appears at the right corner of the slot in C l and the right corner o f the contour in C 2 . The fibers are l ikely to staple where there is a stagnation point. The 1mm rigid libers staple the same way as for the C l slot except that the height where fibers staple is at the exit layer which is a little lower than in the C l slot (Figure 6.31). Sti l l 1mm flexible fibers seldom staple in the C2 slot for the high slot velocity, but staple close to the wall or at the exit layer height for low slot velocity (Figure 6.40). For low slot velocity, 3mm rigid fibers staple at 6>n!2, which is different than in the C l slot. For high slot velocity, the stapled orientation o f 6 is the same as that in C l , but <j> is not close to nil above half o f the exit layer height (Figure 6.36). The orientation o f stapled 3mm flexible fibers has no preference in 6 and 4> (Figure 6.44). 0 1 2 3 _ 1 5.3 -2 -1 0 X (mm) X (mm) Figure 6.49 Fiber stapling in C2 slot. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 128 6.4.4 S t a p l i n g i n L2 s l o t Figure 6.50 shows some examples of the fibers stapled in the L2 slot. The 1mm rigid fibers seldom staple in the L2 slot as shown in Figure 6.32. The 1mm flexible fibers seldom staple in the L2 slot with low slot velocity but staple in this slot with high slot velocity (see Figure 6.41). The 3mm rigid and flexible fibers staple in about the same way (see Figure 6.37, 6.45), which is similar to the case of 3mm rigid fibers in the C l slot. L2 slot:jcfd5 1mm flexible fiber i h 1 / / / - // / L2 slot: cfd5 3mm flexible fiber i \ Figure 6.50 Fiber stapling in L2 slot. 6.5 Fiber passage ratio in different slot shapes This section uses the passage ratio formula (2.4) to calculate the passage ratio in several slot shapes. In order to use (2.4), two important variables need to be found: C(z) is the concentration profile upstream of the slot at height z, which is described in chapter 5. For simplicity, Olson's Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 129 fitting curve (2.5) is used for all the passage ratio calculations, although the calculated results are closer to Olson's original experimental data. For p(z), Olson (1996) assumed that the fibers initially within the exit layer come into the slot and the fibers initially beyond the exit layer pass away from the slot. So Olson approximated p(z) by (2.6). From current calculations, this approximation is about right for 1mm fibers in all geometries investigated. For 3mm fibers, this approximation is not accurate, e.g. for 3mm flexible fibers in the Cl slot (Figure 6.43), almost all the fibers initially below the exit layer enter the slot and there are also quite a number of fibers beyond the slot which pass through the slot. Lawryshry's formula is even less accurate according to the current calculations. For present purposes, we will not use any previous assumptions but will calculate p(z) directly as the ratio of the number of fibers going into the slot to the number of fibers at each height z . For this calculation, it should be mentioned that the stapled fibers in the smooth, Cl , LI and L2 slots are assumed to go into the slot. But the fibers stapled in the C2 slot are dealt with differently. If the center of the stapled fiber is close to the downstream corner of the contour, we assume that it passes away from the slot. 6.5.1 Passage ratio for 1mm rigid fiber Figure 6.51 shows the comparison of the currently calculated passage ratio of 1mm rigid fibers with other available results in smooth slot. The current calculation is a little under Olson's prediction. It shows that not all the fibers initially within the exit layer go into the slot. The fact that Lawryshyn's simulation is much lower than all the other results shows that his calculation with only one initial orientation is not accurate. The current calculation is very close to Kumar's experimental results for low slot velocities but is much lower than Kumar's results for high slot velocities. This Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 130 may be because that the turbulence effect is bigger for higher slot velocities and also the fiber-fiber interaction is important. These factors are not considered in the current simulation. 1.2 1 h .2 0 8 h ro or §,0 .6 (0 0. 0.4 0.2 - B - Smooth slot: current simulation Smooth slot: Kumar's experiment Smooth Slot: Olson's prediction Smooth Slot: Lawryshyn's simulation Fiber: 1mm, rigid <*• 0.2 0.4 0.6 0.8 1 Slot Velocity/Channel Velocity 1.2 Figure 6.51 Passage ratio of 1mm rigid fibers in smooth slot. Figure 6.52 shows the current results of passage ratio in Cl and C2 slots compared with Kumar's experiments. The current results are lower than the experimental results for low slot velocities, especially for the C2 slot. Stapling in C2 appears at two places: the downstream corner of the contour and the downstream corner of the slot. In defining the passage ratio, we treat the fibers stapling at the downstream corner of the slot as "going into" the slot and some fibers stapling at the downstream corner of the contour as "passing away". This may be inaccurate. In real cases, the stapled fiber may never go into the slot or pass away from the slot. Figure 6.53 is a comparison of the simulated passage ratio for 1mm rigid fibers in five different slot shapes: smooth, C l , C2, LI and L2. It is seen that slot shape has an important influence on the passage ratio. Not all of the contour slots have higher passage ratios than the smooth slot, e.g. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 131 C2 has a lower passage ratio than the smooth slot. For the step-step contour C2 and slope-slope contour L2, which have about the same height, L2 has a higher passage ratio than C2. Also, the passage ratio for LI is also higher than that for C l . This shows that the slope-slope contour has higher passage ratio than the same height step-step contour. For the same kind of contour C l and C2, the shallow step-step contour C l has a higher passage ratio than C2. And for the same kind of contour LI and L2, the shallow slope-slope contour LI has a higher passage ratio than L2. 1 4 1 .2 o | 0.8 • oi ro ro IL 0 .4 0 .2 C 1 s l o t : c u r r e n t s i m u l a t l o n C 2 s l o t : c u r r e n t s i m u l a t i o n C 1 s l o t : K u m a r ' s e x p e r i m e n t C 2 s l o t : K u m a r ' s e x p e r i m e n t F i b e r : 1 m m , r i g i d 0.2 0.4 0.6 0.8 1 S l o t V e l o c i t y / C h a n n e l V e l o c i t y 1 .2 Figure 6.52 Passage ratio of 1mm rigid fibers in C l , C2 slot. 1 .2 0.2 r S m o o t h s l o t : c u r r e n t s i m u l a t i o n C 1 s l o t : c u r r e n t s i m u l a t i o n C 2 s l o t : c u r r e n t s Im u l a t i o n L 2 s l o t : c u r r e n t s i m u l a t i o n L 1 s l o t : c u r r e n t s i m u l a t i o n 0.2 0.4 0.6 0.8 1 S l o t V e l o c i t y / C h a n n e l V e l o c i t y 1 .2 Figure 6.53 Passage ratio of 1mm rigid fibers in different slot shapes. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 132 Passage ratio of 1mm rigid fibers is also obtained in C l smooth slot, C l slope slot and Cltwoslope slot which have the same contour height as C l as shown in Figure 6.16. Table 6.3 shows the passage ratio for one slot velocity 5.4 m/s in those slots together with the results in the smooth slot, step-step slots (Cl and C2), and slope-slope slots (LI and L2). From this Table, we can make the following conclusions (Contour shapes are shown in Figure 6.2 and 6.16): Contour slots with contour height 0.5mm gives the biggest passage compared to smooth slot (0 height), C2 (1.2mm) and L2 (1.3 mm). For the slots with same contour height (e.g. 0.5mm) and same width (e.g. 3.0mm), downstream sloped side contours (Clslope, Cltwoslope) give nearly the same passage as C l . But, LI has more sloped downstream side contour and gives better passage than Clslope, Cltwoslope. Clsmooth has smaller width than C l and so provides lower passage. So when design a contour slot for better fiber passage, the focus should be on contour height and the angle of downstream side of the contour. Table 6.3 Passage ratio of 1mm rigid fiber in slots with different shapes. Contour height (mm) Slot shapes Passage ratio (Vs =5.4 m/s) 0 Smooth 0.6919 C l smooth 0.8145 C l 0.9099 0.5 Clslope 0.9133 Cl_two_slope 0.9111 LI 0.9821 1.2 C2 0.6840 1.3 L2 0.7742 Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 133 6 . 5 . 2 P a s s a g e r a t i o f o r 3 m m r i g i d f i b e r Similar to Figure 6.51, Figure 6.54 is the passage ratio of 3mm rigid fibers in the smooth slot obtained from available sources. Again, the current calculation is very close to Olson's prediction. Both current calculations and Olson's prediction are higher than Lawryshyn's simulation. But Kumar's experimental result keeps an "S" shape while all the others have a linear relationship. Because Olson's concentration curve was already verified in chapter 5, so the differences may result from differences in the function p(z). Olson's assumption is that the fibers initially within the exit layer will enter into the slot. The current calculation shows that not all the fibers within the exit layer enter the slot. If the "S-shape" appears, there would have to be many fibers beyond the exit layer which enter the slot, contrary to the results found in the current calculation, Olson's prediction and Lawryshrn's simulation. Only mean velocities were used in the passage ratio calculation and the fibers were calculated one by one, without affecting each other. These are simplifications of the actual experimental conditions. In the experiment, the flow is turbulent flow and the turbulence may draw the fibers beyond the exit layer into the slot. Stapling of many fibers was observed in Kumar's experiment for dilute suspension but this is impossible to observe in the current simulation because only a single fiber was calculated each time and no fiber-fiber interaction was considered. Figure 6.55 shows the passage ratio for 3mm rigid fibers in the Cl and C2 slots compared with Kumar's experimental results. The result for the C2 slot is close to the experimental results. But for Cl, the passage ratio for low slot velocities is higher than the experimental result and for high slot velocities it is lower than the experimental results. This may be explained in the same way as for the 3mm rigid fibers in the smooth slot. Comparing Figure 6.53 for 1 mm fibers in various slot shapes with Figure 6.56 for 3 mm fibers in the same slot shapes, it is apparent that the same trends are found for the passage ratio for Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 134 both fiber lengths, but that all the passage ratio amplitudes are lower for the 3 mm fibers. This is true because short rigid fibers go into the slots more easily than the long rigid fibers. 1 .2 0.2 0.4 0.6 0.8 1 S lot Velocity/C hannel Velocity 1 .2 Figure 6.54 Passage ratio of 3mm rigid fibers in the smooth slot. Figure 6.55 Passage ratio of 3mm rigid fibers in the C l , C2 slots. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 135 0.8 o 0.6 • u> 0) O) ro u> to ro 0.4 • Q-0.2 - H --e-Smooth slot: current simulation C 1 slot: current simulation C 2 slot: current simulation L 2 slot: current simulation Fiber: 3mm, rigid 0.2 0.4 0.6 0.8 1 Slot Velocity/Channel Velocity 1.2 Figure 6.56 Passage ratio o f 3mm rigid fibers in different slot shapes. 6 . 5 . 3 P a s s a g e r a t i o s f o r f l e x i b l e f i b e r s Figure 6.57 shows the passage ratio for 1mm flexible fibers in the smooth slot, compared with other available results. The 1mm flexible fiber has a higher passage ratio than the 1mm rigid fiber. The curve of the current simulation for 1mm flexible fibers (green with circles) lies between Kumar 's experimental result and Lawryshyn (1997)'s simulation result and is closer to the experiment. This again shows Lawryshyn (1997)'s assumption about p(z) is not suitable. Figure 6.58 compares the passage ratio o f 1mm flexible fibers in the C l and C2 slots with Kumar 's experimental results. The simulated results are closer to the experimental results for high slot velocities and are generally lower than the experimental results for low slot velocities. The reason may the same as that for 1mm rigid fiber. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 136 Figure 6.59 includes the passage ratio o f 3mm flexible fibers. They are much higher than that o f 3mm rigid fibers. Differences with Kumar 's experimental results occur for high slot velocities and the result is lower than the experiment. This may result from the fact that the fiber-fiber interaction effects may be bigger for long flexible fibers for high slot velocities. Figure 6.60 shows the passage ratio for 3mm flexible fibers in the C l a n d C2 slots. They are also greater than that o f 3mm rigid fibers, although the increased amplitude is smaller than that in Figure 6.59 o f smooth slot. There is no comparison with other results because no experiments or simulation data are available for the C l or C2 slots. Figure 6.57 Passage ratio of 1mm flexible fibers in the smooth slot. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 1.4 1 .2 rS 0.8 a cn co 0.6 co C L 0.4 0.2 -e-C1 slot: current simulation C2 slot: current simulation _ _ _ C1 slot: Kumar's experiment C2 slot: Kumar's experiment Fiber: 1mm, flexible 0.2 0.4 0.6 0.8 1 Slot Velocity/Channel Velocity 1.2 1.4 Figure 6.58 Passage ratio o f 1mm flexible fibers in the C l , C2 slots. 1.2 1 F .2 0 8 h <0 or §>0 .6 | -n CO 0- 0.4 |-0.2 h - B S m o o t h s l o t : 3 m m r i g i d f i b e r , c u r r e n t s i m u l a t i o n ~Q S m o o t h s l o t : 3 m m f l e x i b l e f i b e r , c u r r e n t s i m u l a t i o n _ — S m o o t h s l o t : K u m a r ' s e x p e r i m e n t f o r 3 m m r i g i d f i b e r S m o o t h S l o t : K u m a r ' s e x p e r i m e n t f o r 3 . 6 m m f l e x i b l e f i b e r S m o o t h S l o t : L a w r y s h y n ' s s i m u l a t i o n f o r 3 m m f l e x i b l e f ibeTj 0.2 0.4 0.6 0.8 Slot Velocity/Channel Velocity Figure 6.59 Passage ratio of 3mm rigid and flexible fibers in the smooth slot. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 138 1.2 .2 0 8 tt g.0.6 0.4 ro 0.2 - B --e-C1 slot: 3mm flexible fiber, current simulation C2 slot: 3mm flexible fiber, current sim ulation C1 slot: 3mm rigid fiber, current sim ulation C2 slot: 3mm rigid fiber, current simulation Fiber: 3mm, rigid and flexible 0.2 0.4 0.6 0.8 1 Slot Velocity/Channel Velocity 1.2 Figure 6.60 Passage ratio of 3mm rigid and flexible fibers in the C l , C2 slots. 6.6 Summary This chapter has examined in detail the flow field in slots with different shapes and the behavior o f fibers with different length and flexibility in those slots. The current results show that slot shape is the most important variable in determining fiber passage ratio. From the above passage ratio plots, slope-slope slot may be the best slot shape for all the slot shapes. For the same kind o f contour slot, the contour height is also very critical, and 0.5mm seems best for both the step-step contour and slope-slope contour. For the contour slots, the upstream side o f contour plays no significant effect on the fiber behavior and only the downstream side is important. Hydraulic resistance increases with the increase of turbulence intensity. Stapling may appear where there are large spatial gradients o f turbulence intensity. Chapter 6 Numerical simulation of the flow and fiber passage in slots with different shapes 139 The passage ratios for 1mm rigid and 3mm rigid fibers agree reasonably well with Olson's predictions. The passage ratios for 1mm rigid or flexible fibers and for 3mm rigid or flexible fibers agree better with Kumar's experimental results than Lawryshyn's simulation, which implies that the current simulation tool is more realistic. Differences between the current results and Kumar's experimental results may be due to the fact that the current simulation considers only the mean velocities and an isolated fiber simulation. Turbulence fluctuations and fiber-fiber interaction may be important in the passage ratio calculation, and further developments need to be made in these areas. Based on this chapter's results, the contour height and the angle of the downstream side of the contour are the two most important factors which should be paid more attention when designing a contour slot. Chapter 7 Other applications 7.1 Introduction The simulation tool developed in chapter 3 and 4 can be easily applied to other pulp and paper equipment. Section 7.2 introduces the application of the simulation tool to model the fiber orientation at the exit of the paper machine headbox. Section 7.3 displays another application of the simulation tool: fiber separation in a hydrocyclone. 7.2 Fiber orientation in the headbox Fiber orientation in the paper determines the quality of the paper and the study of fiber orientation in the paper machine headbox becomes more and more important. Zhang (2001) measured the fiber orientations along the centerline of an asymmetric headbox. He also predicted the fiber orientation in both symmetric and asymmetric headboxes using the simulated tool developed in this thesis, but only using the mean flow field obtained from the K-e model. His comparison between the experimental and predicted results showed that the simulated orientations tend to align more with the flow than the experimentally observed results because the turbulence was not taken into consideration for his calculation. Dong et al. (2002) added a turbulence effect in the prediction. Figure 7.1 is the mean flow field in a symmetrical headbox. Figure 7.2 is the preliminary results predicted by a "fixed" LES method, together with the results from the K-e model and Zhang's experimental measurements. A "fixed" LES method means that a one time step flow field is used for the calculation of the entire fiber orientation. It is shown that the result by the "fixed" LES is much closer to Zhang's experiment, which means that the turbulence is very important for the fiber orientation. This treatment of using "fixed" LES can dramatically decrease 140 Chapter 7 Other applications 141 the computation time, compared to the L E S method used in chapter 5 for the concentration calculation. Whether it can fully represent the turbulence effects needs further examination. 0 0.1 0.2 0.3 X (m) Figure 7.1 Grids and the mean flow field in a symmetry headbox. i n x z p l a n e x = 0 . 2 6 m a K - E m o d e l f i x e d L E S / \ Z h a n g ' s e x p e r i m e n t / \ 1 r ^ i ^ C T - . - p i f r - | a n g l e ( r a d ia n s ) i n x y p l a n e x = 0 .2 6 m K - E m o d e l f i x e d L E S S — Z h a n g ' s e x p e r i m e n t • ' 1 : L a n g le ( r a d ia n s ) Figure 7.2 Orientation in two planes near the headbox exit. Chapter 7 Other applications 1 4 2 7.3 Fiber separation in hydrocyclone The application of the developed methods to the fiber separation in hydrocyclone was included in the work of Wang (2002) and Wang et al. (2002). Hydrocyclones (another name is centrifugal cleaners) use centrifugal forces to separate solid particles from the suspending fluid. The centrifugal force is generated by swirling fluid motion. Generally, a hydrocyclone is composed of a cylindrical section and a conical section (Figure 7.3). It has two exits: an overflow opening and an underflow orifice. The feed suspension is injected tangentially near the top of the hydrocyclone. The flow moves toward the underflow exit along the inside wall of the hydrocyclone and rotates around the central axis. The flow velocity increases as the cone diameter decreases. When the flow approaches the underflow orifice, some of it moves out through the underflow orifice (rejects orifice), some flow rotates and moves upward in a smaller-diameter inner vortex and finally leaves the hydrocyclone from the overflow opening (accept orifice). Chapter 7 Other applications 143 Wang (2002) and Wang et al. (2002) studied the fiber motion and separation in two different hydrocyclones: Dabir's Hydrocyclone (Dabir 1983) and Bauer's Hydrocyclone (Ho 2001). A centrifugal force was added to the simulation tool developed in this thesis to model the fiber motion in the hydrocyclone. Three fiber fates are predicted: fibers moved along the wall downward and came out of the outflow orifice; fibers first moved downward and then upward and left through the overflow exit; or fibers moved up and down without going out from any exit. Figure 7.4 shows a comparison of the predicted feed, accepts and rejects coarseness to experimental observations (Ho 2001) for the same flow and fiber conditions. Coarseness is defined as the ratio of the total fiber mass to the total fiber length. A total of 50 fiberl and 50 fiber2 are used for this simulation. The characteristics of fiberl and fiber2 are listed in Table 7.1. It is shown that the simulation results agree well with experimental values for the accepts coarseness, but small differences appear for the rejects coarseness at lower flow rates. This may result from fiber flocculation which is more likely to happen at lower flow rates. The current fiber model does not consider fiber-fiber interaction. The essential tools developed in this thesis were also used by Wang (2002) to investigate the influence of fiber properties (density, diameter, length and flexibility) on separation efficiency. Chapter 7 Other applications 1 4 4 .30 0.75 ? E 2 0 5 a) c ro o O 0.25 40 - i — | i i 50 60 — r — i 1 — i 1 I i i i i 70 8 0. —I 1 1 1 1 — i 1 — a r e j e c t s - e x p e r i m e n t a l — b r e j e c t s - m o d e l a c c e p t s - e x p e r i m e n t a l B a c c e p t s - m o d e l R e j e c t s A c c e p ts 3 0 40 1 1 1 1 1 1 1 1 1 1 0.75 0.5 0.25 50 60 f l o w r a t e ( k g / m i n ) 70 80 Figure 7.4 Bauer's Hydrocyclone: comparison of mean coarseness in inlet and outlet streams (Wang 2002). Table 7.1 Properties of nylon fibers used in the experiments (Ho 2001) Fiber Fiber Fiber Density Length Diameter Coarseness (kg/m3) (mm) (um) (mg/m) Fiberl 1145 1 13.75 0.17 Fiber2 1145 1 27.50 0.68 Chapter 8 Summary and conclusions In this thesis, a complete numerical tool to simulate fiber motion in pulp and paper equipment has been developed. Three completed models are needed in this simulation: a fluid model to predict the flow field; a fiber model to predict the fiber motion; and a wall model to deal with the situation when a fiber touches the equipment. This tool was especially used to examine the flow and fiber behavior near a screen slot. The major findings of this thesis are summarized below: 8.1 Three models The fluid field in the screen was examined using two numerical methods: Reynolds Averaged Navier-Stokes (RANS) method for the slot shapes calculation (chapter 6) and Large Eddy Simulation (LES) for the concentration calculation upstream of the screen slot (chapter 5) . The standard k - e model was used to close the RANS equations. Smagorinsky's subgrid model was used for the LES method. A three-dimensional flexible fiber model was employed to track the fiber trajectory in the screen. The fiber was modeled as a chain of spheroids so that the computational time was lower than that for the previous model of the fiber composed of spheres. The present fiber model allows for the fiber flexibility, which is necessary for simulations in some pulp and paper equipment. A very general wall model was developed to deal with the fiber-wall interaction. The friction force at the wall is assumed proportional to (0.2 times) the normal force on the fiber. This wall model can be used when the fiber touches the wall of any geometry, without specifying the wall geometry, except through the flow calculation. 145 Chapter 8 Summary and conclusions 146 8.2 Fiber concentration upstream of the slot The distribution of fiber concentration was obtained using the flow computed in two ways: a Random Walk method was developed to use the k - E predictions, and a direct calculation of fiber motion was made, based on the detailed flow field computed by the LES method. The Random Walk method uses experimental RMS distributions for the three components of velocity close to the wall, as these components are not predicted separately by the k - e turbulence model. These two methods have been applied to predict the concentration of fibers in a fully developed three-dimensional channel flow which is assumed to be similar to the flow upstream of the screen. From the Random Walk method, with fibers of 1mm, 2mm, 3mm lengths, the concentration results scaled by the fiber length show a linear region near the wall and an almost constant region above a height of about one-half fiber length. Similar results were observed in the experiments of Olson (1996). A peak between the linear and constant concentration regions is observed in the simulation, a feature which also appears (with somewhat reduced magnitude) in Olson's (1996) data. Using the LES method with different length fibers, the concentration curves of the 2mm and 3mm fibers are similar to the experiments. The 1mm fiber results do not fit well with other computed or measured results, probably because the grid size for the LES calculation in the x direction (0.6mm) is too coarse to model the turbulence adequately at the scale of the smallest fibers. A finer grid could not be used for the current LES calculation because of computer memory limitations. The LES method is more general, as it does not rely on experimental data and gives a better representation of the turbulence characteristics, but it is computationally much more expensive. 8.3 The flow and fiber passage in slots with different shapes Only the mean flow field obtained from the RANS method and K-e model was used to examine the flow field in slots with different shapes and the behavior of fibers with different length and flexibility in those slots. Chapter 8 Summary and conclusions 147 Hydraulic resistance of the flow increases with the increase of turbulence intensity. Stapling may appear where there are large spatial gradients of turbulence intensity. Slot shape has the most important influence on fiber passage ratio. For the same kind of fiber, same flow condition and same slot width, the slope-slope slot provides the best passage for the fiber among all the slot shapes tested. For the same kind of contour slot, the contour height is also very critical, and 0.5mm seems best for both the step-step contour and slope-slope contour in the cases examined. The predicted passage ratios of 1mm rigid and 3mm rigid fibers agree reasonably well with Olson's predictions. The passage ratios for 1mm rigid or flexible fibers and for 3mm rigid or flexible fibers agree better with Kumar's experimental results than Lawryshyn's simulation, which implies that the current simulation tool is more realistic. Differences between the current results and Kumar's experimental results may be due to the fact that the current simulation considers only the mean velocities and an isolated fiber simulation. Turbulence fluctuations and fiber-fiber interaction may be important in the passage ratio calculation, and further developments need to be made in these areas. In conclusion, the present simulation tool allows us to predict fiber concentration and fiber passage through a screen slot. It could also be used to study the behavior of fibers in other complex turbulent flows. The optimization of the screen plate geometry using the present computation tool is now possible. Chapter 9 Recommendations for future work This thesis has developed a complete simulation tool to model fiber concentration and fiber behavior in the screen with different slot shapes. However, the tool still needs further development so that it can be used as a design tool. Some recommendations are as follows: The fiber-fiber interaction has not yet been considered in the current simulation. This is reasonable for only dilute pulp suspension. From this thesis, the results of passage ratio for longer fibers do not agree well with the experiment. This may result from that fact that longer fibers can easily interact with each other, even at low consistency. It would be better to include a fiber-fiber interaction in the present model to observe the behavior of long fibers at dilute suspension. Fiber-fiber interaction is certainly necessary for the application of higher feed consistency. Currently, only the RANS method was used to examine the flow and fiber performance in slots of different shapes. Therefore, only the effect of the mean flow field was taken into account. Turbulence effects must be examined by other methods. LES or Random Walk methods could be developed further to be applied to the simulation of the flow and fibers through the slots. In the Random Walk method, only one time scale was considered. Because turbulent flow is anisotropic near a wall, it has three unique time scales. The current technique is used only for simplicity, and could be improved by the use of different time scale in the three different directions. The flow field assumed in the present work is fully mixed upstream of the screen. This needs to be compared with experiments in an actual screen, and the effect of the rotor, not included here, should be added as an unsteady flow force. 148 Nomenclature a long axis of spheroid Aj(h) resistance tensors b short axis of spheroid C orbit constant Ci turbulence constant C2 turbulence constant Ca concentration of particles in the flow through the aperture Ci ( h ) resistance tensors CM turbulence constant Cu concentration upstream of the aperture Cf,c total contaminant concentration Cf.c.prob contaminant concentration by probability screening Cf.c.bar contaminant concentration by barrier screening C j a body-fixed connectivity vector C(z) upstream concentration of fibers at height z di unit directional vector along the fiber axis e spheroid eccentricity e. coordinate base E ( o o ) rate of ambient fluid strain tensor E removal efficiency 149 Nomenclature elastic modulus (fx,fy,fz) unit normal vector to the wall ft magnitude of the wall normal force F f fiber flexibility F normal force FT wall force Fi external force acting at spheroid i F.(g) body force acting at spheroid i F.o» hydrodynamic force acting at spheroid i F.(p) interparticle force acting at spheroid i g gravitational acceleration direction cosine matrix h time step Hex exit layer height H{*> resistance tensors I moment of inertia I, turbulence intensity K pressure drop coefficient L fiber length L, approximate size of local size of eddy Mi external torque acting at spheroid i m, mass of spheroid i M ; ( h ) hydrodynamic torque acting at spheroid N number of ellipsoids Nomenclature px, p2 average pressure at plane 1 and 2 P mean fluid pressure p filtered pressure p(z) probability of fiber passage at height z Pc passage ratio of contaminants Pe penetration number Pe penetration number based on the rotor tip velocity Vt Pf passage ratio of fiber Ap pressure drop q n , q n , qH, q n Euler parameters re aspect ratio R reject rate Rb Reynolds number based on bulk velocity RT Reynolds number based on friction velocity i-j position of spheroid i, measured from a fixed reference frame S, Stokes number SsliJj- stiffness Sy connectivity matrix SS mean strain-rate tensor (tx, ty, tz) unit tangential vector along the wall t time T residence time of the fiber T tangential friction force Nomenclature T, turbulence time scale Tc connectivity matrix T orbit period TN transforms fiber coordinate to the world reference TR transforms the world reference to fiber coordinate u( filtered velocity u instantaneous fluid velocity vector ur, vr, wr relative velocity Ui W ambient fluid translational velocity U mean velocity v(z) fluid velocity at height z V volume of the particle Vx, V2 average velocity at plane 1 and 2 V r ei relative velocity of the fiber with respect to the fluid Vs average flow velocity through the slot Vt rotor tip velocity Vu average upstream flow velocity W slot width x, y, z distance parameter xc,yc, zc fiber center position in the fixed reference coordinate X b internal constraint forces in joint b X c internal constraint forces in joint c XA, X c resistance function Nomenclature y+ non-dimensional distance yc half channel height Y b internal torques in joint b Y c internal torques in joint c YA, Y c , Y H resistance function Z,, Z 2 elevation at plane 1 and 2 Greek Letters a constant screen index ax, « 2 kinetic energy correction factor P constant screen index g. relative angular velocity S,- Kronecker delta e dissipation rate of turbulence eijk alternating tensor T] coefficient of wall friction cp angle between the Y-axis and the XY-projection of the fiber K kinetic energy of turbulence ju dynamic viscosity |IT eddy viscosity p density Nomenclature 154 9 angle between the fiber's major axis and the vorticity axis (Z axis o k turbulence constant oE turbulence constant rtj Reynolds-stress tensor T f fluid time scale rp particle response time 00| absolute angular velocity of spheroid i Q <°°) rate of ambient fluid vorticity Overstrikes vector Bibliography Adeniji-Fashola, A.A., Chen, C. P. and Schafer, C. F., Numerical predictions of two-phase gas-particle flows using Eulerian and Lagrangian schemes, NASA Technical Paper, 1988. Allison, B J . and Olson, J.A., Optimization of multiple screening stages for fibre length fractionation: two-stage case, J. Pulp Paper Sci., 26(3): 113-119, 2000. Boettcher P.C., "Results from a new design of contoured screen plate", Proc. 1986 Tappi Pulping Confi, 279-283, 1987. Brooke, J. W., Kontomaris, K., Hanratty, T. J. and McLaughlin, J.B., Turbulent deposition and trapping of aerosols at a wall, Phys. Fluids A 4, 825-834, 1992. Brooke, J. W., Hanratty, T.J. and McLaughlin, J.B., Free-flight mixing and deposition of aerosols, Phys. Fluids 6 (10), 3404-3415, 1994. Cash, J. R. and Karp, Alan H., A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides, A C M transactions on mathematical software, 16 (3): 201-222,1990. Chen, M. , Kontomaris, K. and McLaughlin, J. B., Dispersion, growth and deposition of coalescing aerosols in a direct numerical simulation of turbulent channel flow, Gas-Particle Flows, ASME-FED, 228, 27, 1995. Comte-Bellot, G., Ecoulement turbulent entre des sparois parallels, Publ. Sci. Tech. Minis. Air, NR, 1965. Crowe, C. T., Gore, R. A. and Troutt, T. R., Particle dispersion by coherent structures in free shear flows, Part. Sci. Tech. 3, 149, 1985. 155 Bibliography 156 Dabir, D., Mean Velocity Measurements in a 3 Inches Hydrocyclone Using Laser Doppler Anemometry, Ph.D. Thesis, Department of Chemical Engineering, Michigan State University, USA, 1983. Dong, S., Feng, X., Salcudean M . and Gartshore, I., Concentration of pulp fibers in 3D turbulent channel flow, to appear in Journal of Multiphase Flow, 2002. Dong, S., Feng, X., Salcudean, M. , Gartshore, I. and Shariati, M. , Turbulence and Fiber Orientation in the Converging Section of a Paper-machine Headbox, The 4th ASME/JSME/KSME Symposium on Computational Techniques for FluidVThermal/Chemical Systems with Industrial Applications, Vancouver, Canada, 2002. Dong, S., Salcudean, M . and Gartshore, I., Fiber motion in single and multiple screen slot, Proceedings of 2000 TAPPI Papermakers Conference and Trade Fair, Vancouver, Canada, 597, 2000. Feng, X., Large Eddy Simulation of the turbulent flow in a headbox, internal report, Department of Mechanical Engineering, University of British Columbia, 2001. Ferziger, J. and Peric, M . Computational methods for fluid dynamics, second edition, Springer, 1996. Fessler, J. R., Kulick, J. D. and Eaton, J. K., Preferential concentration of heavy particles in a turbulent channel flow, Phys. Fluids 6, 11, 3742-3749, 1994. Frejborg, F., "Profile screens improve TMP system efficiency and increase pulp quality", Pulp and Paper, 61(6), 99-101, 1987. Gooding, R. W., The passage of fibres through slots in pulp screening, M . A. Sc. thesis, University of British Columbia, 1986. Gooding, R. W., Flow Resistance of Screen Plate Apertures, Ph.D. thesis, University of British Columbia, 1996. Bibliography 157 Gooding, R. and Kerekes, R. J., The motion of fibres near a screen slot. J. Pulp paper Sci., 15 (2): 59-62, 1989. Gooding, R. W., Kerekes, R. J. and Salcudean, M., The Flow Resistance of Slotted Apertured in Pulp Screens, The Science of Papermaking, 1, 287-338, 2001. Gosman, A. D. and Ioannides, E., Aspects of computer simulation of liquid fueled combustors, AIAA paper, 81-0323, 1981. Haarlem, B. v., Boersma, B. J. and Nieuwstadt, F. T. M., Direct numerical simulation of particle deposition onto a free-slip and no-slip surface, Physics of Fluids, 10 (10), 2608-2620, 1998. Halonen L., Ljokkoi R. and Peltonen K, "Improved screening concepts", Proc. 1989 Patti Pulping Conf, 61-66, 1989. Ho, S-L, Synthetic Fiber Fractionation in Hydrocyclons, Ph.D. Thesis, University of British Columbia, Vancouver, B.C., 2001. Hopper, A.W., "The screening of chemical pulp", Pulp and Paper Manufacture, 3rd ed., vol. 5, ch. 12, TAPPI/CPPA, 1989. Jeffery, G. B., The motion of ellipsoidal particles immersed in a viscous fluid. Proc. Royal Soc, A102, 161-179, 1922. Kagermann, H. and Kohler, W. E., On the motion of nonspherical particles in a turbulent flow. Physica, 116A, 178-198, 1982. Kallio, G. A. and Reeks, M. W., A numerical simulation of particle deposition in turbulent boundary layers, Int. J. Multiphase Flow, 15 (3), 433-446, 1989. Kerekes, RJ. and Schell, C.J., Characterization of Fiber Flocculation Regimes by a Crowding Factor, J. Pulp Paper Sci., 18(1): 32-38, 1992. Bibliography 158 Kim, J. and Moin, P., Application of fractional-step method to incompressible Navier-stokes equations, J. Comp. Phys. 59, 308-323, 1985. Kim, S. and Karrila, S. J., Microhydrodynamics: Principles and selected applications, Butterworth-Heinemann, Boston, 1991. Kroger, C. and Drossinos, Y. A random-walk simulation of thermophoretic particle deposition in a turbulent boundary layer, Int. J. Multiphase Flow, 26, 1325-1350, 2000. Krushkal, E. M . and Gallily, I., On the orientation distribution function of non-spherical aerosol particles in a general shear flow-ii. The turbulent case. Journal of Aerosol Science, 19 (2): 197-211, 1988. Kubat, J. and B. Steenberg, Screening at low particle concentration, Svensk Papperstid., 58 (9), 319-324,1955. Kumar, A., Passage of fibers through screen apertures. Ph.D. thesis, University of British Columbia, 1991. Launder, B. E. and Spalding, D. B., The numerical computation of turbulent flows, Computer Methods in Applied Mechanics and Engineering, 3, 269-289, 1974. Lawryshyn, Y. A., Statics and dynamics of pulp fibers, Ph.D. thesis, University of Toronto, 1997. Liu, Z. C , Landreth, C. C , Adrian, R. J. and Hanratty T. J., High Resolution measurement of turbulent structure in a channel with particle image velocimetry, Experiments in Fluids, 10(6), 301-312,1991. Mason, S.G., "Fiber Motions and Flocculation", TAPPI Journal, 37(11), 494-501, 1954. McLaughlin, J. B., Aerosol particles deposition in numerically simulated channel flow, Phys. Fluid A 1, 1211, 1989. Bibliography 159 Moin, P. and Kim J., Numerical investigation of turbulent channel flow, J. Fluid Mech., 118, 341-377, 1982. Nelson, G.L., The Screening Quotient: A Better Index for Screening Performance, Tappi J., 64(5), 133-134, 1981. Nowak, P., A Multi-grid and Multi-block Method, Technical Report, The University of British Columbia, 1992. Olson, J. A., The effect of fiber length on passage through narrow apertures, Ph.D. thesis, University of British Columbia, 1996. Olson, J.A., A Lecture on Pressure Screening, Pulp and paper Center, Mechanical Engineering Department, University of British Columbia, 2000. Olson, J.A., Fibre length fractionation caused by pressure screening, slotted screen plates, J. Pulp Paper Sci., 27(8): 255-261, 2001. Olson, J.A., Allison, B J . and Roberts, N. , "Fibre length fractionation caused by pressure screening, smooth-hole screen plates", J. Pulp Paper Sci., 26(1): 12-16, 2000. Olson, J.A. and Kerekes, R.J., Fiber passage through a single screen aperture, Appita J. 51(2): 122-126, 1997. Olson, J.A., Roberts, N. , Allison, B.J. and Gooding, R.W., Fibre length fractionation caused by pulp screening, J. Pulp Paper Sci., 24(12): 393-397,1998. Olson, J.A. and Wherrett, G , A model of fibre fractionation by slotted screen apertures, J. Pulp Paper Sci., 24(12): 398-403, 1998. Pedinotti, S., Mariotti, G. and Banerjee, S., Direct numerical simulation of particle behaviour in the wall region of turbulent flows in horizontal channels, Int. J. Multiphase Flow, 18 (6): 927-941, 1992. Bibliography 160 Roberson J. A. and Crowe. C. T., Engineering Fluid Mechanics 3 r d ed., Houghton Mifflin, 1985. Rogallo, R. S. and Moin, P., Numerical simulation of turbulent flows. Ann. Rev. Fluid Mech. 16, 99-137, 1984. Ross, R. F. and Klingenberg, D. J., Dynamic simulation of flexible fibers composed of linked rigid bodies, J. Chem. Phys. 106 (7), 2949-2960, 1997. Rouson, D. W. I. and Eaton, J. K., Direct numerical simulation of turbulent channel flow with immersed particles, The 1994 ASME Fluids Engineering Division Summer Meeting, Lake Tahoe, Nevada, 1994. Shuen, J. S., Solomon. A. S. P., Zhang, Q. F. and Faeth, G. M. , A theoretical and experimental study of turbulent particle-laden jets, NASA CR 168293, 1983. Shuen, J. S., Solomon. A. S. P., Zhang, Q. F. and Faeth, G. M. , Structure of particle-laden jets: measurements and predictions, AIAA J. 23, 396-404, 1985. Skjetne, P., Ross, R. F. and Klingenberg, D. J., Simulation of single fiber dynamics, J. Chem. Phys. 107(6), 2108-2121, 1997. Smagorinsky, J., General circulation experiments with the primitive equations. Mon. Weather Rev. 91 (3), 99-164, 1963. Spalding, D. B., A single formula for the law of the wall, J. Appl. Mech., 28,455-457, 1961. Stockie, J. M . and Green, S. I., Simulating the motion of flexible pulp fibers using the immersed boundary method, J. Comp. Phys. 147, 147-165, 1998. Tangsaghasaksri W. Uber die sortierung von fasersuspensionen mittels geschlitzter siebe. Ph.D. thesis. Technische Hochschule Darmstadt, 1994. Bibliography 161 Tappi, "Screening Symbols, Terminology and Equations", Technical Information Sheets, TIS 0604-04, 1989. Thomas A. S. W. and Cornelius K. C , Investigation of a laminar boundary-layer suction slot, AIAA J. 20(6), 790-796, 1982. Wang, Q. and Squires, K. D. Large eddy simulation of particle-laden turbulent channel flow, Physics Fluids 8, 1207-1223, 1996. Wang, Z. 2002 Numerical simulation of fiber separation in Hydrocyclones, M . A. Sc. thesis, University of British Columbia. Wang, Z., Dong, S., Saucudean, M . and Gartshore, I. 2002 Fiber separation in hydrocyclones, 2002 TAPPI Fall Technical Conference and Trade Fair (Engineering and Pulping Conference), San Diego, CA. Wei, T. and Willmarth, W. W., Reynolds-number effects on the structure of a turbulent channel flow, J. Fluid Mech., 204, 57-95, 1989. Wherrett, G., Fiber motion in shear flow, M . A. Sc. thesis, University of British Columbia, 1996. Wherrett, G., Private communication, 1999. Wilcox, D. C. and Alber, I. E., A turbulence model for high speed flows, Proc. of the 1972 Heat Trans. & Fluid Mech. Inst., Stanford Univ. Press, 231-252, 1972. Wilcox, D. C. Reassessment of the scale determining equation for advanced turbulence models, AIAA Journal, 26(11), 1299-1310, 1988. Wittenburg, J., Dynamics of systems of rigid bodies, B. G. Teubner Stuttgart, 1977. Yamamoto, S. and Matsuoka, T., A method for dynamic simulation of rigid and flexible fibers in a flow field, J. Chem. Phys., 98 (1), 644-650,1993. Bibliography 162 Yamamoto, S. and Matsuoka, T., Viscosity of dilute suspensions of rodlike particles: A numerical simulation method, J. Chem. Phys., 100 (4), 3317-3324, 1994. Yu C. J. and Defoe R. J., Fundamental study of screening hydraulics. Part 1: Flow patterns at the feed-side surface of screen baskets: mechanism of fiber-mat formation and remixing, Tappi J. 77(8), 219-226, 1994a. Yu C. J. and Defoe R. J., Fundamental study of screening hydraulics. Part 3: Model for calculating effective open area, Tappi J. 77(9), 125-131, 1994b. Zhang, X., Fiber orientation in Headbox, M . A. Sc. thesis, University of British Columbia, 2001.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Modeling of fiber motion in pulp and paper equipment
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Modeling of fiber motion in pulp and paper equipment Dong, Suqin 2002
pdf
Page Metadata
Item Metadata
Title | Modeling of fiber motion in pulp and paper equipment |
Creator |
Dong, Suqin |
Date Issued | 2002 |
Description | This thesis describes a computational tool developed to simulate the fiber motion in pulp and paper equipment with particular emphasis on pulp screens. It is focused on modeling the fiber concentration upstream of a screen slot and examining the effects that different slot shapes have on the passage of fibers through the narrow apertures. This provides insight into the screen performance and has been used in this research to predict the characteristics of the screen. This simulation tool includes three coupled models: the flow model to predict the flow field in the equipment, the fiber model to trace the fiber trajectory in the equipment and the wall model to deal with the case when a fiber touches the equipment wall. The fluid field in the screen was examined using two numerical methods: the Reynolds Averaged Navier-Stokes (RANS) method for the calculation of flow through various slot shapes and Large Eddy Simulation (LES) for the fiber concentration calculation upstream of the screen slot. A three-dimensional flexible fiber model was employed to track the fiber trajectory in the screen. A very general wall model was developed to deal with the fiber-wall interaction. This wall model could be used for the case where a fiber touches a wall of any geometry, while specifying the wall geometry only through the fluid flow calculation routine. The simulated fiber concentration distribution increases linearly near the wall and becomes approximately constant farther from the wall, in reasonable agreement with experimental observations reported elsewhere. Slot shape has the most important influence on fiber passage ratio. For the same kind of fibers, same flow conditions and same overall slot width, a slope-slope slot provides the best passage for the fiber among all the slot shapes tested. For these contour slots, the contour height is also very critical, and a contour depth of 0.5mm seems best for both the step-step contour and slope-slope contour for fiber lengths of 1 to 3 mm. For the contour slots, the upstream side of contour plays no significant role in the fiber behavior; only the downstream side is important. |
Extent | 37944671 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-11-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0081004 |
URI | http://hdl.handle.net/2429/14756 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2002-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2003-792145.pdf [ 36.19MB ]
- Metadata
- JSON: 831-1.0081004.json
- JSON-LD: 831-1.0081004-ld.json
- RDF/XML (Pretty): 831-1.0081004-rdf.xml
- RDF/JSON: 831-1.0081004-rdf.json
- Turtle: 831-1.0081004-turtle.txt
- N-Triples: 831-1.0081004-rdf-ntriples.txt
- Original Record: 831-1.0081004-source.json
- Full Text
- 831-1.0081004-fulltext.txt
- Citation
- 831-1.0081004.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.831.1-0081004/manifest