UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Three-dimensional flow through forming fabrics and the motion of flexible fibres in flow Vakil, Ali 2010

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2010_spring_vakil_ali.pdf [ 7.63MB ]
Metadata
JSON: 24-1.0069055.json
JSON-LD: 24-1.0069055-ld.json
RDF/XML (Pretty): 24-1.0069055-rdf.xml
RDF/JSON: 24-1.0069055-rdf.json
Turtle: 24-1.0069055-turtle.txt
N-Triples: 24-1.0069055-rdf-ntriples.txt
Original Record: 24-1.0069055-source.json
Full Text
24-1.0069055-fulltext.txt
Citation
24-1.0069055.ris

Full Text

Three-dimensional Flow through Forming Fabrics and the Motion of Flexible Fibres in Flow  by Ali Vakil B.Sc., The Sharif University of Technology, Iran, 2002 M.Sc., The Sharif University of Technology, Iran, 2004  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY  in  The Faculty of Graduate Studies (Mechanical Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  January 2010 © Ali Vakil, 2010  Abstract This dissertation addresses the 3D flow through forming fabrics. In the first part of this thesis, the single phase 3D flow through certain specific fabrics was modeled. In practice, the Reynolds number of the fabric flow, based on paper-side filament diameter, is around 100. Consequently, it is not too surprising that the permeability of the fabrics was found to vary approximately as Re-0.4, which is intermediate between the expected low Re (Re-1) and high Re (Re0) limits. The resistance of a multilayer fabric was found to be nearly equal to the sum of the resistance of each layer considered in isolation. The effect of filament-scale and weave-scale flow non-uniformity on the fiber distribution in the finished paper was considered. For one specific fabric, there was 3 times more chance for short fibers to accumulate initially over openings than blockages of the fabric. Jet-to-wire speed ratio was found to have an insignificant effect on permeability results, but a marked effect on the Machine Direction shear stress in the vicinity of the paper-side filaments. In an attempt to model sheet formation, numerical simulations of the motion of a single fiber in the flow field of a cylinder was carried out as a canonical test case of the fiber/filament interaction system. Seven dimensionless groups govern the problem. A range of dimensionless groups were found for which the fiber hung up on the cylinder, whereas for other values of the dimensionless groups the fiber slid over the cylinder. In general, longer and more flexible fibers had a greater likelihood to be caught by the filament. Yawed finite aspect ratio cylinders at moderate Reynolds numbers are approximate representations of fibers in the flow field of forming fabrics. No analytical solution or any experimental data are reported in the literature to predict the drag and lift force on such particles in such flows. Computational Fluid Dynamics (CFD) simulations were conducted to find the drag and lift coefficients of inclined finite circular cylinders at Reynolds numbers in the range 1-40. The simulations showed that the Independence Principle was highly inaccurate for low inclination angles.  ii  Table of Contents Abstract ............................................................................................................................... ii Table of Contents............................................................................................................... iii List of Figures ..................................................................................................................... x Acknowledgments........................................................................................................... xvii Dedication ...................................................................................................................... xviii Co-Authorship Statement................................................................................................. xix Chapter 1............................................................................................................................. 1 Thesis Introduction ............................................................................................................. 1 1.1 Introduction............................................................................................................... 1 1.1 Literature Review...................................................................................................... 3 1.2.1 Literature Review on Forming Fabrics (Chapters 2 and 3) ............................... 3 1.2.2 Literature Review on Fiber Modeling (Chapter 4) ............................................ 7 1.2.3 Literature Review on Circular Cylinders (Chapter 5)...................................... 11 1.2 Problem Analysis .................................................................................................... 12 1.3.1 Forming Fabrics Modeling (Chapters 2 and 3)................................................ 15 1.3.2 Fiber Modeling (Chapter 4) ............................................................................. 17 1.3.3 Flow Around a Circular Cylinder (Chapter 5)................................................. 20 1.4 Thesis Objectives .................................................................................................... 21 1.5 References............................................................................................................... 22 Chapter 2........................................................................................................................... 29 Three-Dimensional Geometry and Flow Field Modeling of Forming Fabrics................. 29  iii  2.1 Introduction............................................................................................................. 29 2.2 Numerical Method .................................................................................................. 31 2.2.1 CAD Modeling of Fabrics ............................................................................... 31 2.2.2 Meshing and Boundary Conditions ................................................................. 38 2.2.3 Mesh Independence ......................................................................................... 40 2.2.4 Discretization Error.......................................................................................... 40 2.2.5 Code Validation ............................................................................................... 41 2.3 Results and Discussion ........................................................................................... 42 2.3.1 Velocity Contours ............................................................................................ 43 2.3.2 Pressure Drop Versus Flow Rate ..................................................................... 45 2.3.3 Velocity Averaging.......................................................................................... 46 2.4 Conclusions............................................................................................................. 53 2.5 Acknowledgments................................................................................................... 53 2.6 References............................................................................................................... 54 Chapter 3........................................................................................................................... 56 The Influence of Machine-Side Filaments and Jet-to-Wire Speed Ratio on the Flow Through a Forming Fabric ................................................................................................ 56 3.1 Introduction............................................................................................................. 56 3.2 Numerical Method .................................................................................................. 58 3.2.1 Meshing and Boundary Conditions ................................................................. 58 3.2.2 Mesh Independence ......................................................................................... 61 3.2.3 Discretization Error.......................................................................................... 62 3.3 Results and Discussion ........................................................................................... 63  iv  3.3.1 Single Layer Versus Multilayer Fabrics .......................................................... 63 3.3.2 Velocity Contours ............................................................................................ 64 3.3.3 Averaged Velocity Contours ........................................................................... 66 3.3.4 Pressure Drop Versus Flow Rate ..................................................................... 67 3.3.5 Oblique Versus Perpendicular Dewatering...................................................... 68 3.3.6 Pressure Drop Versus Flow Rate ..................................................................... 69 3.4 Conclusion .............................................................................................................. 71 3.5 Acknowledgement .................................................................................................. 71 3.6 References............................................................................................................... 72 Chapter 4........................................................................................................................... 74 Flexible Fiber Motion in the Flow Field of a Cylinder..................................................... 74 4.1 Introduction............................................................................................................. 74 4.2 Computational Methods.......................................................................................... 77 4.2.1 Single Phase Flow............................................................................................ 77 4.2.1.1 Flow Validation ........................................................................................ 78 4.2.1.2 Pure Water Flow Around a Cylinder ........................................................ 79 4.2.2 Fiber Modeling ................................................................................................ 80 4.2.2.1 Fiber in an Unbounded Field .................................................................... 83 4.2.2.2 Fiber in Contact with a Wall..................................................................... 83 4.2.2.3 Fiber Model Validation............................................................................. 86 4.2.2.3 Fiber Friction ............................................................................................ 90 4.3 Results and Discussion ........................................................................................... 91 4.3.1 Subset of Space State with Π 3 =  χ EI , Π5 = fc , Π 6 = as Variables... 95 3 D μU ∞ L v  4.3.2 Subset of Space State with Π 3 =  EI , Π 5 = f c , Π 7 = θ as Variables ... 96 μU ∞ L3  4.3.3 Subset of Space State with Π 1 =  ρU ∞ D EI χ , Π3 = , and Π 6 = as 3 D μ μU ∞ L  Variables ................................................................................................................... 98  ρU ∞ D EI , Π3 = , and Π 7 = θ as μ μU ∞ L3  4.3.4 Subset of Space State with Π 1 =  Variables ................................................................................................................... 99 4.3.4 Subset of Space State with, Π 2 =  χ EI L , Π3 = , and Π 6 = as Variables 3 D D μU ∞ L  ................................................................................................................................ 100 4.4 Summary and Conclusion................................................................................. 101 4.5 References......................................................................................................... 103 Chapter 5......................................................................................................................... 105 Drag and Lift Coefficients of Inclined Finite Circular Cylinders at Moderate Reynolds Numbers.......................................................................................................................... 105 5.1 Introduction........................................................................................................... 105 5.2 Numerical Method ................................................................................................ 108 5.2.1 Meshing and Boundary Conditions ............................................................... 108 5.2.2 Mesh Independence ....................................................................................... 111 5.2.3 Discretization Error........................................................................................ 112 5.2.4 Code Validation ............................................................................................. 113 5.3 Results and Discussions........................................................................................ 116 5.3.1 Flow Visualization ......................................................................................... 116  vi  5.3.2 Drag and Lift Coefficient............................................................................... 120 5.3.3 Normalized Drag Coefficient......................................................................... 123 5.3.4 Ratio of Lift to Drag Coefficient ................................................................... 124 5.3.5 Curve Fitting to Data ..................................................................................... 125 5.4. Concluding Remarks............................................................................................ 129 5.5. Acknowledgments................................................................................................ 129 5.6. References............................................................................................................ 130 Chapter 6......................................................................................................................... 132 Closing Remarks............................................................................................................. 132 6.1 Introduction........................................................................................................... 132 6.2 Summary of Chapter 2 .......................................................................................... 132 6.2.1 Further Comments ......................................................................................... 133 6.3 Summary of Chapter 3 .......................................................................................... 134 6.3.1 Further Comments ......................................................................................... 134 6.4 Summary of Chapter 4 .......................................................................................... 135 6.4.1 Further Comments ......................................................................................... 136 6.5 Summary of Chapter 5 .......................................................................................... 136 6.6 Recommendations and Future Work .................................................................... 137 6.7 References............................................................................................................. 139  vii  List of Tables Table 1.1 Effect of fabric design on drainage and paper properties ................................... 4 Table 1.2 Screen/grids resistance-flow rate relationships from experiments ..................... 6 Table 1.3 Fabric resistance to the flow through CFD simulations. .................................... 7 Table 1.4 Slender particle in a creeping flow of a Newtonian fluid ................................... 8 Table 1.5 Main methods to simulate flexible fibers in a flow field.................................... 9 Table 2.1 Air permeability comparison between simulations and experimental data. The applied pressure drop across the fabric is 125 Pa. Fabrics 2 and 3 were modeled with the new technique described here. cfm × 0.00508 = 1m s . ..................................................... 42 Table 3.1 General specification of Fabric A used in simulations. .................................... 60 Table 3.2 Comparison of the pressure drop through a full fabric and its layers............... 68 Table 4.1 Comparison between small deflection theory and the fiber model simulations of the free-end deflection of a cantilevered fiber in a channel flow. .................................... 90 Table 4.2 Salient variables that affect the interaction of a slender flexible particle with a rigid cylinder in a flow...................................................................................................... 91 Table 4.3 Dimensionless parameters for a slender particle in the flow of a cylinder...... 92 Table 5.1 Comparison of the lift and drag coefficient between different solvers, for L D = 10 , α = 40 o at Re = 40 ...................................................................................... 111  Table 5.2 Comparison of drag coefficient between asymptotic value of exponential curve fit to simulations and previous two-dimensional simulations [25] and experiments [24]. The experimental results were read from a logarithmic graph with limited resolution.. 114  viii  Table 5.3 Maximum and RMS difference between simulations and curve fit. .............. 126 Table 5.4 Comparison of C L and C D between simulations and curve fits for a specific simulation case that had not been used to derive the curve fit ( α = 45o , L D = 7 , Re = 7 ) ........................................................................................................................... 126  Table 5.5 Ratio of cylinder drag force predicted by simulations to that from the Independence Principle, for cylinders with various values of α, L/D, and Re................ 127 Table 5. 6 Lift coefficient curve fit................................................................................. 127 Table 5.7 Drag coefficient curve fit................................................................................ 128  ix  List of Figures  Figure 1.1: Schematic of fabric and impinging jet on a Fourdrinier machine.................. 13 Figure 1.2: (left) rows of circular cylinders as a simplified model of forming fabrics and relevant velocity and length scale, (right) different forces exerted on a fiber governing the fiber transient motion. The hydrodynamic interactions between particles are neglected in the force diagram. ............................................................................................................. 14 Figure 1.3: Computational geometry indicating the inlet, yarn and outlet volumes......... 16 Figure 1.4: Representation of a fiber as a whole in continuous form (left) and as a chain of N linked rigid bodies in discrete form (right) ............................................................. 17 Figure 1.5: Computational domain indicating the inlet, outlet, and slip walls ................. 20 Figure 2.1: Modeling of a forming fabric sample. (a) Prototype of a fabric (b) CAD model of the fabric. ........................................................................................................... 31 Figure 2.2: Typical x-ray cross section images (MD horizontal and ZD vertical) taken at varying depths along the CMD. ........................................................................................ 32 Figure 2.3: CMD-ZD image planes extracted from DATA VIEWER. ............................ 32 Figure 2.4: Identifying cross sections by the roundness metric........................................ 33 Figure 2.5: Radius of curvature is computed continuously along the border by tangent circles. ............................................................................................................................... 34 Figure 2.6: The concentration of identical tangent circles in the regions marked in blue indicates a detected yarn strand. ....................................................................................... 34  x  Figure 2.7: Projection of a 3D curve from raw data onto YZ plane. (a) Trace of filament in YZ plane (unprocessed) (b) Smooth 2D spline in YZ plane (c) Filament model viewed in YZ plane. ...................................................................................................................... 35 Figure 2.8: Projection of a 3D curve from raw data onto XZ plane. (a) Trace of filament in XZ plane (unprocessed) (b) Smooth 2D spline in XZ plane (c) Filament model viewed in XZ plane. ...................................................................................................................... 36 Figure 2.9: Individual filament that was modeled in Figure 2.7....................................... 37 Figure 2.10: Assembled complete model of a forming fabric. ......................................... 37 Figure 2.11: Computational geometry indicating the inlet, yarn and outlet volumes....... 38 Figure 2.12: Hybrid computational mesh scheme. ........................................................... 39 Figure 2.13: Upstream average velocity against mesh density for a pressure drop of 125Pa................................................................................................................................. 40 Figure 2.14: Top view of Integra structure, fabric 1. The units of the axes are metres. The dark elliptical area represents a typical fiber 1mm in length and 40 μm in diameter....... 42 Figure 2.15: Top view of Integra T structure, fabric 2. The units of the axes are metres. The dark elliptical area represents a typical fiber 1mm in length and 40 μm in diameter.43 Figure 2.16: Velocity contours of ZD, MD, and CMD velocities. Velocities in meters per second. (a) Integra (b) Integra T. .................................................................................... 45 Figure 2.17: C P − Re curve for fabric 1 and fabric 2....................................................... 46 Figure 2.18: Circular averaging over an opening and a knuckle, on a plane d/4 in the upstream of the Integra T fabric. (a) above an opening, (b) above a knuckle. ................. 48 Figure 2.19: Spatial variation of circular averaging over r = 0.2mm for Integra T fabric. ........................................................................................................................................... 49  xi  Figure 2.20: Spatial variation of ZD velocity when averaged over an ellipse with semimajor axes = (0.2,0.4,0.5,0.6)mm and semi-minor axis = 0.02mm . The ellipse is aligned with MD. Velocities in meters/second.............................................................................. 50 Figure 2.21: Spatial variation of ZD velocity when averaged over an ellipse with semimajor axes = (0.2,0.4,0.5,0.6)mm and semi-minor axis = 0.02mm . The ellipse is aligned with CMD. ........................................................................................................................ 51 Figure 2.22: Effect of the aspect ratio on the averaged maximum and minimum velocities. The ellipse is aligned with (a) MD (b) CMD. ................................................. 52 Figure 3.1: Schematic of fabric and impinging jet. .......................................................... 58 Figure 3.2: Hybrid computational mesh scheme. ............................................................. 60 Figure 3.3: Computational geometry indicating the inlet, yarn and outlet volumes......... 61 Figure 3.4: Upstream average velocity against mesh density for a pressure drop of 125Pa. Flow is normal to the fabric. ............................................................................................. 62 Figure 3.5: Top view of Fabric A ..................................................................................... 63 Figure 3.6: Top view of layers composing Fabric A. ....................................................... 64 Figure 3.7: Pointwise velocity contours of ZD, MD, and CMD velocities. Velocities in meters per second. (a) full Fabric A (b) Paper-side layer of Fabric A.............................. 65 Figure 3.8: Spatial variation of the difference between ZD velocities when averaged over an ellipse with semi-major axis = 0.2mm and semi-minor axis = 0.02mm. The difference is normalized by the mean velocity. The ellipse is aligned with MD............................... 66 Figure 3.9: Spatial variation of the difference between ZD velocities when averaged over an ellipse with semi-major axis = 0.2mm and semi-minor axis = 0.02mm. The difference is normalized by the mean velocity. The ellipse is aligned with CMD. ........................... 67  xii  Figure 3.10: C P − Re curve for full Fabric A and its comprising layers. ........................ 68 Figure 3.11: Spatial variation of non-dimensional shear stress when averaged over an ellipse with semi-major axis = 0.2mm and semi-minor axis = 0.02mm. The ellipse is aligned with MD. Owing to the unphysical imposed zero shear stress wall condition applied to MD=2.5mm and MD=1.4mm, only the region of the fabric away from these walls is shown. .................................................................................................................. 69 Figure 3.12: C P − Re curve for Fabric A with impingement angle as a parameter. ........ 70 Figure 4.1: Schematic of linting precursor in forming section. Cross section through pulp mat and forming fabric...................................................................................................... 74 Figure 4.2: (a) Computational domain indicating the boundary conditions. (b) Close-up view of structural mesh around the cylinder..................................................................... 78 Figure 4.3: Effect of the blockage ratio on the drag coefficient of a cylinder at Re = 0.5. ........................................................................................................................................... 79 Figure 4.4: Flow streamlines and field variables around a single cylinder at Re = 1....... 80 Figure 4.5: Representation of a fiber as a series of N linked rigid bodies. ..................... 81 Figure 4.6: A prolate spheroid representing each segment of a fiber. .............................. 81 Figure 4.7: Schematic of a fiber in contact with a solid surface....................................... 84 Figure 4.8: A statically indeterminate stationary system which is split into two subsystems as shown in (b) and (c). In each subsystem, the segments in contact with the wall are assumed to be stationary, but the rest may hang around them............................ 86 Figure 4.9: Comparison between Jeffery’s orbital periods and the fiber model simulations in a simple shear flow. ...................................................................................................... 87  xiii  Figure 4.10: A cantilever beam under a uniform transverse loading and its counterpart in the discrete form. .............................................................................................................. 88 Figure 4.11: Ratio of the free end deflection from the small deflection theory to that predicted by the fiber model simulations.......................................................................... 88 Figure 4.12: A suspended fiber is a channel flow with a parabolic velocity profile. ....... 89 Figure 4.13: A belt winding around a circular object and its counterpart in discrete form. ........................................................................................................................................... 90 Figure 4.14: Typical configurations of a fiber in the proximity of a cylinder. The angle θ and offset χ are measured far upstream of the cylinder. ................................................. 92 Figure 4.15: Motion of a nearly rigid fiber near a cylinder. The fiber is shown at unequal time increments. Π 1 = 1 , Π 2 = 40 , Π 3 = 1.25 , Π 4 = 8 , Π 5 = 0.2 , Π 7 = 5 o ................. 94  Figure 4.16: Fiber hang-up on a cylinder for different values of flexibility with f c = 0.2 , and Re = 1 . The fiber stiffness in figures (a), (b), and (c) is respectively EI = [80,25,5]× 10 −12 N .m 2 . .............................................................................................. 95  Figure 4.17: The hung-up regions of the space state as a function of Π 3 , Π 5 , and Π 6 . . 96 Figure 4.18: Fiber hung-up on a cylinder for different values of flexibility with f c = 0.2 , and Re = 1 . The fiber stiffness from left to right is EI = [140,40,2]× 10 −12 N .m 2 , and  θ = [13o ,16 o ,40 o ] .............................................................................................................. 96 Figure 4.19: The hung-up region for as a function of Π 3 , Π 5 , and Π 7 . ......................... 97 Figure 4.20: Snapshots of a fiber sliding on a cylinder (a) rigid fiber (b) flexible fiber. Flow Reynolds number is 1, and EI = [80,2]× 10 −12 N .m 2 from left to right. .................. 98 Figure 4.21: The hung-up region as a function of Π 1 , Π 3 , and Π 6 . .............................. 99  xiv  Figure 4.22: The hung-up region as a function of Π 1 , Π 3 , and Π 7 . .............................. 99 Figure 4.23: The hung-up region as a function of Π 2 , Π 3 , and Π 6 . ............................ 100 Figure 5.1: Computational domain indicating the inlet, outlet, and slip walls. .............. 108 Figure 5.2: Velocity profile along the dashed line at Re = 40 , showing the negligible influence of the wall........................................................................................................ 109 Figure 5.3: Velocity profile along the dashed line at Re = 1 , showing the negligible influence of the wall........................................................................................................ 109 Figure 5.4: Hybrid computational mesh used in this research........................................ 110 Figure 5.5: Unstructured mesh used close to the cylinder. ............................................. 110 Figure 5.6: Drag coefficient versus number of mesh volumes for L D = 2 , Re = 40 ,  α = 90 o . Note the highly expanded ordinate axis........................................................... 112 Figure 5.7: Drag coefficient versus number of mesh volumes for L D = 5 , Re = 1 ,  α = 10o . Note the highly expanded ordinate axis........................................................... 112 Figure 5.8: Drag coefficient versus Reynolds number ( 0 p Re ≤ 40 ) for a cylinder normal to the flow with L D as a parameter. ................................................................ 114 Figure 5.9: Comparison of drag coefficient between simulations and Cox, Re = 1 . ..... 115 Figure 5.10: Streamlines around a cylinder with L D = 20 at different inclination angles. The Reynolds number is 40. ........................................................................................... 117 Figure 5.11: Streamlines around a cylinder with L D = 2 at different inclinations angles. The Reynolds number is 40. ........................................................................................... 118 Figure 5.12: Streamlines around a cylinder with different aspect ratios and inclination angles. Re = 1.................................................................................................................. 119  xv  Figure 5.13: Nondimensional pressure distribution along a cylinder with L D = 20 for different inclination angles at Re = 40 . (a) Front pressure (b) Rear pressure................ 120 Figure 5.14: Forces on a finite cylinder with ( L D = 2, 5, 10, 20) versus the angle of inclination with Re as a parameter. a) Drag coefficient b) Lift coefficient................... 121 Figure 5.15: Ratio of pressure drag to viscous drag as a function of inclination angle at two different aspect ratios at (a) Re = 1 , (b) Re = 40 .................................................... 122 Figure 5.16: Nondimensionalized torque as a function of the inclination angle for different aspect ratios. (a) Re = 1 , (b) Re = 40 . ............................................................. 123 Figure 5.17: Non-dimensional drag coefficient ( L D = 2, 5, 10, 20) versus the angle of inclination with Re as a parameter................................................................................. 124 Figure 5.18: Lift-to-drag ratio on a finite cylinder with ( L D = 2, 5, 10, 20) versus the angle of inclination with Re as a parameter................................................................... 125  xvi  There is no duty more indispensable than that of returning a kindness. Cicero  Acknowledgments I want to pay special tribute to my supervisor Prof. Sheldon I. Green for his technical guidance and exceptional patience during my research. I am grateful to my committee members Prof. Richard J. Kerekes, Prof. Mark Martinez and Prof. Carl Ollivier-Gooch for their helpful input and valuable discussions. The administrative support of Pulp and Paper Centre staff, especially Ken Wong, Tim Paterson, and Lisa Hudson and Mechanical Department staff, especially Alan Steeves, Jay Zhao and Yuki Matsumura are highly appreciated. Thank you to all the members in my research group and my friends, innumerous to be counted, for being a source of inspiration and motivation for me. AstenJohnson Inc., NSERC, and Du Pont Canada Inc. are gratefully acknowledged for their financial supports. The helpful hands and critical comments of my beloved brother Dr. Mohammad Vakil and my friend Dr. Poorya Ali Ferdowsi have been a treasure in my research. I thank to them for the many conversations which have enhanced my understanding in all aspects of my research. I came to know that ever-increasing frustrations and aspirations are an indispensable constituent of any research. The never-ending encouragements and words of wisdom of my brother and my mother and the memory of my father kept me paddling my spinning will through the river of self-doubt.  xvii  Dedication To the relentless encouragement and support of my mother and brother, and to the loving memory of my father.  Swilling wine and being cheerful are my morality. Disengaging from blasphemy and religion is my formality. I asked the bride of fortune: What is her heart demand? She replied: Abandoning thy confident heart is my prosperity! Omar Khayyám  xviii  Co-Authorship Statement Authorship of Chapter 2 The authors of chapter 2 are: Ali Vakil, Arash Olyaei, Sheldon Green. My contributions to this paper are described below.  Identification and design of research program  Green was originally interested in studying the hydrodynamics of pulp suspension filtration through forming fabrics. He previously conducted simulations of 2D models of fabrics and was interested in simulating 3D flow through fabrics. The manufacturer of the fabrics and founder of the project, AstenJonhson Inc., was interested in Rapid Prototype Models (RPM) of their fabrics. Green’s research group had developed a method to build 3D CAD models based on micro-machining. With new development in the research group, we were able to produce highly accurate 3D models of fabrics. Green was also interested in characteristics of the flow field through fabrics. I designed the computational domain and conducted the numerical simulations of the flow.  Performing the research  Olyaei developed a novel method to build 3D CAD models of forming fabrics from 2D X-ray images. I produced computational domains for the simulations using Gambit and preformed the simulations using the Computational Fluid Dynamics (CFD) package Fluent. I ran the flow simulations for the different sets of parameters needed to gather the data for further analysis.  Data Analysis  With direction from Green, I studied the flow field results and analyzed the data to present the flow characteristics. Green also stated a simple theory to correlate the initial fiber mass distribution on fabrics to the flow field upstream of the paper-side layer. I conducted the studies based on the proposed theory and produced the whole figures.  xix  Manuscript preparation  Olyaei wrote the CAD modeling section of the paper and I wrote the numerical modeling and discussion sections. The whole manuscript was edited and revised by Green.  Authorship of Chapter 3 The authors of chapter 3 are: Ali Vakil, Arash Olyaei, Sheldon Green. My contributions to this paper are described below.  Identification and design of research program  An interesting question for the papermaking industry, since the development of the multilayer fabrics, was the effect of the machine-side layer on the flow field through the paper-side layer. The effect of the Jet-to-Wire speed ratio on flow characteristics has also been a matter of interest to papermakers. Green simplified the problems into a tractable one, and I conducted the simulations using methods similar to those presented in chapter 2.  Performing the research  I designed the computational domains, based on the 3D CAD models produced by Olyaei, and ran the simulations for different fabrics. I gathered the data to study the questions described above.  Data Analysis  With direction from Green, I analyzed the data to compare flow fields of multilayer fabrics with those of its comprising layers. I conducted an analysis of the flow regarding the effect of Jet-to-Wire speed ratio with Green’s advice.  Manuscript preparation  I wrote the entire manuscript with revisions and suggestions from Green. xx  Authorship of Chapter 4 The authors of chapter 4 are: Ali Vakil, Sheldon Green. My contributions to this paper are described below  Identification and design of research program  The fiber deposition and dewatering process occurring on a forming machine in the sheet formation process is a highly complex phenomenon. Fibers in papermaking vary in dimensions and stiffness. Therefore, a method should be adopted to model both flexible and rigid fibers. Green proposed using the Particle Level Simulation (PLS) method, originally developed by Klingenberg et al. I followed the PLS approach and modified the method to be used in the fiber deposition problem. We considered a highly diluted pulp suspension and only considered the interaction of one fiber with one filament of fabric.  Performing the research  I rearranged the equations to be suitable for both unbounded and bounded flow fields. Based on the modified equations, I wrote a FORTRAN code to track flexible fibers in given flow fields. I ran the flow simulations in Fluent and formatted the data to be fed into the FORTRAN code. The problem is governed by several dimensionless groups. I conducted the fibre simulations by varying different groups to classify a set of “hung-up space”. These spaces are to identify whether a fiber, upon impact on a cylinder, slides or sticks to the cylinder.  Data Analysis  With direction from Green, I analyzed the data based on the 7 dimensionless groups involved, and presented the fiber trapping regions in graphs.  Manuscript preparation  I wrote the entire manuscript with revisions and suggestions from Green.  xxi  Authorship of Chapter 5 The authors of chapter 5 are: Ali Vakil, Sheldon Green. My contributions to this paper are described below.  Identification and design of research program  Fibers, in papermaking conditions, move relative to the background water with Reynolds numbers as high as 40. Surprisingly, the drag and lift coefficients of inclined, finite cylinders at low to moderate Reynolds numbers did not exist in the literature. Green proposed the idea of conducting simulations to document such results. I conducted the numerical simulations of the flow.  Performing the research  I set up the computational domains and ran the simulations for numerous parameters involved.  Data Analysis  With instruction from Green, I analyzed the flow field and presented the drag and lift coefficient graphs as a function of the desirable dimensionless groups. I also found the best-fit curves for the data, with help from Green.  Manuscript preparation  I wrote the entire manuscript with revisions and suggestions from Green.  xxii  Chapter 1 Thesis Introduction 1.1 Introduction  The pulp slurry exiting a headbox is a mixture of fibers, fines, fillers, and additives in water. In the forming section, forming fabrics are employed to remove the water mechanically from this slurry to form a wet sheet of paper. This sheet is then pressed and dried to make paper. Forming fabrics typically consist of woven polyester filaments about 150 μm in diameter, with a pitch of 300 − 400μm . Typical pulp fibers have a diameter of 20 − 40μm and are 100 − 3000μm long, and pulp fines are fiber fragments less than 100 μm in length. Thus, paper is an intricate composite of slender particles formed on a complex network of strands. Historically, the art of papermaking dates back to A.D. 105 in China. The practice of papermaking as a continuous process, as we know it today, was invented by Nicholas Louis Robert (1798), and financed by the Fourdrinier brothers (1807) [1]. The process of transforming pulp into paper has undergone many changes and modifications during the past two centuries. Before the 1970’s, all forming sections (Fourdrinier machines) used one horizontal wire to drain water from the slurry. Drainage in such Fourdrinier formers was accomplished by inertial impingement, gravity, and drainage elements (table rolls, foils, suction boxes) [2]. These first generations of forming sections due to their limitations [3,4], including lack of shear affecting the top surface and excessive free surface instability, were superseded by twin-wire/gap formers. In these formers dewatering takes place through both the outer and inner fabrics. The invention of twinwire/gap formers fulfilled the need for increased speed, reduction of two-sidedness, and improved formation. The old single-layer forming fabrics were also supplanted by new multilayer fabrics to increase (i) structural integrity on the machine-side, in contact with dewatering 1  devices, by using a coarse mesh and thick filaments, (ii) sheet smoothness, by using finer mesh and filaments on the fabric side in contact with the pulp suspension. The main requirements of forming fabrics are high drainage capacity, superior retention and better paper formation. However, there is a compromise between the last two, as Fourdrinier machines benefit from good retention, whereas twin-wire/gap formers may have less retention. This thesis considers the influence of the non-uniform flow entering the paperside of fabrics on non-uniformity in papers. “Uniformity is an indication of good formation”[5]. Computational Fluid Dynamics (CFD) simulations of the flow through three-dimensional CAD models of real forming fabrics are first presented (chapters 2,3). Surprisingly good agreement between the predicted flow and bulk flow experimental measurements demonstrate the utility of the simulations. In each of these chapters probable correlations between the flow non-uniformity and initial fiber/fine mass over different areas of fabrics are presented. The effect of machine-side layer on the upstream flow of the paper-side layer and the jet-to-wire speed ratio are dealt with in chapter 3. In chapter 4, as an initial step to simulate the sheet formation process, an attempt is made to predict the motion of a single flexible fiber in the flow of a circular cylinder, as a simplified model of one forming fabrics filament. Seven dimensionless groups governing the fiber-cylinder interaction are defined. These groups are then used to classify the likelihood of fibre interception by the cylinder into hung-up and non-hung-up regions. Pulp fibres are approximately cylindrical in shape, with an aspect ratio in the range 2-150, and move at Reynolds numbers of 1-40 relative to the water. There is surprisingly little information for the flow around a cylinder in the moderate Reynolds number range. Chapter 5 fills this gap in the literature.  2  1.1 Literature Review  As this thesis is comprised of several stand-alone, yet interrelated papers, it is appropriate to present the literature review of each section separately as follows.  1.2.1 Literature Review on Forming Fabrics (Chapters 2 and 3)  In the forming section of a papermachine, pulp is spread over a moving fabric and water is drained mechanically through the fabric, leaving the fibers behind to form a wet sheet of paper. The mechanical dewatering which commences in the forming section plays an important role in the final paper quality as any non-uniformities arising in this section cannot be corrected later in the papermachine [6]. Therefore, understanding the flow field through the forming section provides insight into the finished paper. Helle summarizes the historical development of forming fabrics from the “bronze age (1914-1951)” to the “synthetic age (1960-1990)” [7]. A summary of fabric design and its effect on drainage and sheet formation is presented in Table 1.1. Good forming fabrics are characterized by a low running resistance and good fiber retention. Multi-layer fabrics meet the retention requirement better than a singlelayer one [8,9]. Multilayer fabrics also optimize operating life and sheet quality. The former is met by means of a coarse filament size on the machine-side layer, and the latter by using finer filaments on the paper-side layer [10]. The non-uniform distribution of the fibers in the sheet, the “formation”, is a quantitative measure of sheet quality. Fibers deposition on and penetration into fabric openings result in a pattern which remains in the finished paper, known as topographical wire mark [11]. It is also postulated that the flow non-uniformity produced by the  structure of a fabric can unevenly distribute the pulp fibers, producing hydrodynamic wire mark [12]. This weak systematic variation in the distribution and orientation of the  fibers can be distinguished from much stronger random variation by examining the frequency spectrum of the printed structure [12,13]. As stated by Helle “the human eye has a pronounced sensitivity to the systematic variations of paper surface associated with wire mark.”  3  It is known that the wire-mark non-uniformity can have a profound effect on the quality of the printed end product [18,19]. Uneven ink penetration resulting from the wire mark can be detected as peaks in the frequency spectrum obtained by FFT [12,13]. The fiber orientation distribution in a sheet arises from two sources: fiber orientation in the headbox and fiber redistribution that occurs on and after jet impingement on the fabrics. It is known that fiber orientation before and upon impingement has a large effect on the final paper properties [5,14]. Paper properties such as strength, smoothness, tensile and tear ratio are correlated with fabric characteristics such as air permeability, drainage capacity, and fiber support [5]. Table 1.1 Effect of fabric design on drainage and paper properties Author Description Conclusion Helle (1978) [14]  Beran (1979) [15] Johnson (1984) [8] Johnson (1986) [9] Adanur (1994) [5]  Helmer et al. (2006) [16]  Xu et al. (2009) [17]  Investigated the effect of long knuckle orientation on drainage and sheet formation. (single-layer fabric) Introduced a fiber support index (FSI) based on mean fiber length and the average number of supports per fiber. (single-layer fabric) Introduced a drainage index (DI) based on fabric air permeability, CD mesh counts and the FSI. Using DI predicted retention and drainage behaviour of multilayer fabrics. Investigated effects of fabric characteristics on sheet formation and properties. Introduced formation index (FI) as a function of “crowding number”, “Froude number” and “aspect ratio of flow depth to fiber length”. Studied the effect of internal structure of triple fabrics on drainage behaviour.  4  Faster drainage and better release when fibers oriented across the long knuckles. FSI can determine optimum fabric designs for a particular paper machine. Higher FSI corresponds with papers with finer grades. DI is used to assess drainage performance. Higher values of DI reflect better drainage. Inverse relationship of retention to MD frame length; proportional relationship of drainage to DI. An increase in DI and FSI resulted in an increase in sheet strength; An increase in mesh count reduced sheet thickness. Characterized dewatering as “early” and “late” drainage velocity; Higher late drainage rate refers to higher FI and higher topside anisotropy. An increase in the center plane resistance increases the total resistance only at lower basis weight.  Monofilament grids and screens may be considered as simplified models of single-layer forming fabrics. There has been a great number of publications relevant to high Re flows through screens [20,21,22]. To a lesser extent, the literature contains work on screens at low to moderate Re. Table 1.2 lists some of the important work done experimentally on screens and their semiempirical resistance-flow rate relationships. All papers attempted to relate the flow resistance to the pore geometry and fabric structural properties as well as fluid properties. It is conventional to define a pressure drop coefficient for fabrics, K , and present the data in the form K − Re . This coefficient, K , is usually defined as the ratio of static pressure drop across the screen against the dynamic pressure of the upstream ( K =  Δp 1 ρU ∞2 2  ). It should be noted that researchers have  used different definitions of Reynolds number for screens: (i)  Based on the filament diameter as Re d = ρU ∞ d μ .  (ii)  Based on the ratio of projected open area of the screen to the total crosssectional area, “porosity” β , as Re p = Re d β .  (iii)  Based on the surface area exposed to the flow per unit volume, S , Re s = Re d (Sd ) .  (iv)  Based on the wetted perimeter of the orifice where the flow is most constricted,  (  W,  and  the  filaments’  center-to-center  distance,  l,  )  Re w = 4l 2 Wd Re d . where ρ is the fluid density, μ is the dynamic viscosity, and d is the screen filament diameter. Rushton and Griffiths showed that the analysis based on the orifice model correlates best with experimental results among the other proposed methods [30]. These models for single-layer fabrics are not applicable for multilayer fabrics. The most relevant CFD simulations of fabrics are presented in Table 1.3. No simulations are reported in the literature to consider the flow through 3D models of multilayer forming fabrics.  5  Table 1.2 Screen/grids resistance-flow rate relationships from experiments Author Description By making an analogy to the flow through orifices and nozzles, Robertson  presented the resistance, discharge coefficient, over a range of  (1950) [23]  20 < Re d < 4000 . In contrast to a single orifice, at high flow rates the discharge coefficients could be greater than unity. By  Wieghardt (1953) [24]  means  of  experimental  measurements  over  the  range  60 < Re s < 6000 , he correlated the pressure drop coefficient, K , to Re s and β . Modeled the flow through screens with the flow through porous media.  Ingmanson et al. (1961) [25]  Using experimental values, they proposed a generalized equation in the form  −1  f = A Re s + B Re s  C  for 3 < Re s < 1500 .  Screen structural  properties are considered in the constants of the correlation. Armour and  Modeled screens as a thin packed bed and proposed the correlation  Cannon  f = A Re e + B for 0.1 < Re e < 1000 , applicable to diverse groups of  (1968) [26]  −1  screens. Compared the flow with orifice flow and presented a two-variable  Pedersen  discharge coefficient analytically, C d = f (ReW ) . He neglected the  (1969)  flow interaction between neighboring orifices and considered  [27]  cylindrical filaments in both the Machine Direction (MD) and Cross Machine Direction (CMD).  Wakeland and  Measured the low-frequency oscillating flow across square-mesh  Keolian  screens in the range 0.002 < Re d < 400 and correlated K , to Re d and  (2003) [28]  β.  6  Table 1.3 Fabric resistance to the flow through CFD simulations. Author Description and conclusions Lu et al. (1996) [34]  Extended the work of Pedersen for cylindrical filaments dissimilar in CMD and MD. They provided C d = A ReWB . Following the work of Pedersen derived a discharge coefficient,  Wang et al.  C d , for filaments of elliptical cross-section (a r = a b ) . They  (2007) [35]  provided  (  )  C d = Aa rB Ln ReW + Ca rD .  For  the  same  Re,  increasing a r decreased C d . He simplified fabrics as rows of dissimilar cylinders studied the Huang (2003) [36]  effect of second row on the single row flow field. He showed that if the offset between cylinder layers was greater than 0.7d, the second layer does not affect the flow upstream of the first layer. His work was experimentally validated in [29].  Green et al. (2008) [37]  Conducted studies on the flow through a single-layer fabric and showed if one filament was shifted laterally, the effect of the flow non-uniformity remained localized.  1.2.2 Literature Review on Fiber Modeling (Chapter 4)  Six degrees of freedom specify completely the position and orientation of a rigid particle in a flow. Analytical descriptions of the flow around particles are limited to the creeping flow regime. Some pioneering analytical work on the motion of a rigid slender particle is presented in Table 1.4. Flexible particles, in contrast to rigid ones, require infinite degrees of freedom to be specified spatially and temporally. There have been approaches to simulate flexible particles in fluid flows, which are tabulated in Table 1.5.  7  Author  Table 1.4 Slender particle in a creeping flow of a Newtonian fluid Description  Jefferey  Found the motion of a neutrally buoyant rigid elliptical particle in a v simple shear flow U ∞ = (γ&y,0,0) and showed the dimensionless period  (1922) [38]  of rotation, Tγ& , depends merely on the aspect ratio of ellipse Tγ& = 2π (rP + 1 rP ) .  Using the Euler’s equation they showed that the minimum value at Forgacs and Mason (1958) [39]  which rod-like particles of axis ratio rp and bending modulus Eb , may be expected to buckle under shear-induced compress is given by  (γ&μ )crit ≈ Eb (Ln 2rP − 1.75)  2rP4 .  Showed Jeffrey’s analysis to be valid for any axisymmetric body with Bretheron (1962) [40]  Tillet (1970) [41] Batchelor (1970) [42] Cox (1970) [43] Khayat and Cox(1989) [44]  fore-aft symmetry, provided that an equivalent aspect ratio re is used in place of the actual spheroid aspect ratio, rp . Computed the force analytically on long slender bodies in a general undisturbed creeping flow. The body’s cross-section was circular and center-line was straight. Modified the Tillett’s work for non-circular shape of the cross-section of the body. Cox’s theory allowed also a curvature of the body comparing to the work of Tillet and Batchelor. They improved the Cox’s theory including an expression in Re to include the inertia effects. The work corrected the earlier results of Chwang and Wu [45].  Subramanian  Studied the inertia effect in more detail and found a gentle O(Re) drift  (2005, 2006)  in closed Jeffrey orbits.  [46, 47]  8  Author  Table 1.5 Main methods to simulate flexible fibers in a flow field Description  Hinch (1976) [48] Gradon et al. (1988-1992) [49,50,51] Yamamoto and Matsuoka (1993-1997) [52-55]  Utilized the slender-body theory of Cox to formulate a flexible inextensible isolated thread in a shearing flow. Following the Hinch’s approach, they analyzed the deposition of a flexible fibrous particle in the creeping flow around a cylindrical obstacle. Any local contact of the particle with the cylinder surface would result in the particle collection by the cylinder. By inspiring from the molecular dynamics, they decomposed a fiber or plate-like particles into a lined-up spheres bonded to each neighbors. To ensure that the no-slip conditions are satisfied, an iterative constraint adjusting the sphere angles is solved to maintain the connectivity.  Lawryshyn and Presented a method to describe the dynamics of a flexible body based Kuhn (1997) [56]  on the equations of motion for a beam in vibration. They examined a flexible fiber motion in a channel flow with a slot. To eliminate the need for iterative constraint to maintain fiber connectivity in Yamamoto’s method, they proposed a particle-level  Klingenberg et al. (1997-2002) [57-60]  simulation method in which fibers are modeled as chains of spheroids/spheres/cylinders connected through ball and socket joints. To reduce the computation time, Skjetne’ method also neglected the hydrodynamic interactions between the spheres comprising a fiber. Schmid and Switzer’ work simplifies the calculation of interfiber separation (in suspensions) in Ross’s method and reduces the computation time.  Stockie and  They applied the Immersed Boundary Method (IBM), based on a  Green  mixture of Eulerian and Lagrangian variables [61], to fiber motion in a  (1998)  shear flow. Due to the two-way coupling between phases, they could  [62]  observe disturbance in the flow caused by the fiber motion.  9  Table 1.5 Main methods to simulate flexible particle in a flow field (cont.) Zhu used the IBM to simulate flapping filaments in a flowing soap Zhu and  film. His results showed the amplitude of flapping increased as the  Peskin  mass of the filament increased, needed for a sustained flapping. In the  (2002,2003)  case of two flapping fibers he found two distinct modes of sustained  [63,64]  synchronized flapping, depending on the ratio of the distance between the two fibers to their length [64].  Dong (2002) [65]  Dong used Ross’s approach in modeling the fiber moving through slots of different geometries in pressure screens.  Tornberg and  Based on a non-local slender body theory, they introduced a method to  Shelley  simulate interactions of flexible fibers in Stokes flow. They used the  (2004)  forces described by Euler-Bernoulli elasticity theory in lieu of the  [66] Qi (2006) [68]  slender body integral equations. He employed the original Lattice-Boltzmann method [67] for fiber simulations. This method can be used for low to finite Reynolds numbers flows considering inertial effects. Due to its computational time fibers with small aspect ratios were used in the simulations.  Lindstorm  In this recent work, Lindstorm extended Schmid’s approach to include  (2008)  the two-way coupling between the fiber and the fluid phase and particle  [69]  inertia was considered  These methods verified that they are able to reproduce the salient features of single fiber dynamics reported by Mason et al. [70,71], and the pole-vaulting behavior of a rigid fiber near a bounding surface, observed by Stover et al. [72]. There is no previous work to have studied flexible fiber motion in forming fabrics flow fields.  10  1.2.3 Literature Review on Circular Cylinders (Chapter 5)  In papermaking, pulp fibers/fines commonly move relative to the surrounding water at speeds up to U ∞ = 1.5 m s . The flow around them is approximately that around a circular cylinder of aspect ratio 2 − 100 with arbitrary yaw angle, and with a Reynolds number, Re = ρU ∞ D μ ( ρ = fluid density, μ = fluid dynamic viscosity), of up to 60.  The flow around a circular cylinder has been the subject of intense studies for a long time. Zdravkovich has collected most of the work [73,74]. In the creeping flow regime there exist analytical solutions for the flow of particles with different shapes [75,76]. For Re above unity, no analytical solution exists and therefore understanding of the flow depends on experimental data and numerical simulations. The flow around a cylinder is steady and laminar up to Re ≈ 40 , and thereafter becomes unsteady and laminar ( 40 < Re < 200 ) and then for Re > 300 the vortices shed in the wake become turbulent. For finite length cylinders, the onset of vortex shedding is found to be strongly dependent on the end configuration used [77-80], and is also a strong function of the boundary layer thickness on the walls to which the cylinder is attached [81]. Norberg has investigated the influence of aspect ratio on the Strouhal number for vortex shedding from a circular cylinder perpendicular to a flow, for a range of Reynolds numbers from about 50 to 4 × 10 4 [82]. He found that for aspect ratios in the range L D < 40 the onset of laminar vortex shedding was delayed. A summary on past work on infinite cylinders at normal and oblique angles to air stream at moderate to high Re is given in [83]. Norberg wrote a review of previous simulations of fluctuating lift and has found that most simulations have covered the range 45 < Re < 500 and 100 < Re < 10 5 for twodimensional and three-dimensional simulations, respectively [84,85]. In a recent work, some new types of periodic vortex shedding have been classified for cylinders of finite lengths at low Reynolds numbers 40 < Re < 300 [86]. When a circular cylinder is yawed to the flow, the cross section parallel to the flow becomes elliptical. The flow behavior around the yawed cylinder cannot be deduced from the flow around elliptical cylinders normal to the flow, as the streamlines bend as they approach the stagnation line, at least for higher velocity flows [87]. Ramberg  11  investigated the flow over yawed cylinders in the Reynolds-number range 160 < Re < 10 3 . He highlighted the sensitivity of the results to the cylinder end conditions [88]. Sears showed that the governing equations in a laminar boundary layer in the normal and spanwise directions are uncoupled, and thus the component of the force normal to a yawed cylinder may be calculated through “cosine laws”; the “Independence Principle” [89]. This simple approach has been used widely due to its practical application but it has a restricted range of validity. Thus, the empirical correlation of Roshko [90] for a cylinder normal to the flow has been extended to yawed cylinders by simply replacing the Reynolds number with the one defined based on the normal component of velocity [91,92]. Although there is a vast amount of information on finite and/or infinite cylinders either normal or oblique to the free stream velocity at high Re, surprisingly there is no previous research on yawed cylinders with finite aspect ratio in the Reynolds number range 1 < Re < 40 .  1.2 Problem Analysis  Figure 1.1 depicts a schematic of phases involved and the jet impingement in the dewatering process on the fabric of a Fourdrinier machine. The incoming jet impinges on the fabric with angles in the range 3o < α < 10 o , and the machine speed may reach up to 20 m s . The jet-to-wire speed ratio, JTW = V j VW , is normally between 0.990 and 1.035, and the pulp suspension consistency is typically 0.5% to 1%. In a typical machine, the drainage through the fabric may happen at speeds of U ∞ = 0.1 − 1.0m s [93, 94]. Thus, the Reynolds number of the flow, Re = U ∞ D υ ( υ is the kinematic viscosity), is in the range 15 to 150, calculated based on a fabric filament diameter on the order of D = 0.15mm . A typical pulp fiber is d = 30μm in diameter and L = 2mm in length.  Consequently, in the given cross flow, fiber Reynolds numbers Re f = U ∞ d υ are on the order of 3 to 30. The first stage of dewatering through fabrics may be conceptualized as follows: (1) free-jet zone in which the suspension jet travels between the headbox and the location where it impinges on the fabric, (2) oblique-dewatering zone where the jet has just 12  impinged on the fabric, (3) perpendicular-dewatering zone where the suspension drains in the Z-direction.  Free-jet zone  Oblique dewatering  Perpendicular dewatering free surface  V j = jet velocity  α VF =  fabric velocity impingement angle  Figure 1.1: Schematic of fabric and impinging jet on a Fourdrinier machine. It is instructive to look at time scales and length scales involved in the early dewatering process (upstream of paper-side filaments). The length scale of the fabric flow may be associated with the jet thickness, which is approximately 1cm. Therefore, the time a fluid particle spends near the fabric (early dewatering time) is on the order of 10 −2 s < t  fluid particle  =  1cm < 10 −1 s . To find particles time scales consider Figure 1.2, which U∞  shows a simplified model of a single-layer forming fabric as a row of cylinders. On the right of the figure, the different forces exerted on a particle are shown. For dilute suspensions the fiber/fiber interaction may be neglected. Moreover, the hydrodynamic interactions of neighboring fibers are also neglected. Using the creeping flow approximation for hydrodynamic forces, the normal and tangential components of the drag force on a slender particle with a semi-major axis a and semi-minor axis b are given [43] FD , parallel =  4πμa 8πμa ut = K t ut , FD ,normal = un = K nun ln (2a b ) − 0.5 ln (2a b ) + 0.5  13  FB  θ 0.1  m m < U ∞ < 0.1 s s  un  θ −α  1cm nˆ W  uf  ut  tˆ  Figure 1.2: (left) rows of circular cylinders as a simplified model of forming fabrics and relevant velocity and length scale, (right) different forces exerted on a fiber governing the fiber transient motion. The hydrodynamic interactions between particles are neglected in the force diagram. Applying the Newton’s Second Law to the particle shown in Figure 2.1 leads to the following equations of motion n-direction : W sin θ − FB sin θ − K n u n = ρ P ∀  du n dt  t-direction : W cos θ − FB cos θ − K t u t = ρ P ∀  du t dt  4 where ρ P ≈ 1100 kg m 3 is the particle density, ∀ = πb 2 a is the particle volume, FB is 3 the buoyancy force, and K n and K t are the drag coefficient in the normal and tangential direction, respectively. Solving the resulting first-order differential equation for the particle motion, the relaxation time of the particles (fibers, fines) would be τ = K ρ P ∀ . For a typical pulp fiber the relaxation time is 10 −4 s < τ fiber < 10 −3 s , and for a typical pulp fine it is 10 −5 s < τ fine < 10 −4 s . Thus, the Stokes number of rigid particles within such forming fabric flows are approximately in the range 10 −2 < Stk fiber = τ fiber t < 10 0 and 10 −3 < Stk fine = τ fine t < 10 −1 . These calculations show that pulp fines are expected to  follow the streamlines more closely than pulp fibers. Particles in the fabric flow fields may bend while rotating and translating. Therefore, their motion cannot be simply calculated as above and a more rigorous approach should be adopted. In the following sections, the description of flow field and flexible particle modeling are presented.  14  A comprehensive study of drainage phenomena and sheet formation in the forming section of a papermachine is a formidable challenge. Therefore, one should make reasonable assumptions before tackling the problem. In chapter 2 and 3, we investigate the pure water flow through fabrics associated with “perpendicular” and “oblique” dewatering. In chapter 4, we reduce the complexity of forming fabrics and consider only one strand of it. Therefore, flexible fiber motion in the flow field of one cylinder is taken as the representation of a more complicated problem. Chapter 5 deals with the drag and lift of finite cylinders inclined to flow at moderate Re, relevant to fiber motion in papermaking conditions. The summary and conclusions of the thesis is presented in Chapter 6.  1.3.1 Forming Fabrics Modeling (Chapters 2 and 3)  The primary objective of this research is to establish a relationship between fabric flow non-uniformity and uneven fiber distribution onto the fabric. The flow non-uniformity may occur on 3 different scales: the filament diameter, the weave pattern, and the edge flow. We are concerned here with non-uniformity on the scale of the filament and the weave pattern. It is known that this non-uniformity can have a profound effect on the quality of the printed end product [13, 14]. After CAD modeling of fabrics, a computational domain is constructed and the CFD software Fluent is then used to obtain the flows. The main assumptions in the fabric flow modeling are: 1) Fibers in the suspension are neglected; therefore we simulate a pure water flow through fabrics. 2) The free surface interaction between the suspension phase and the surrounding air is neglected. 3) Owing to the low Reynolds number of the flow, the laminar viscous model and steady solver are chosen. In cases where the Reynolds number is higher, we ran the simulations with an unsteady solver as well. The discrepancy between the results of the two solvers was less than 2%.  15  The impingement angle on fabrics in the fabric frame of reference, θ , can be related to that of the stationary frame, α , as tan θ =  JTW × sin α , where JTW is the jet-toJTW × cos α − 1  wire speed ratio. For simulations of oblique dewatering, the inlet volume is sloped with an angle θ with slip boundary condition applied on its walls. For simulations with inlet flow perpendicular to the fabric, both slip wall and periodic wall boundary conditions have been used; the difference in computed quantities between these boundary conditions was less than 2%. We either choose velocity/velocity or pressure/velocity boundary conditions for the inlet/outlet boundaries. The incoming free jet is turbulent, but in our modeling we neglect turbulent fluctuations. A typical computational domain for inclined simulations is shown in Figure 1.3.  Figure 1.3: Computational geometry indicating the inlet, yarn and outlet volumes.  16  When the flow fields are obtained, a representative plane a quarter of filament diameter away from the paper-side layer is considered. The pointwise velocities and average velocities are calculated on this plane. It is then assumed that the area-to-area fiber/fine distribution is proportional to flow non-uniformity associated with those areas. Therefore, initial particle mass deposition on the fabric can be inferred from the flow fields. Fabrics pressure drop-flow rate are also presented in a form C p = a Re b , where C p is the dimensionless pressure drop across the fabric and Re is the flow Reynolds number.  1.3.2 Fiber Modeling (Chapter 4)  Fibers in papermaking vary in length and flexibility. In chapter 4, we perform particle level simulations, proposed originally by Ross and Klingenberg [57], to model flexible fibers in fabric flows as well as their interaction with a model fabric. We briefly indicate the methodology employed here and more details of the method can be found in [57]. Figure 1.4 shows a flexible fiber in a general XYZ coordinate system in both v continuous and discrete form. The differential element of length ds in the continuous form is replaced by a segment in the discrete one. The center of mass of segment i is r denoted by vector ri , and the local coordinate system attached to the segment is xi y i z i . Each segment is a prolate spheroid with a major axis 2a and a minor axis 2b, connecting to its neighboring segments through ball and socket joints. The torsional springs at joints, depending on fiber flexibility, tend to bring a bent fiber back to its equilibrium position. i=N  zi yi  ds xi  i =1 Z  Z  Y  Y X  r ri  X  Figure 1.4: Representation of a fiber as a whole in continuous form (left) and as a chain of N linked rigid bodies in discrete form (right)  17  Isolating segment i from the whole chain, identifying forces and torques exerted on it, and solving the Newton’s second law and the law of moment of momentum give the evolution of the segment in the flow. The connectivity condition at joints provides the external equations to close the dynamics of the problem. The main assumptions in the fibre modeling are 1) Fibers segments are rigid, and the flexibility properties of the fiber are concentrated only at joints. 2) Fibers are inextensible. This can result in system indeterminacy. 3) Fibers are sufficiently thin that they do not disturb the background flow. Therefore, fibers move in a prescribed flow field. 4) All external forces pass through each segments’ center of mass. Therefore, the moment of tangential component of contact forces around the segment center of mass is neglected. 5) A linear dependence of drag and velocity is used for hydrodynamics forces. 6) Short-range and long-range hydrodynamics interactions between segments are neglected. 7) Fiber inertia is also neglected. 8) When fiber is in contact with a wall, unknown contact forces should be found simultaneously with velocities from the equations of motion. Considering the above assumptions, and rearranging the terms in the equation of motion, the resulting equations for translational and rotational velocities of segment i are r r r r&i = M ij1ω j + M 2 F jc + M 3 r r r N il1 ωl = N 2j F jc + N 3Ti c + N i4 The matrices  M ij1 , M 2 , M 3 , N il1 , N 2j , N 3 , N i4  (1) (2)  combines the dimensions, location,  flexibility, and orientation of segments with the flow properties ( U i∞ , Ω i∞ , E i∞ ) at the r r current location of the segment. F jc refers to contact forces and Ti c to contact torque if the fiber is fixed to a wall. The solution algorithm for a fiber in an unbounded flow is different than when it is in contact with a wall.  18  In an unbounded flow, there are no contact forces and torques. This decouples the rotational and translational equations of motion. Therefore, Eq.2 is first solved independently for angular velocities. The translational velocities are then found from Eq.1. When a fiber is in contact with a wall, the segments’ velocities and wall contact forces should be calculated simultaneously. Two cases arise here: (1) Fiber staples to the wall and remains stationary, (1) Fiber slides off the wall. The stationary condition requires that the two reaction forces are independent and r r only need to satisfy ∑ Fi t < μ f ∑ Fi n , where μ f is the coefficient of friction and the summation operator is performed over all segments in contact with the wall. The extra constraints to be solved along side the equation of motion for contacting fibers are r ui = 0 . This condition assures all components of velocity are brought to rest. If a contacting segment does not satisfy the stationary condition, then the fiber should slide over the cylinder. In this case, the normal and tangential components of r r contact forces are related by Fi t = μ f Fi n . The extra constraint to be solved with the r equation of motion is only u i .nˆ i = 0 ; segments should not be pass through the wall.  The fiber model developed here was compared against different benchmarks: (I) a rigid fiber motion in an unbounded shear flow, (II) a clamped fiber in a channel flow, and (III) a straight fiber sitting on a planar surface exposed to a uniform flow. One of the shortcomings of the fiber model manifested itself when we compared our results with the Capstan equation. We were able to match the fiber model simulations with the Capstab equation for low values of wrap angles or coefficients of friction, comparisons were more acceptable than higher values. Owing to the inextensible chain model, the system dynamics may not necessarily match with the geometrical constraints, resulting in system indeterminacy. We employ the particle level simulation to predict fiber motion in a fabric flow field. We consider only one filament of the fabric as a circular cylinder and simulate the fiber based on the proposed method which incorporates fiber/cylinder interaction. This enables us to calculate the contact forces at each time step and subsequently to define whether the fiber staples to the cylinder or slides off. The trapping region is presented  19  based on several dimensionless groups including fiber properties, flow characteristics, fiber/filament coefficient of friction, and the initial orientation and position of fiber.  1.3.3 Flow Around a Circular Cylinder (Chapter 5)  From a papermaking viewpoint, pulp fiber may be approximated by a circular cylinder moving in water with various inclinations. As already motioned, in the range of papermaking interest, 1 < Re < 40 , the flow around a finite circular cylinder inclined to the flow has not been studied before. In chapter 5 we document CFD simulations of the flow around finite cylinders with 2 ≤ L D ≤ 20 oriented to the flow at angles between 0 ≤ α ≤ 90  for  1 ≤ Re ≤ 40 . We present the  C L (α , L D , Re ) ,  C D (α , L D , Re )  relationships as best curve fits to computational data. The code validation is presented by comparing the simulations results with existing numerical and experimental data. The simulations are also used to test the "Independence Principle”, which assumes the drag coefficient for a cylinder is dependent only on the normal component of the velocity. Figure 1.5 shows the generic computational domain used for the simulations.  Figure 1.5: Computational domain indicating the inlet, outlet, and slip walls  20  1.4 Thesis Objectives  This thesis is intended to obtain some insights into: A. Forming Fabric 1) The flow micro non-uniformity caused by the forming fabrics structure in the initial stage of water removal from pulp suspensions. 2) Fiber/fine uneven distribution due to uneven fabric flows on a filament-based scale and a weave-based scale. 3) The influence of machine-side layer of forming fabrics on the upstream flow of the paper-side layer and on pressure drop-flow rate characteristics of fabrics. 4) The effect of speed difference between the oncoming jet and the machine on the pressure drop through the fabric and the shear stresses in the flow. B. Fiber Motion 5) A single fiber motion in the flow of one straightened strand of fabrics, a circular cylinder. This results in the classification of fiber/cylinder interaction as hung-up and not-hung-up regions in terms of several dimensionless groups. C. Finite Circular Cylinder Flow 6) Drag and lift of a finite cylinder 1 < L D ≤ 20 with an arbitrary angle to flow, ( 0 o ≤ α ≤ 90 o ) at low to moderate Reynolds numbers, 1 ≤ Re ≤ 40 .  21  1.5 References  [1] K.J. Haunreiter. 200th Anniversary of the paper machine, the second hundred years. Tappi Journal, 80(9):86-97, 1997.  [2] A. Roshanzamir. Hydrodynamics of Blade Gap Former. PhD Thesis, Department of Mechanical Engineering, The University of British Columbia, 2000. [3] F. Reddiough and A. Heinen. Sheet formation and dewatering in the Single Fourdrinier machine, Technical Assistance, Service and Know-how, 2001. [4] D.C. Miller. Practical aspects of water management (the operation and design of wet end drainage elements). Tappi, Wet End operations, short course, 1996. [5] S. Adanur. Effect of forming fabric structural parameters on sheet properties. Tappi Journal, 77(10):187-195, 1994.  [6] S. Adanur. Paper Machine Clothing. Asten, Inc, Basel, Switzerland, 1997. [7] T. Helle. Paper forming wires over 75 years. Pulp and Paper Canada. 91(6):107-114, 1990. [8] D.B. Johnson. Retention and drainage of forming fabrics. Pulp and Paper Canada. 85(6):126-132, 1984: [9] D.B. Johnson. Retention and drainage of multi-layer fabrics. Pulp and Paper Canada, 87(5):56-59, 1986. [10] R. Danby. The impact of multilayer fabrics on sheet formation and wire mark. Pulp and Paper Canada, 87(6):69-73, 1986.  [11] T. Helle. Influence of wire structure on sheet forming, Paper Technology and Industry, 123-131, 1980.  [12] T. Helle. Analysis of wire mark in printing paper. Journal of Pulp and Paper Science. 14(4): J91-J95, 1988.  [13] C. Antoine. Wire marking and its effect upon print-through perception of newsprint. Appita Journal. 60(3): 196-200, 2007.  [14] T. Helle. How forming fabric design affects drainage and release. Pulp and Paper Canada. 79(11):91-98, 1978.  [15] R.L. Beran. The evaluation and selection of forming fabrics. Tappi Journal, 62(4): 39-44, 1979.  22  [16] R.J.N. Hermer, G.H. Covey, W.D. Raverty, N. Vanderhoek, and G. Tan. Forming fabrics. drainage rates and paper properties. Appita Journal, 59(3):202-206, 2006. [17] J. Xu, R. Danby, D. Johnson, and J. VandeKolk. The effect of center plane resistance on the drainage property of forming fabric. 95th EXFOR and Annual meeting, 343-3.49, 2009. [18] R. Danby. The impact of forming fabric structures on print quality. Pulp and Paper Canada, 95(1): 48-51, 1994.  [19] R. Danby and P. Plouffe. Print quality improvements through forming fabric design changes. Pulp and Paper Canada, 101(9):66-69, 2000. [20] G.B. Schubauer, W.G. Spangenberg and P.S. Klebanoff. Aerodynamic Characteristics of Damping Screens. NACA TN 2001, Jan. 1950 [21] P. Bradshaw. The Effect of Wind-Tunnel Screens on Nominally Two-Dimensional Boundary Layers. Journal of Fluid Mechanics, 22(4):679-687, 1965. [22] J. Scheiman. Comparison of Experimental and Theoretical Turbulence Reduction Characteristics for Screens, Honeycomb, and Honeycomb-Screen Combinations. NASA Technical Paper 1958, 1981. [23] A.F. Robertson. Air porosity of open-weave fabrics, Part I: Metallic meshes. Textile Research Journal. 20(12):838-844, 1950.  [24] K.E.G. Wieghardt. On the Resistance of Screens. The Aeronautical Quarterly, 4: 186-192, 1953. [25] W.L. Ingmanson, S.T. Han, H.D. Wilder and W.T. Myers. Resistance of wire screen to flow of water. Tappi Journal, 44(1):47-54, 1961. [26] J.C. Armour and J.N. Cannon. Fluid flow through woven screens. American Institute of Chemical Engineers, 14(3):415-420, 1968.  [27] G.C. Pedersen. Fluid flow through monofilament fabrics. Paper presented at 64th national meeting of AIChE, New Orlean, 1969.  [28] R.S. Wakeland and R.M. Keolian. Measurements of resistance of individual squaremesh screens to oscillating flow at low and intermediate Reynolds numbers. Journal of Fluid Engineering, 125: 851-862, 2003.  [29] S. Chilchrist. Experimental investigation of a model forming fabric. Master Thesis, Department of Mechanical Engineering, The University of British Columbia, 2006.  23  [30] A. Rushton and P. Griffiths. Fluid flow in monofilament filter media. Transactions of the Institution of Chemical Engineers. 49:49-59, 1971.  [31] T. Miyagi. Viscous flow at low Reynolds numbers past an infinite row of equal circular cylinders. Journal of the physical society of Japan. 13(5):493-496, 1958. [32] J. Happel. Viscous flow relative to arrays of cylinders. American Institute of Chemical Engineers, 5(2):174-177,1959.  [33] V.C.L. Hutson.The resistance of gauze filters to slowly flowing fluid. Chemical Engineering Journal, 232:362-364, 1969.  [34] W.M. Lu, K.L. Tung, and K.J. Hwang. Fluid flow through basic weaves of monofilament filter cloth. Textile Research Journal, 66(5), 311-323, 1996. [35] Q.Wang, B.Maze, H.V. Tafreshi, and B.Pourdeyhimi. On the pressure drop modeling of monofilament-woven fabrics. Chemical Engineering Science, 62: 48174821, 2007 [36] Z. Huang. Numerical Simulations of flow through Model Paper Machine Forming Fabric. Master Thesis, Department of Mechanical Engineering, The University of British  Columbia, 2003. [37] S.I. Green, Z. Wang, T.Waung, and A. Vakil. Simulation of the flow through woven fabrics. Computers and Fluids, 37:1148 -1156, 2008. [38] G.B. Jeffery. The motion of Ellipsoidal Particles Immersed in a viscous Fluid. Proc. Roy. Soc. London Ser. A, 102:161-179, 1922.  [39] O.L. Forgacs and S.G. Mason. Particle suspension in sheared suspensions, IV. Spin and deformation of threadlike particles. Journal of colloid science, 14:457-472. 1959. [40] F.P. Bretheron. The motion of rigid particles in a shear flow at low Reynolds number. Journal of Fluid Mechanics, 14:284-304, 1962. [41] J.P.K. Tillett. Axial and transverse Stokes flow past slender axisymmetric bodies. Journal of Fluid Mechanics 44:401-417, 1970  [42] G.K. Batchelor. Slender-body theory for particles of arbitrary cross-section in Stokes flow. Journal of Fluid Mechanics 44:419-440, 1970. [43] R.G. Cox. The motion of long slender bodies in a viscous fluid, part 1, general theory. Journal of Fluid Mechanics, 44:791-810, 1970.  24  [44] RE. Khayat and RG. Cox. Inertia effects on the motion of long slender bodies. Journal of Fluid Mechanics. 209:435-462, 1989.  [45] A.T. Chwang and T.Y. Wu. Hydromechanics of low-Reynolds-number flow. Part 4. Translation of spheroids. Journal of Fluid Mechanics, 75: 677-689, 1976. [46] G. Subramanian and D.L. Koch. Inertial Effects on fiber motion in simple shear flow. Journal of Fluid Mechanics, 535:383-414, 2005. [47] G. Subramanian and D.L. Koch. Inertial Effects on orientation of nearly spherical particles in simple shear flow. Journal of Fluid Mechanics, 557:257-296, 2006. [48] E.J. Hinch. The distortion of a flexible inextensible thread in a shearing flow. Journal of Fluid Mechanics, 74:317-333, 1976.  [49] L. Gradon, P. Grzybowksi and W. Pilacinski. Analysis of motion and deposition of fibrous particle on single filter element. Chemical Engineering Science, 43(6): 12631259, 1988. [50] L. Gradon, A. Podgorski. Flexible fibrous particle behaviour in the carries gas flow around cylindrical obstacle. Chemical Engineering Science, 45(12):3435-3441, 1990. [51] P. Grzybowski and L. Gradon. Analysis of motion and deposition of fibrous aerosol particle flowing around a single filter element: The electrostatics effect. Chemical Engineering Science, 47(6):1453-1459, 1992.  [52] S. Yamamoto and T. Matsuoka. A method for dynamic simulation of rigid and flexible fibers in a flow field. Journal of Chemical Physics, 98(1):644-650, 1993. [53] S. Yamamoto and T. Matsuoka. Viscosity of dilute suspensions of rodlike particles: A Numerical simulation method. Journal of Chemical Physics, 100(4):3317-3324, 1994. [54] S. Yamamoto and T. Matsuoka. Dynamic Simulation of fiber suspension in shear flow. Journal of Chemical Physics, 102(5):2254-2260, 1995. [55] S. Yamamoto and T. Matsuoka. Dynamic Simulation of a platelike particle dispersed system. Journal of Chemical Physics, 107(8):3300-3308, 1997. [56] Y.A. Lawryshyn. Statics and Dynamics of Pulp fibers. Ph.D. Thesis, University of Toronto, 1997. [57] R.F. Ross amd D.J. Klingenberg. Dynamic Simulation of flexible fibers composed of linked rigid bodies. Journal of Chemical Physics, 106(7):2949-2960, 1997.  25  [58] P. Skjetne, R.F. Ross and D.J. Klingenberg. Simulation of single fiber dynamics. Journal of Chemical Physics. 107(6):2949-2960, 1997.  [59] C.F. Schmid, L.H. Switzer, D.J. Klingenberg. Simulations of fiber flocculation: Effects of fiber properties and interfiber friction. Journal of Rheology, 44(4):781-809, 2000. [60] L.H. Switzer. Simulating systems of Flexible Fibers. Ph.D. Thesis, University of Wisconsin-Madison, 2002. [61] C.S. Peskin. The Immersed Boundary Method, Acta Numerica, Cambridge University Press, 2002. [62] J.M. Stockie and S.I. Green. Simulating the Motion of Flexible Pulp Fibers Using the immersed Boundary Method. Journal of computational Physics, 147:147-165, 1998. [63] L. Zhu and C.S. Peskin. Simulation of a flapping flexible filament in a flowing soap film by the Immersed Boundary Method. Journal of Computational Physics, 179:452468, 2002. [64] L. Zhu and C.S. Peskin. Interaction of two flapping filaments in a flowing soap film. Physics of Fluids, 15(7):1654-1960, 2003.  [65] S. Dong. Modeling of fiber motion in pulp and paper equipment. Ph.D. Thesis, University of British Columbia, 2002. [66] A.K. Tornberg and M.J. Shelley. Simulating the dynamics and interactions of flexible fibers in Stokes flows. Journal of Computational Physics, 196:8-40, 2004. [67] S. Chen and G. Doolen. Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech, 30:329-364, 1998.  [68] D. Qi. Direct simulation of flexible cylindrical fiber suspension in finite Reynolds number flows. Journal of Chemical Physics, 125:114901, 2006. [69] S.B. Lindstrom. Modeling and simulation of paper structure development. Ph.D. Thesis, Mid Sweden University, 2008. [70] O.L. Forgacs and A.A. Robertson and S.G. Mason. The Hydrodynamic Behavior of Paper-Making Fibers. Pulp and Paper Magazine of Canada. 41(5):117-128, 1958. [71] O.L. Forgacs and S.G. Mason. The Flexibility of Wood-Pulp Fibers. Tappi, 41(11):695-704, 1958.  26  [72] C.A. Stover and D.L. Koch and C. Cohen. Observations of fiber orientation in simple shear flow of semi-dilute suspensions. Journal of Fluid Mechanics, 238:277296,1992. [73] M.M. Zdravkovich. Flow Around Circular Cylinders, Vol 1: Fundamentals. Oxford University Press, New York, 1997. [74] Zdravkovich MM . Flow Around Circular Cylinders, Vol 2: Applications. Oxford University Press, New York, 2002. [75] J.Happel and H.Brenner. Low Reynolds number hydrodynamics : With special applications to particulate media. Noordhoff International Publishing , Leiden, 1973.  [76] S. Kim and S.J. Karrila. Microhydrodynamics: Principles and Selected Applications. Dover Publications, New York, 2005. [77] F.H. Shair, A.S. Grove, E. Peterson and A. Acrivos. The effect of confining walls on the stability of the steady wake behind a circular cylinder. Journal of Fluid Mechanics, 17: 546-550, 1965. [78] S. Szepessy and P.W. Bearman. Aspect ratio and end plate effects on vortex shedding from a circular cylinder. Journal of Fluid Mechanics, 234:191-217, 1992. [79] A. Slaouti, J.H. Gerrard. An Experimental investigation of end effects on the wake of a circular cylinder towed through water at low Reynolds numbers. Journal of Fluid Mechanics, 112:297-314, 198.  [80] R. Stager and H.Eckelmann. The effect of endplates on the shedding frequency of circular cylinders in the irregular range. Physics of Fluids A, 3(9):2116-2121, 1991. [81] D. Sumner and J.L. Heseltine. Tip vortex structure for a circular cylinder with a free end. Journal of Wind Engineering and Industrial Aerodynamics. 96:1185-1196, 2008. [82] C. Norberg. An experimental investigation of the flow around a circular cylinder: influence of aspect ratio. Journal of Fluid Mechanics. 258:287-346, 1994. [83] J.D. Ju and R.L. Shambaugh. Air drag on fine filaments at oblique and normal angles to the air stream. Polymer Engineering and Science, 34(12):958-964,1994. [84] C. Norberg. Flow around a circular cylinder: aspect of fluctuating lift. Journal of Fluids and Structures, 15:459-469, 2001.  [85] C. Norberg. Fluctuating lift on a circular cylinder: review and new measurements. Journal of Fluids and Structures, 17:57-96, 2003.  27  [86] Q. Inoue and A. Sakuragi. Vortex shedding from a circular cylinder of finite length at low Reynolds numbers. Physics of Fluids, 20:033601-033601-12. 2008. [87] M. Shirakashi, Y. Ishida and S. Wakiya. Higher Velocity Resonance of Circular Cylinder in Cross Flow. Transactions of ASME: Journal of Fluids Engineering. 107:392– 396, 1985. [88] S.E. Ramberg. The Effect of yaw and finite length upon the vortex wakes of stationary and vibrating circular cylinders. Journal of Fluid Mechanics. 128:81-107, 1983. [89] W.R. Sears. The Boundary layer of yawed cylinder. Journal of Aerospace Science. 15:49-52, 1948. [90] A. Roshko. On the development of turbulent flow wakes from vortex streets, NACA report 1191, 1954. [91] C.W. Van Atta. Experiments on vortex shedding form yawed circular cylinders. American Institute of Aeronautics and Astronautics Journal, 6:931-933, 1968.  [92] A.R. Hanson. Vortex shedding from yawed cylinders. American Institute of Aeronautics and Astronautics Journal, 4:738-740, 1966.  [93] B. Dalpke, R.J. Kerekes, S.I. Green. Modeling Jet Impingement and the Initial Drainage Zone in Roll Forming. Journal of Pulp and Paper Science, 30(3):65-70, 2004. [94] A. Roshanzamir, S.I. Green and R.J. Kerekes. Two dimensional simulation of pressure pulses in blade gap formers. Journal of Pulp Paper Science, 24(11):364-368, 1998.  28  Chapter 2 Three-Dimensional Geometry and Flow Field Modeling of Forming Fabrics1 2.1 Introduction  In the process of making paper, a very dilute suspension of pulp fibers is spread over a forming fabric. While the forming fabric moves along the machine water is drained mechanically from the suspension and the deposited fibers form a wet sheet of paper on the grid [1]. The main requirements of a forming fabric are low running resistance and good fiber retention. Having a single layer fabric can enhance the water removal from the suspension, but at the cost of reduced retention. In contrast, a multi-layer fabric exhibits superior fiber retention, but an inferior drainage rate [2,3]. The mass distribution of the fibers in the sheet, the “formation” or macro-scale non-uniformity, is highly affected by the weave structure and number of layers [4]. It is postulated that the flow non-uniformity produced by the structure of a fabric can unevenly distribute the pulp fibers, producing hydrodynamic wire mark. For a normal fabric, flow non-uniformity may occur on 3 different scales: the filament diameter, the weave pattern, and the edge flow. We are concerned here with non-uniformity on the scale of the filament and the weave pattern. It is known that this non-uniformity can have a profound effect on the quality of the printed end product [4,5,6]. Uneven ink penetration resulting from the wire mark can be detected as peaks in the frequency spectrum obtained by FFT [7,8]. The fiber orientation distribution in a sheet, which in part results from the fabric geometry, has a tremendous effect on paper properties [9]. The fiber orientation distribution in a sheet arises from two sources: the fiber orientation in the headbox and fiber redistribution that occurs on and after jet impingement on the fabrics. Therefore, 1  A version of this chapter has been published. A. Vakil, A. Olyaei, and S.I. Green. Three-Dimensional Geometry and Flow Field Modeling of Forming Fabrics. Nordic Pulp and Paper Research Journal, 24:338346, 2009.  29  fabric characteristics such as air permeability, drainage capacity, and fiber support are highly connected with paper properties such as strength, smoothness, tensile and tear ratio. For example, as the air permeability and drainage capacity increases, fibers can interact more closely, which potentially results in a stronger sheet [10]. Previous investigations of forming fabrics considered pure water flow through simplified 2D model forming fabrics [11,12]. In these investigations, the machine direction filaments were removed and the fabric was idealized as banks of dissimilar cylinders. The authors showed that if the offset between cylinder layers was greater than 0.7d (d being the first layer filament diameter), the second layer does not affect the flow upstream of the first layer. Following these 2D investigations, studies on the flow through single layer woven fabrics were conducted [13] to examine the enhancement in the flow non-uniformity when one filament is shifted laterally. The results showed that the effect of the non-uniformity is highly localized. Flow through adjacent unperturbed apertures between filaments was affected by only a few percent by the displacement of one filament. The subject of the present work is a continuation of the aforementioned previous studies on forming fabrics. Here, we first present a new method to produce accurate CAD modeling of forming fabrics. To validate our model we simulated the flow through real industrial fabrics. The surprisingly good agreement between the experimental data and our simulations confirms that the method is accurate. We then present a simple theory that describes the fiber motion just above the top (paper-side) surface of fabrics, on a plane d/4 upstream of the fabric. The organization of the article is as follows: the numerical methods are presented in section 2. They include a novel method to build up a 3D CAD model from 2D x-ray images of fabrics. The paper continues by describing the meshing and computational domains of CFD simulations and presents some code validation. The third section is devoted to comparisons between multilayer fabrics in the form of air permeability and flow non-uniformity on the upstream side, and this section is followed by conclusions.  30  2.2 Numerical Method  2.2.1 CAD Modeling of Fabrics  A technique has been developed to produce accurate models of forming fabrics (shown in Figure 2.1. The accuracy of this technique is mainly due to the precision of x-ray technology and an image analysis code that is written in MATLAB software. As the first step in model production, a forming fabric sample of area 5 × 5 mm is micro C-T scanned yielding 1000 cross sections of equal depths, i.e. each slice is 5 micrometers thick. Typical slices are shown in Figure 2.1. Each resulting image is a machine direction (MD) by Z-direction (ZD), slice through the fabric. These 2D images are input into the DATA VIEWER software, which automatically constructs a 3D model of the sample. While the surface of the model generated by DATA VIEWER is not smooth, and can not be used in CFD, it does enable us to digitally extract another 1000 cross sections perpendicular to the original view, i.e. in the cross machine direction (CMD)-ZD plane. Refer to Figure 3.1.  (b)  (a)  Figure 2.1: Modeling of a forming fabric sample. (a) Prototype of a fabric (b) CAD model of the fabric.  31  Figure 2.2: Typical x-ray cross section images (MD horizontal and ZD vertical) taken at varying depths along the CMD.  Figure 2.3: CMD-ZD image planes extracted from DATA VIEWER. The micro C-T scan is more accurate than the micromachining approach adopted by Adebar et al. [14], and therefore so is the resulting 3D CAD model. Each slice from the two sets of MD and CMD 1000 image planes is fed into a custom code that processes each set of data in two stages: image feature extraction and feature processing. These stages produce three dimensional curves describing each filament in the model. Image feature extraction is a combination of two techniques that detect geometrical details of filaments from cross sectional images. First, all cross sectional Xray pictures are thresholded in MATLAB, yielding crisp boundaries between the filaments (white) and the background (black), as seen in Figure 2.4. Then, a geometry detection algorithm based on boundary search and centroid computation is implemented 32  to detect round geometries with isolated and discrete boundaries. In this algorithm, a detected boundary is tested with a roundness factor defined by 4πA P 2 (where A is the area and P is the perimeter), as shown in Figure 2.4. If the metric exceeds a threshold value of 0.8, the object is qualified as a circular or an ellipsoid cross section of a filament. For such objects the centroid of the cross section will be computed and added to the collection of data so the filament can be traced through its cross sections. Machine and cross-machine filaments are detected and processed from the orthogonal image sets, and then assembled to complete the model.  MD filament  Roundness Metric = 0.93  Roundness Metric = 0.08  CMD filament  Figure 2.4: Identifying cross sections by the roundness metric. This first algorithm is insufficient to generate good CAD models because in many places crossing filaments are in contact and hence the cross section of a yarn is not isolated in the image. Therefore a second algorithm, based on detection of partial circular boundaries, was developed to tell the blended entities apart. Since a forming fabric is often constructed from several different diameters of yarn, the diameters of all yarn filaments are measured far from where they interact with other filaments. This information is used as explained below. White elements identified by the original image processing that fail the centroid test are tested for the local radius of curvature. The boundary of each such white element is fitted with a series of circles that are locally tangent to the boundary (Figure 2.5). When many successive points along the tangent share the same local tangent circle, and that tangent circle diameter is about equal to that of the filaments, those points are 33  assumed to be part of a filament. The combination of these two algorithms yields a complete mapping of the machine and cross-machine direction filaments in 3D space.  Figure 2.5: Radius of curvature is computed continuously along the border by tangent circles.  Figure 2.6: The concentration of identical tangent circles in the regions marked in blue indicates a detected yarn strand. Once all cross sections have been detected, each yarn is associated with a 3D curve connecting centroids of the cross sections. A data processing algorithm sorts the collection of detected cross sections yarn by yarn so that these curves can be constructed and processed individually. The grouping is done by examining the proximity of centroids in adjacent image slices. The smoothness of these curves plays an important role in the smoothness of CAD surfaces that are to be created. The next step in model construction is therefore to smooth the raw centroid data collected from the CT images. To do this smoothing we connect the centroids and project the unsmoothed curve onto the YZ and XZ planes, as shown in Figure 2.7a and Figure 2.8a, respectively.  34  (a)  (b)  (c)  Figure 2.7: Projection of a 3D curve from raw data onto YZ plane. (a) Trace of filament in YZ plane (unprocessed) (b) Smooth 2D spline in YZ plane (c) Filament model viewed in YZ plane.  35  (a)  (b)  (c)  Figure 2.8: Projection of a 3D curve from raw data onto XZ plane. (a) Trace of filament in XZ plane (unprocessed) (b) Smooth 2D spline in XZ plane (c) Filament model viewed in XZ plane. From each projection a two-dimensional spline under tension is generated with a carefully chosen tolerance. Then, along the Z (shared) axis of the two orthogonal planes, and at equal increments, X and Y values are extracted, as shown in Figure 2.7b and Figure 2.8b. These values form control points of a smooth curve in three-dimensional space. This curve describes the shape of a yarn, as in Figure 2.9.  36  Figure 2.9: Individual filament that was modeled in Figure 2.7. Finally, these curves, MD and CMD sets, are imported into the GAMBITTM meshing software. Circles and ellipsoids, detected from the CT image cross sections, are swept along the curves creating volumes with smooth surfaces, as shown in Figure 2.10.  Figure 2.10: Assembled complete model of a forming fabric.  37  2.2.2 Meshing and Boundary Conditions  In a thorough modeling of the flow field through forming fabrics one must take into account the fibers suspended in the water and the free surface interaction between the suspension phase and the surrounding air. In this article we neglect the water-air interface. Furthermore, we neglect fibers in the suspension and therefore we simulate a pure water flow through the forming fabrics. We also only consider dewatering perpendicular to the fabric, i.e., we neglect fabric “rush” or “drag”. To produce a 3D simulation of the fabric, the computational domain first needs to be created and meshed. Then, appropriate boundary conditions for the inlet, outlet and lateral walls should be applied. The inlet-volume and outlet-volume were chosen to have periodic walls. The slip boundary condition was applied to the yarn-volume side walls, but no-slip was applied at the yarn surfaces. Finally, pressure inlet and pressure outlet boundary conditions were used to define the inlet and outlet (Figure 2.11).  Figure 2.11: Computational geometry indicating the inlet, yarn and outlet volumes.  38  As the flow gradients near the fabric are greater than elsewhere, this volume requires fine mesh elements. A yarn volume is created for areas around the fabric and this volume is meshed with unstructured tetrahedral volume elements. To save computational time, the inlet-volume and outlet-volume are meshed with structured wedge volume elements. Figure 2.12 shows a cross section of the fabric and the hybrid mesh used for it. The meshed geometry was then imported into the CFD package FluentTM to create a CFD simulation of the flow through fabrics.  Figure 2.12: Hybrid computational mesh scheme. The Reynolds number of a flow, Re=VD/ν, is a function of the fluid velocity V, fabric filament diameter D, and the fluid kinematic viscosity ν. Fabric manufacturers commonly test fabrics by driving air through the fabric with an applied pressure differential of 125Pa. This corresponds with a bulk air speed of about 2.5m/s for typical fabrics. For a filament diameter of 0.15 mm (also typical), the Reynolds number is 21. This Reynolds number is typical of those found on operating paper machines. For example, in blade and roll forming sections, water is driven through the fabric at speeds of perhaps 0.1-1.0 m/s [15,16]. These speeds correspond with Reynolds numbers between 15 and 150. Re=21 is within this range. We have done simulations at Reynolds numbers up to 125 (Figure 2.17), and the qualitative trends are independent of Re over this range. Owing to the low Reynolds numbers of the fabric flow, and in view of the qualitative trends mentioned above, the steady laminar viscous model was used in FLUENT. For certain cases with Re > 40 , unsteady simulations were also done. These simulations produced results within 2% of the steady computations. 39  2.2.3 Mesh Independence  Any CFD simulation should be run with meshes of different density to assure that the results are mesh independent. We have conducted simulations for one fabric geometry by gradually refining the mesh. Figure 2.13 shows the average velocity upstream of the fabric as a function of the number of computational cells. Simulations are run for a given pressure drop of a 125Pa across the fabric. The average velocity on a mesh with 2.5 million cells differs by less than 0.05% from the velocity computed with a mesh with 3.5 million cells. Therefore, most of the simulations reported here were run on a grid of approximately 2.5 million cells.  Figure 2.13: Upstream average velocity against mesh density for a pressure drop of 125Pa.  2.2.4 Discretization Error  As the mesh size approaches zero, a CFD solution should produce a solution to the continuum PDEs. In order to quantitatively estimate the size of the discretization error, we have followed the procedure proposed in (Roache 1997). On a systematically refined grid, the grid refinement ratio ( r ), is first defined as the ratio of discrete spacings on a  40  coarse grid ( h2 ) to that of the fine grid ( h1 ). As integral quantities are used for error estimations, these representative cells are chosen as: ⎡1 h=⎢ ⎣N  1  3 (ΔVi )⎤⎥ ∑ i =1 ⎦ N  where ΔVi is the volume of the i th cell and N is the total number of cells used for the computations. Then, the difference between a fine-grid solution and a coarse-grid solution is normalized by the grid refinement factor, based on the order of convergence of the solution. In this way, the maximum numerical uncertainty in the fine-grid solution for the velocity through the fabric was found to be 0.05%, consistent with 0.05% difference computed in section 2.3.  2.2.5 Code Validation  The manufacturer of the fabrics provided us with their experimental measurements of fabric flow. In their experiments they supply air at a gauge pressure of 125Pa on one side of a forming fabric and measure the flow rate (in cubic feet per minute) through a 7cm ( 2.75′′ ) diameter portion of the fabric. For the different fabrics used in this study, the air permeability from the simulations is compared against the experimental data of the manufacturer, as tabulated in Table 2.1. For fabric 1 the air permeability is 13% lower than the experimental value, whereas for both fabrics 2 and 3 this discrepancy is reduced to less than 1%. The good agreement between simulation and experiments validates the CFD code and the way fabrics are modeled. These results, in addition, confirm that we have improved our fabric modeling: fabric 1 was modeled based on the method presented in (Adebar et al. 2007), but fabrics 2 and 3 were modeled using the methodology presented in section 2.1 of this article. The extremely close agreement between the pressure drops predicted with the new modeling method, and the experimental pressure drops, provides compelling evidence that the 3D CAD models are highly accurate.  41  Table 2.1 Air permeability comparison between simulations and experimental data. The applied pressure drop across the fabric is 125 Pa. Fabrics 2 and 3 were modeled with the new technique described here. cfm × 0.00508 = 1m s . Simulations (cfm) Experiments (cfm) Difference % “Integra” Fabric (1) 424 485 12.5 “Integra T” Fabric (2) 436 440 0.9 “Sample 469” Fabric (3) 354 352 0.6  2.3 Results and Discussion  This section compares different forming fabrics flows. We compare the air permeability of different fabrics considering their structures. We also use a data analysis method to study the preferential orientation of fibers based on average velocity on the top surface of fabrics. The two fabrics considered in this study, the Integra and Integra T, are manufactured by AstenJohnson Inc. and are multi-layer fabrics whose layers are bound together with “binding filaments”. In general, these filaments repeat every 3 or 4 sheds. Top views of the fabric structures are shown in Figure 2.14 and Figure 2.15.  Figure 2.14: Top view of Integra structure, fabric 1. The units of the axes are metres. The dark elliptical area represents a typical fiber 1mm in length and 40 μm in diameter.  42  Figure 2.15: Top view of Integra T structure, fabric 2. The units of the axes are metres. The dark elliptical area represents a typical fiber 1mm in length and 40 μm in diameter.  2.3.1 Velocity Contours  Figure 2.16 shows the ZD, MD, and CMD components of the fluid velocity through the two fabrics, at a height d/4 above the fabric. This velocity is computed for air with an applied pressure differential of 125Pa across the fabric. One might speculate that the spatial variation of velocity contours upstream of the top layer arises principally from the top layer structure. This is indeed the case. The more nearly “square” top layer of Integra manifests itself as “checkerboard” velocity contours, whereas the velocity contours of Integra T are more elongated. Referring to Figure 2.16 b.1, it might be surprising that the prominent vertical gap in the filaments, centred near the MD coordinate=0.002, does not manifest itself as a region of very high velocity. The reason for this surprising behaviour is that the two CMD-filaments in close proximity to MD=0.002 are slightly lower (in the Z-direction) than the other filaments in the section. As the extraction plane is taken from the highest point of the fabric and parallel to the fabric, the reduction of the velocity that results from the locally larger distance between the filaments and the extraction plane compensates for the increased velocity caused by the local gap in the fabric. Because 43  fines and fillers very nearly follow the flow (neglecting their interaction with fibres), the factor of 3 variation in local velocity implies that there are areas in the fabric over which the initial accumulation of fines and/or filler content can be 3 times higher than in adjacent areas. Once fibres, fines, and fillers start to accumulate on the fabric, the “healing effect” will reduce the magnitude of this variation.  a.1) ZD velocity contours  b.1) ZD velocity contours  a.2) MD velocity contours  b.2) MD velocity contours  44  a.3) CMD velocity contours  b.3) CMD velocity contours  Figure 2.16: Velocity contours of ZD, MD, and CMD velocities. Velocities in meters per second. (a) Integra (b) Integra T.  2.3.2 Pressure Drop Versus Flow Rate  Simulations were run for a number of different pressure drops and the dimensionless pressure drop was thus determined as a function of the Reynolds number. Here, Re = VD /ν where V is the superficial velocity of the fluid through the fabric, D is the  filament diameter for filaments on the paperside (= 0.15mm for fabric 1 and 0.125mm for fabric 2), and ν is the kinematic viscosity of the fluid. Power law best curve fits to the results are shown in Figure 2.17. In both cases, the coefficient of determination, R 2 , is very high, indicating the quality of the function fit. Furthermore, for both fabrics C P is approximately proportional to Re −0.4 . Interestingly, the exponent on the Reynolds number is almost exactly half way between the creeping flow limit ( C P ∝ Re −1 ) and the fully developed turbulence limit ( C P ∝ Re 0 ). Since Δp ∝ C PV 2 and C P ∝ Re −0.4 and Re ∝ V , therefore Δp ∝ V  8  5  . In other words, if 125Pa of pressure difference drives  440cfm through an area of the Integra T fabric, 250Pa of pressure would be expected to ⎛ 250 ⎞ drive 440 × ⎜ ⎟ ⎝ 125 ⎠  5  8  = 616 cfm through the same fabric.  45  Figure 2.17: C P − Re curve for fabric 1 and fabric 2.  2.3.3 Velocity Averaging  The velocity contours displayed in section 2.3.1 are pointwise velocities. However, individual pulp fibers, and even large fines, have length scales comparable to, or larger than, many of the length scales that are salient in Figure 2.16. A pulp fiber will not respond to a local velocity, but rather to some average velocity over its length. To gain insight into this averaging, the average velocity on a circle above the fabric has been computed as a function of the circle radius. All averaging was carried out over an area parallel to the fabric surface. For circles drawn above a fabric opening, the average velocity should be highest for smaller circles. As the radius of averaging increases the average velocity decreases to finally reach the overall mean velocity, as shown in Figure 2.18a. In contrast to what happens over an opening, the circular averaging over a knuckle starts from a minimum and rises until it reaches the overall mean velocity (Figure 2.18b). In both examples, the ZD velocity component is within 5% of the mean value for r > 0.4mm . In broad terms, this means that a disk of radius greater than 0.4mm would be  pulled onto the fabric at about an equal rate, irrespective of the lateral location of the disk. 46  To obtain the spatial variation of circular averaging, a disk of radius 0.2mm was considered and averaging was performed over the disk area as the disk was moved parallel to the fabric. Figure 2.19 displays the isolines of such averaging. The contours show there are some areas over which a particle with r < 0.2mm , with an arbitrary orientation to the MD direction, is sucked through the fabric up to 20% faster than its neighboring areas. To account for the effect of particle orientation in the averaging, we performed averaging over an ellipse with semi-major axes of (0.2, 0.4, 0.5, 0.6) mm and a semiminor axis of 0.02mm (i.e., roughly corresponding with (0.4, 0.8, 1, 1.6) mm long fibers 40μm in diameter). The results of the averaging are shown in Figure 2.20 for α = 0 o (aligned with MD) and in Figure 2.21 for α = 90 o (aligned with CMD). The contours of  α = 0 o , if compared with circular averaging, show stronger localized variations for smaller aspect ratios. This implies that if small isolated fibers are oriented in the MD, their accelerating motion towards the fabric depends more highly on their spatial location than would be true of isolated disks of comparable maximum size. The isolines of  α = 90 o show an even greater variation spatially, so small fibers that are oriented in the CMD would be even more affected by flow non-uniformity than fibers in the MD. As the aspect ratio of the fiber increases, the ratio of maximum to minimum averaged velocities decreases (Figure 2.22). Therefore, as anticipated, the velocity variations are smoothed out for longer fibers.  47  (a)  (b)  Figure 2.18: Circular averaging over an opening and a knuckle, on a plane d/4 in the upstream of the Integra T fabric. (a) above an opening, (b) above a knuckle.  48  Figure 2.19: Spatial variation of circular averaging over r = 0.2mm for Integra T fabric.  49  (a) semi-major axis = 0.2mm  (b) semi-major axis = 0.4mm  (c) semi-major axis = 0.5mm  (d) semi-major axis = 0.6mm  Figure 2.20: Spatial variation of ZD velocity when averaged over an ellipse with semimajor axes = (0.2,0.4,0.5,0.6)mm and semi-minor axis = 0.02mm . The ellipse is aligned with MD. Velocities in meters/second.  50  (a) semi-major axis = 0.2mm  (b) semi-major axis = 0.4mm  (c) semi-major axis = 0.5mm  (d) semi-major axis = 0.6mm  Figure 2.21: Spatial variation of ZD velocity when averaged over an ellipse with semimajor axes = (0.2,0.4,0.5,0.6)mm and semi-minor axis = 0.02mm . The ellipse is aligned with CMD.  51  (a)  (b)  Figure 2.22: Effect of the aspect ratio on the averaged maximum and minimum velocities. The ellipse is aligned with (a) MD (b) CMD.  52  2.4 Conclusions  A novel method to build accurate 3D models of forming fabrics is presented. The accuracy of the method mostly relies on the high resolution of 2D C-T scans of the fabrics and the MATLAB code we developed to detect the boundaries of the filaments’ cross sections. By means of a roundness metric, the code identifies the circular and elliptical cross sections amongst the objects in each slice. A three dimensional spline is fit through the different elliptical cross sections. Filament volumes are constructed by sweeping ellipses along each fitted spline. We ran flow simulations of pressure-driven flow through the 3D fabric models. The predicted flow rates of the simulations differed by less than 1% from experimental measurements thereof. This method for generating 3D models has therefore been verified to be accurate. The computed velocity fields through each fabric can be interpreted to yield insight into the behaviors of fibers in the flow. Specifically, if the velocity field is averaged spatially over areas comparable to the projected area of a fiber, the resulting averages are related to the amount of fiber that should accumulate in different regions of the fabric during the initial stages of dewatering. For a particular industrial fabric studied, the ratio of zenith to nadir averaged ZD velocity is about 3; up to 3 times more fiber will initially accumulate on one area of the fabric than on other areas in the forming section.  2.5 Acknowledgments  The authors appreciate the financial support of the Natural Science and Engineering Research Council (NSERC) and AstenJohnson Inc.  53  2.6 References  [1] S. Adanur. Effects of forming fabric structural parameter on sheet properties, Tappi  Journal. 77(10):187-195,1994. [2] S. Adanur. Paper Machine Clothing, Asten, Inc, Basel, Switzerland, 1997. [3] T.K. Adebar, S.I. Green, R.J. Kerekes and S.M. Gilchrist. Imaging the ThreeDimensional Structure of Forming Fabrics. Tappi Journal, 6(8):11-15, 2007. [4] C. Antoine. Wire marking and its effect upon print-though perception of newsprint,  Appita Journal, 60(3):196-199, 2007. [5] B. Dalpke, R.J. Kerekes and S.I. Green. Modelling Jet Impingement and the Initial Drainage Zone in Roll Forming, Journal of Pulp and Paper Science, 30(3):65-70, 2004. [6] R. Danby. The impact of multilayer fabrics on sheet formation and wire mark, Pulp  and Paper Canada, 87(6):69-74, 1986. [7] R. Danby. The impact of forming fabric structures on print quality, Pulp and Paper  Canada, 95(1):48-51, 1994. [8] R. Danby. Print quality improvements through forming fabric design changes, Pulp  and Paper Canada, 101(9):66-69, 2000. [9] S. Gilchrist. Experimental investigation of a model forming fabric, M.A.Sc Thesis, The University of British Columbia, 2006. [10] S.I. Green, Z. Wang, T. Waung and A Vakil. Simulation of the flow through woven fabrics, Computers and Fluids, 37(9):1148-1156, 2008. [11] T. Helle. How forming fabric design affects drainage and release, Pulp and Paper  Canada, (79)11: 91-98, 1978. [12] T. Helle. Analysis of wire mark in printing paper. Journal of Pulp and Paper  Science, 14(4):J91-J95. 1988. [13] Z. Huang. Numerical simulations of flow through model paper machine forming  fabrics, M.A.Sc Thesis, The University of British Columbia, 2003. [14] D.B. Johnson. Retention and drainage of forming fabrics, Pulp and Paper Canada. 85(6):126-132, 1984. [15] D.B. Johnson. Retention and drainage of multi-layer fabrics, Pulp and Paper  Canada, 87(5):56-59, 1986.  54  [16] P.J. Roache. Quantification of uncertainty in computational fluid dynamics. Annu.  Rev. Fluid. Mech. 29:123-160, 1997. [17] A. Roshanzamir, S.I. Green, and R.J. Kerekes. Two dimensional simulation of pressure pulses in blade gap formers. Journal of Pulp and Paper Science, 24(11);364368, 1998.  55  Chapter 3 The Influence of Machine-Side Filaments and Jet-to-Wire Speed Ratio on the Flow Through a Forming Fabric 2 3.1 Introduction  In the process of making paper, a very dilute suspension of pulp fibers is spread over a forming fabric to mechanically drain water from the suspension [1]. Good forming fabrics are characterized by a low running resistance and good fibers retention. Having a single layer fabric can enhance the water removal from the suspension, but at the cost of decreased retention. In contrast, a multi-layer fabric exhibits superior fiber retention, but an inferior drainage rate [2, 3]. The mass distribution of the fibers in the sheet, the “formation” or macro-scale non-uniformity, is affected by the fabric weave structure and number of layers [4]. The systematic variation in the distribution and orientation of the fibers, fines, and fillers in finished paper is known as wire mark [5]. It is postulated that the flow non-uniformity produced by the structure of a fabric can unevenly distribute the pulp fibers, producing a hydrodynamic wire mark. The wire mark can be distinguished from random variation by examining the frequency spectrum of the printed structure [5, 6]. The flow non-uniformity may occur on 3 different scales: the filament diameter, the weave pattern, and the edge flow. We are concerned here with the non-uniformity on the scale of the filament and the weave pattern.  It is known that this non-uniformity  can have a profound effect on the quality of the printed end product [7, 8]. The fiber orientation distribution in a sheet arises from two sources: the fiber orientation in the headbox and fiber redistribution that occurs on and after jet impingement on the fabrics. It is known that fiber orientation before and upon 2  A version of this chapter has been submitted for publication. Ali Vakil, Arash Olyaei, and Sheldon Green. The Influence of Machine-Side Filaments and Jet-to-Wire Speed Ratio on the Flow Through a Forming Fabric.  56  impingement has a tremendous effect on the final paper properties [9, 10]. Paper properties such as strength, smoothness, tensile and tear ratio are tightly connected with fabric characteristics such as air permeability, drainage capacity, and fiber support [11]. Due to the complexity of the flow and fiber/fabric interaction much previous work has been limited to pure water flow through 2D model forming fabrics. In these investigations the forming fabrics were modeled as banks of non-uniform cylinders [12, 13]. It was found that if the surface spacing between the two rows was greater than 0.7d (d being the paper-side filament diameter), the upstream flow field remained fairly unaffected by the presence of the second layer. 3D simulations of flow through single layer fabrics have also been done [14]. These authors investigated the effect of lateral displacement of filaments on the non-uniformity of the flow through the adjacent unperturbed apertures. In the companion article [15], we described a novel method to accurately model forming fabrics. Based on an analytical technique we predicted the particle distribution in the flow fields. The flow non-uniformity and its probable effect on particles were dealt with both on the scale of filaments and a larger-scale weaved-based non-uniformity. The subject of this paper is a continuation of the studies in [15] on forming fabrics. Old single-layer forming fabrics were supplanted by a new generation of multilayer fabrics to meet two major requirements: increasing mechanical stability and lifespan on the machine-side layer by means of coarse mesh size; and improving formation on the top layer by using finer mesh size [4]. Although a multilayer fabric improves fiber retention and reduces the flow non-uniformities, it generally leads to lower drainage rates. It has been postulated that machine-side fabric filaments may adversely affect the flow field on the paper-side filaments. As one part of the present study we investigate this effect. We simulate the flow through a real industrial fabric and also its machine-side and paper-side layers separately, and infer from the flow field differences the effect of machine-side on paper-side flow. The jet to wire speed ratio is also modeled in this study by inclining the inlet volume of the simulations, and the spatial distribution of the velocity is compared between normal and inclined incoming flow to the fabric.  57  The organization of the article is as follows: The numerical methods are presented in section 2, including meshing and boundary conditions for the simulations. The next section is devoted to comparisons between single layer and multilayer fabrics and normal/inclined incoming flow to a multilayer fabric. The paper closes with conclusions.  3.2 Numerical Method 3.2.1 Meshing and Boundary Conditions  In comprehensive modeling of the flow field through forming fabrics one must take into account the fibers suspended in the water and the free surface interaction between the suspension phase and the surrounding air. In this article we neglect the water-air interface. Furthermore, we neglect fibers in the suspension and therefore we simulate a pure water flow through forming fabrics, which owing to the low Reynolds number is laminar. Figure 3.1 shows these phases and a schematic of jet impingement on the fabric. The first stage of dewatering may be conceptualized as follows: (1) free-jet zone in which the suspension jet travels between the headbox and the location where it impinges on the fabric (2) oblique-dewatering zone where the jet has just impinged on the fabric (3) perpendicular-dewatering zone where the suspension drains in the Z-direction. In the previous article [15] we directed our attention only to the third region. Therefore, the computational domain was constructed to simulate a flow normal to the fabric. A comparison of spatial differences in velocities and shear stresses between perpendicular and oblique dewatering is one subject of this paper. Free-jet zone  Oblique dewatering  Perpendicular dewatering free surface  V j = jet velocity  α VF =  fabric velocity impingement angle  Figure 3.1: Schematic of fabric and impinging jet. 58  In papermaking, the incoming jet impinges on the fabric with angles in the range 3o < α < 10 o , as observed in the frame of reference of the papermachine. In terms of the fabric reference frame, this impingement angle, θ , is: tan θ = in which JTW =  Vj  VW  JTW × sin α JTW × cos α − 1  is the jet-to-wire speed ratio, which is normally between 0.990  and 1.035 [16]. Therefore, for an impingement angle of α = 5 o with a JTW = 1.035 , the corresponding angle in the fabric reference frame would be θ ≈ 70 o . The three-dimensional modeling of fabrics based on 2D X-ray micro-tomography was described thoroughly in [15]. In brief, a MATLAB code identifies circular and elliptical cross sections, either isolated or blended, amongst the objects in each image. The resulting data are smoothed by generating a three dimensional spline through them. Filament volumes are then constructed by sweeping ellipses along each fitted spline. When a 3D model of the fabrics is built, a “yarn volume” is created around it. This “yarn-volume” contains tetrahedral elements around the fabric. The yarn-volume is then extended downstream to create the outlet volume. Upstream of the yarn volume is an inlet volume. For simulations of oblique flow, the inlet volume is sloping. The inlet/outlet volumes are structurally meshed to reduce computational costs. Slip wall boundary conditions are applied to inlet walls. For simulations with inlet flow perpendicular to the fabric, both slip wall and periodic wall boundary conditions have been used; the difference in computed quantities between these boundary conditions was less than 2%. Velocity inlet and pressure outlet boundary conditions are chosen to define the inlet and outlet boundaries. Figure 3.2 shows a representative cross section of the hybrid mesh used for simulations. A typical computational domain for inclined simulations is shown in Figure 3.3. The general specification of Fabric A used in simulations is given in Table 3.1.  59  Table 3.1 General specification of Fabric A used in simulations. Mesh/knock (1/inch) Total  Paper side  Yarn Size (mm)  Air  Machine Permeability (cfm) Paper side side MD CD  146 × 186 73 × 93 73 × 62  440  0.13  0.13  Figure 3.2: Hybrid computational mesh scheme.  60  Machine side MD  CD  0.21  0.20  Figure 3.3: Computational geometry indicating the inlet, yarn and outlet volumes. The meshed geometry was then imported into the CFD package FluentTM to create a CFD simulation of the flow through fabrics. Owing to the low Reynolds numbers for the flow, generally less than 40 based of filament diameter and bulk velocity, the steady laminar viscous model is used to specify the flow. In cases where the Reynolds number is higher, we ran the simulations with an unsteady solver as well. The discrepancy between the results of the two solvers was less than 2%.  3.2.2 Mesh Independence  Any CFD simulation should be run with meshes of different density to ensure that the results are independent of the mesh. We have conducted simulations for one fabric geometry by gradually refining the mesh. Figure 3.4 shows the average velocity upstream of the fabric as a function of the number of computational cells. Simulations are run for a 61  given pressure drop of a 125Pa across the fabric. The average velocity on a mesh with 2.5 million cells differs by less than 0.05% from the velocity computed with a mesh with 3.5 million cells. Therefore, most of the simulations reported here were run on a grid of approximately 2.5 million cells.  Figure 3.4: Upstream average velocity against mesh density for a pressure drop of 125Pa. Flow is normal to the fabric.  3.2.3 Discretization Error  As the mesh size approaches zero, a CFD solution should produce a solution to the continuum PDEs. In order to quantitatively estimate the size of the discretization error, we have followed the procedure proposed in [17]. On a systematically refined grid, the grid refinement ratio ( r ), is first defined as the ratio of discrete spacings on coarse grid ( h2 ) to that of the fine grid ( h1 ). As integral quantities are used for error estimations, these representative cells are chosen as: ⎡1 h=⎢ ⎣N  1  ⎤3 ( ) V Δ ∑ i ⎥ i =1 ⎦ N  where ΔVi is the volume of the i th cell and N is the total number of cells used for the computations. Then, the difference between a fine-grid solution and a coarse-grid 62  solution is normalized by the grid refinement factor, based on the order of convergence of the solution. In this way, the maximum numerical uncertainty in the fine-grid solution for the velocity through the fabric was found to be 0.05%, consistent with 0.05% difference computed in the previous section.  3.3 Results and Discussion  3.3.1 Single Layer Versus Multilayer Fabrics  It is worthwhile to analyze in detail the effects that a machine-side layer produces on the flow characteristics of a single layer fabric. To achieve this, the CAD model of Fabric A was cut on a plane between the layers, and CFD simulations were run for each layer separately. A comparison between the flow fields are given in the following sections. Figure 3.5 depicts the top view of Fabric A. The constituent layers of the fabric are shown in Figure 3.6. The binding filaments may be regarded as a center layer between the paper side and machine side of the fabric [18], but here we consider them as part of both the top and bottom layers.  Figure 3.5: Top view of Fabric A  63  (b) machine-side layer  (a) paper-side layer  Figure 3.6: Top view of layers composing Fabric A.  3.3.2 Velocity Contours  Figure 3.7a shows the ZD, MD and CMD components of fluid velocity through the complete Fabric A. Figure 7b shows the same velocity components for its paper-side layer computed in isolation. The contours are extracted on a plane d/4 upstream away of the highest point of the paper-side layer (d = paper-side filament diameter). There is surprisingly little difference between the corresponding velocity contours. This finding implies that the paper-side layer plays a dominant role on the flow in the pulp suspension.  64  a.1) ZD velocity contours  b.1) ZD velocity contours  a.2) MD velocity contours  b.2) MD velocity contours  a.3) CMD velocity contours  b.3) CMD velocity contours  Figure 3.7: Pointwise velocity contours of ZD, MD, and CMD velocities. Velocities in meters per second. (a) full Fabric A (b) Paper-side layer of Fabric A.  65  3.3.3 Averaged Velocity Contours  Typical pulp fibers are 300 − 3000μm long, and should respond to some average velocity over their length rather than a local velocity. We followed the procedure described in [15] to perform velocity averaging. We performed elliptical averaging over an area comparable to the cross-sectional area of a fiber for both the full Fabric A and its paperside layer only, and then subtracted the resulting contours. The result of such averaging for an ellipse with semi-major axis of 0.2mm and semi-minor axes of 0.02mm, are shown in Figure 3.8 (ellipse major axis in MD) and Figure 3.9 (ellipse aligned with CMD). The difference between “paper-side only” and “full fabric” velocities is up to approximately ± 15 % of the mean velocity. In approximate terms, the difference means that had the  fabric been comprised of only paper side filaments, the initial fiber mass distribution would be changed by up to 15%, relative to the “full fabric” case. %  %  Figure 3.8: Spatial variation of the difference between ZD velocities when averaged over an ellipse with semi-major axis = 0.2mm and semi-minor axis = 0.02mm. The difference is normalized by the mean velocity. The ellipse is aligned with MD.  66  %  %  Figure 3.9: Spatial variation of the difference between ZD velocities when averaged over an ellipse with semi-major axis = 0.2mm and semi-minor axis = 0.02mm. The difference is normalized by the mean velocity. The ellipse is aligned with CMD. 3.3.4 Pressure Drop Versus Flow Rate  Simulations were run for a number of different inlet velocities for Fabric A and its layers separately. The dimensionless pressure drop was thus determined as a function of the Reynolds number. Here, Re = VD /ν where V is the bulk velocity of the fluid through the fabric, D is the filament diameter for filaments on the paperside (= 0.125mm ), and ν is the kinematic viscosity of the fluid. The C P − Re curves are well described by a power law, with high coefficient of determination, as shown in Figure 10. Surprisingly, the power on the Reynolds number for all curve fits is about -0.36. As seen in Table 3.2, the pressure drop through the entire forming fabric is very nearly equal to the sum of the pressure drops through the two layers. Note further that because extra surface area is exposed when one cuts through the binder filaments, one would anticipate that the sum of the two “half” fabric resistances would slightly exceed that of the full fabric. This effect alone may account for the majority of the (already small) differences seen in Table 3.2.  67  Figure 3.10: C P − Re curve for full Fabric A and its comprising layers. Table 3.2 Comparison of the pressure drop through a full fabric and its layers. Re 14 22 33 40 62 88 C P , paper  27  21  18  17  15  14  C P , machine  42  33  29  27  24  22  69  54  47  44  39  36  C P , both  64  52  44  42  37  34  % difference  7  2  2  2  5  3  ∑  CP  3.3.5 Oblique Versus Perpendicular Dewatering  In all previous simulations, the approach flow was taken normal to the fabrics, i.e. we neglected the speed ratio between the jet and the fabric. In this section we consider approach flows 58o ≤ θ ≤ 90 o , corresponding to impingement angles 3o ≤ α < 15o and 0.99 ≤ JTW ≤ 1.035 . In this case, slip walls are assumed for the inlet walls rather than  periodic walls, as were used in previous simulations. Moreover, simulations were run on a Fabric A ∗ , which has the same weave pattern as the Fabric A but with somewhat thinner strands.  68  The most marked characteristic of rush or drag is expected to be observed in the shear stresses. To illustrate this, we averaged the MD shear stress over elliptical areas aligned with the MD. The results of averaging were then normalized with the dynamic pressure, calculated based on the approach flow. Figure 3.11 shows the contours of the elliptically averaged, normalized shear stress for different approach flow angles. There is an apparent similarity between the contours, but the mean value of the normalized shear stress for the 58˚ case was 0.011 whereas for 90˚ impingement it was -0.0013. As expected, the shear stress increases as the incoming flow becomes more parallel to the fabric. For flow perpendicular to an isotropic medium, one would expect zero mean shear stress. As seen in Figure 3.5, the piece of fabric chosen for the flow analysis is far from isotropic, which explains why the mean shear stress for 90˚ impingement is non-zero.  a ) 58˚  a ) 70˚  c) 90˚  Figure 3.11: Spatial variation of non-dimensional shear stress when averaged over an ellipse with semi-major axis = 0.2mm and semi-minor axis = 0.02mm. The ellipse is aligned with MD. Owing to the unphysical imposed zero shear stress wall condition applied to MD=2.5mm and MD=1.4mm, only the region of the fabric away from these walls is shown. 3.3.6 Pressure Drop Versus Flow Rate  We ran simulations for different inlet velocities while changing the impingement angle. The calculated pressure drops across the fabric were converted into pressure coefficients ( C P = ΔP  1 ρV 2 , where V is the bulk velocity). The pressure coefficient is plotted 2  69  against the flow Reynolds number (based on paper-side filament size of 0.125mm). In Figure 3.12 it is evident that there is an insignificant influence of the impingement angle on the bulk flow behaviour. In other words, the fabric permeability varies only slightly between rush and drag conditions.  Figure 3.12: C P − Re curve for Fabric A with impingement angle as a parameter.  70  3.4 Conclusion  We ran simulations of the flow through 3D models of forming fabrics, constructed by the procedure described in [15]. The CFD simulations of the models were shown to predict the fabric permeability within 1% of experimental measurements. A specific fabric, Fabric A (and a variant thereof), was used here to study two quantities of real forming fabrics: the influence of machine-side filaments and the jet-to-wire speed ratio on the flow through a forming fabric. We found that the machine-side filaments have a small influence on the flow field upstream of the paper-side filaments, and consequently on the initial fiber mass distribution of pulp fibers. The pressure drop of a fabric was found to be approximately the sum of that through the paper and machine side layers. We therefore conclude that the machine-side layer plays only a minor role on the upstream flow although it does, of course, increase the overall resistance of the fabric to the flow. The rush or drag of the impingement jet were simulated by inclining the inlet volume of simulations. The air permeability varied only slightly from 58˚ to 90˚ impingement. Rush or drag had a more pronounced effect on the MD shear stress, which will tend to orient fibres in the machine direction.  3.5 Acknowledgement  The authors appreciate the financial support of the Natural Sciences and Engineering Research Council of Canada and AstenJohnson Inc.  71  3.6 References  [1] S. Adanur. Paper Machine Clothing, Asten, Inc, Basel, Switzerland, 1997. [2] D.B. Johnson. Retention and drainage of forming fabrics, Pulp and Paper Canada. 85(6):126-132, 1984. [3] D.B. Johnson. Retention and drainage of multi-layer fabrics, Pulp and Paper Canada. 87(5): 56-59,1986. [4] R. Danby. The impact of multilayer fabrics on sheer formation and wire mark, Pulp  and Paper Canada. 87(6):69-74,1986. [5] T. Helle. Analysis of wire mark in printing paper. Journal of Pulp and Paper Science. 14(4):J91-J95, 1988. [6] C. Antoine. Wire marking and its effect upon print-through perception of newsprint,  Appita Journal, 60(3):196-199, 2007. [7] R. Danby. The impact of forming fabric structures on print quality, Pulp and Paper  Canada. 95(1):48-51, 1994. [8] R. Danby. Print quality improvements through forming fabric design changes, Pulp  and Paper Canada. 101(9):66-69, 2000. [9] T. Helle. How forming fabric design affects drainage and release, Pulp and Paper  Canada. 79(11):91-98, 1978. [10] S. Adanur. Effects of forming fabric structural parameter on sheet properties, Tappi  Journal, 77(10):187-195, 1994. [11] T. Helle, How forming fabric design affects drainage and release, Pulp and Paper  Canada, 79(11):91-98, 1978. [12] Z. Huang. Numerical simulations of flow through model paper machine forming  fabrics, M.A.Sc Thesis, The University of British Columbia, 2003. [13] S. Gilchrist. Experimental investigation on a model forming fabric, M.A.Sc Thesis, The University of British Columbia, 2006. [14] S.I. Green and Z. Wang, T. Waung and A.Vakil. Simulation of the flow through woven fabrics, Computers and Fluids. 37(9):1148-1156, 2008. [15] A. Vakil, A. Olyaei, and S.I. Green. Three-Dimensional Geometry and Flow Field Modeling of Forming Fabrics. Accepted at Nordic Pulp and Paper Research Journal.  72  [16] A.R. Guesalaga, F.M. Cabezas, C. Gottschalk. Measurement science & technology, 8(5):574-580, 1997. [17] P.J. Roache. Quantification of uncertainty in computational fluid dynamics, Annu.  Rev. Fluid. Mech. 29:123-160, 1997. [18] J. Xu, R. Danby, D. Johnson and J. VandeKolk. The effect of center plane resistance on the drainage property of forming fabric, 95th EXFOR and Annual meeting, 343-349, 2009.  73  Chapter 4 Flexible Fiber Motion in the Flow Field of a Cylinder3 4.1 Introduction  In papermaking, a dilute suspension of fibers in water is pumped either onto a single forming fabric, or in the gap between two forming fabrics. These fabrics typically consist of woven filaments about 150 microns in diameter, with a pitch of 300-400 microns. When a pressure differential is applied between the pulp suspension and the exterior of the forming fabric, the pulp suspension begins to flow through the fabric. The majority of the fibers in the suspension, which are on the order of 1000 microns in length, are trapped by the fabric while the water flows through the fabric pores. The resulting thickened wet pulp mat is then pressed and dried to make paper. During the pressing and drying stages of papermaking, a phenomenon known as linting often occurs whereby individual fibers detach from the pulp mat. Linting is believed to be caused by the mechanical detachment of those fibers that had partially penetrated the forming fabric upstream in the papermachine (Figure 1)  Machine Direction fibre mat  fibre likely to lint  Forming fabric filaments (cross-machine direction)  Figure 4.1: Schematic of linting precursor in forming section. Cross section through pulp mat and forming fabric. 3  A version of this chapter will be submitted for publication. Ali Vakil and Sheldon Green. Flexible fiber motion in the flow field of a cylinder.  74  As implied by the foregoing description, an essential part of papermaking involves the interaction of fibers with the filaments of a forming fabric. Filaments within a forming fabric are three-dimensional, approximately sinusoidal shapes, which are formed when flexible strands of material are woven into a mesh. The number of geometric variables required to describe such a mesh is very large. Owing to the geometric complexity of real fabrics, it is instructive to consider a simpler canonical problem. First, we consider only a single filament of the fabric. Second, we straighten the filament so it is a cylinder of infinite aspect ratio (i.e., 2D). Finally, although in papermaking a fiber can be at an arbitrary orientation relative to the cylinder, here we consider only cases where the fiber lies in the plane perpendicular to the cylinder axis. The motion of a rigid elliptical particle immersed in a creeping Newtonian fluid was computed analytically by Jeffrey [1]. Jeffrey’s analysis was shown to be valid for any axisymmetric body with fore-aft symmetry, provided that an equivalent aspect ratio  re was used in place of the actual spheroid aspect ratio, rp [2]. The force and torque acting on long slender bodies in a general undisturbed creeping flow were described analytically by Cox [3]. Cox later improved the theory by including a term in the Reynolds number to account for inertia effects [4]. Many approaches have been used to simulate flexible particles in fluid flows. The original Lattice-Boltzmann method, based on microscopic models and mesoscopic kinetic equations, has been employed for fiber simulations [5]. Due to its computational time only fibers with small aspect ratios were used in these simulations. The immersed boundary method based on a mixture of Eulerian and Lagrangian variables, related to one another by the Dirac delta function, has also been applied to this class of problems [6,7]. Hinch [8] has utilized Cox’s slender-body theory to model a flexible inextensible isolated thread in a shearing flow. Gradon et al. analyzed the motion and deposition of a flexible fibrous particle in the flow around a cylindrical obstacle at low Re [9,10]. The fiber-surface interaction was not considered in their method. In other words, any local contact of the particle surface with the cylinder surface resulted in particle adhesion. Kuhn et al. [11] presented a method to describe the dynamics of a flexible body for a beam in vibration based on the equations of motion.  75  Yamamoto and Matsuoka proposed a method to simulate fibers as a chain of spheres that can stretch, bend, and deform [12]. To eliminate the need to maintain fiber connectivity through iterative constraint, Ross and Klingenberg [13] proposed a particlelevel simulation method in which fibers are modeled as chains of prolate spheroids connected through ball and socket joints. Their approach, which we have used, is briefly described in section 2.2 of this paper. Their method can be used to describe large-aspect ratio fibers using relatively few spheroids. Skjetne [14] modeled flexible fibers as rigid spheres connected by ball and socket joints. To reduce the computation time, his method also neglected the hydrodynamic interactions between the spheres comprising a fiber. This method verified that the linked, rigid bodies are able to reproduce the salient features of single fiber dynamics reported by Mason et al. [15], as well as the pole-vaulting behavior of a rigid fiber near a bounding surface, observed by Stover et al. [16]. Based on the particle-level simulation technique, Schmid and Switzer [17,18] constructed a method in which fibers are modeled as rigid cylinders with hemispherical end-caps, connected by ball and socket joints. This method simplifies the calculation of interfiber separation used in Ross’s method, and reduces the computation time. Dong also used Ross’s approach as a way of modeling fibers moving through slots of different geometries [19]. In Lindstorm’s recent work, Schmid’s approach was extended to include  a two-way  coupling between the fiber and the fluid phase, and particle inertia was considered [20]. In this paper we employ a particle level simulation similar to that of Ross in order to predict flexible fiber motion in the flow field of a single cylinder, which represents one strand of a forming fabric. To consider the fiber/filament interaction, we rearrange the terms in the fiber equation of motion, and thus calculated the contact forces at each time step. We subsequently check the stationarity condition to determine whether fibers stick to the cylinder (“stapling”), or whether they slide off. In the second section of this article we present the single phase water flow around a cylinder. The code is validated by comparing the drag coefficient per unit length of a cylinder normal to the flow. We then describe the equations of motion for a flexible fiber. The fiber model is validated for both rigid and flexible fibers. The fiber friction has also been validated by different benchmarks. In the third section, we explain the physical parameters governing the fiber/filament system, and reduce the parameters to 7  76  dimensionless groups. By creating simulations that involve different sets of dimensionless groups we are able to draw a demarcation line between regions in the 7dimensional space, where fibers are “stapled” to the cylinder, and where they slide off the cylinder. The article closes with a discussion of results and conclusions.  4.2 Computational Methods  In our simulations we assume that fibers are sufficiently thin to the point where they do not disturb the background flow field as they move. Therefore, all flow fields are first regarded as single phase water flow over a cylinder. The commercial CFD software, FLUENT was used to solve the flow field around the cylinder at low Reynolds numbers. The single phase velocity field, vorticity field and rate of strain tensor are exported and reformatted to provide the input for a separate program that tracks fibers. Then, isolated fibers are added to the flow field, where their motion and interaction with the cylinder is tracked. It should be reiterated that a one-way coupling occurs between phases, and that the fiber is moved passively within the flow. This section explains the methods that we have employed for each phase in this article.  4.2.1 Single Phase Flow  Flow over a cylinder has been the subject of extensive study over in the past century [21]. There exists analytical work describing the flow when Re << 1 [22], but for Re ≥ 1 , there is no analytical work that describes the flow in detail. Therefore, in anticipation of future studies where Re >> 1 , CFD was used to obtain the cylinder flow field. A structured CFD mesh was generated in GAMBIT, as shown in Figure 4.1b. Once the CFD simulations were preformed in the CFD package FLUENT, the desirable field variables were exported and reformatted for later use. Figure 4.2a shows the computational domain including the boundary conditions chosen. The lateral walls were considered as zero-shear walls to minimize their effect on the flow.  77  Lup  Ldown  W  (a)  (b)  Figure 4.2: (a) Computational domain indicating the boundary conditions. (b) Close-up view of structural mesh around the cylinder. 4.2.1.1 Flow Validation  To validate the flow code, we compared the drag coefficient of our simulations with the numerical and experimental value of a cylinder in an infinite stream. The best-curve fit to the numerical data is expressed as [23] C D ≈ 1.18 +  0.0004 Re D 6.8 1.96 + 0.5 − 0.89 Re D Re D 1 + 3.64 × 10 −7 Re 2D  We conducted our simulations in a channel to save computational effort. A uniformity of the x-component of velocity was attained, typically for Y > 10 D , whereas the ycomponent of velocity was zero only on the wall. This component of velocity decreases noticeably as the distance increases from the side wall. We ran simulations with several blockage ratios in the range 0.005 < λ < 0.1 . The C D − λ curve at Re = 0.5 is shown in Figure 4.3. The y-intercept of the linear curve-fit to the simulations is 16.66±0.1. For a cylinder in an unbounded flow at Re = 0.5 the C D value from numerical data in [23] is  78  16.55 and from experimental measurements in [24] is 16.10±0.3. Therefore, we can conclude that we are in acceptable agreement with both numerical and experimental data. On the basis of these results, we selected a computational domain size Lup = 10 D ,  Ldown = 20 D and W = 10 D for future work.  Figure 4.3: Effect of the blockage ratio on the drag coefficient of a cylinder at Re = 0.5. 4.2.1.2 Pure Water Flow Around a Cylinder  Figures 4.4a-f depict the pure water flow around a single cylinder at Re = 1, based on the cylinder diameter. Each graph is self explanatory in showing the penetration depth of the fields. The flow field is uniform in the upstream, whereas the velocity gradients are substantially large while approaching the cylinder surface. When a slender particle is introduced into this flow field, the strain rate field tends to orient the particle with its axes parallel to the principle axes of distortion of the surrounding fluid, while the vorticity tensor tends to make the particle adopt the same rotation as the surrounding fluid [22].  79  a  b  c  d e e Figure 4.4: Flow streamlines and field variables around a single cylinder at Re = 1. a) X-velocity contours [m/s], b) Y-velocity contours [m/s], c) velocity magnitude contours [m/s], d) pressure field [Pa], e) strain rate field [1/s], f) vorticity magnitude [1/s].  4.2.2 Fiber Modeling  In order to model a flexible particle, one needs to consider a system with infinite degrees of freedom. There are several approaches suggested in the literature to attack such a problem, which are briefly mentioned in Section 4.1. Following the original method proposed by Ross and Klingenberg, we model each fiber as a chain of rigid bodies, connected through ball and socket joints [13]. Figure 4.5 illustrates a flexible fiber represented as a continuum and its counterpart in discrete form. The center of mass of r each segment i is located at ri , measured from a fixed reference frame. Each segment is a prolate spheroid with a major axis 2a and a minor axis 2b , as illustrated in Figure 4.6. r The orientation vector of each segment is Pi .  80  i=N  zi  yi  ds xi  i =1 Z  Z  r ri  Y  Y X  X  Figure 4.5: Representation of a fiber as a series of N linked rigid bodies. 2a  r Pi  i  Joint  Z  Joint i + 1  2b  r ri  Y X  Figure 4.6: A prolate spheroid representing each segment of a fiber. The motion of each segment is determined by solving the translational and rotational equations derived from Newton’s second law, and the law of moment of momentum. The assumptions made in the fiber model are: (1) Fibers are sufficiently thin; they do not disturb the background flow. (2) Hydrodynamic interactions between fiber segments are neglected. (3) Hydrodynamic forces on fiber segments are calculated using creeping flow theory. (4) The total hydrodynamic force on a fiber is the summation of hydrodynamics forces on its comprising segments. (5) Fiber inertia is neglected. (6) All external forces pass through the center of mass of the segments. 81  The resultant external forces on each segment are hydrodynamics forces and contact forces. The latter exists only when a segment touches a solid wall. For segment i with r r the angular velocity ω and the translational velocity U , the hydrodynamic torque and force are given by [22]  (  Fi h = Aij U ∞j − U j  (  )  (4.1)  )  Ti h = Cij Ω ∞j − ω j + H ijkl E ∞jk  (4.2)  where U ∞ , Ω ∞ , E ∞ are the velocity, vorticity, and rate of strain, respectively, of the imposed flow field at the center of mass of spheroid i . The resistance tensors are defined by  [ = 8πμa [X  ]  Aij = 6πμa X A d i d j + Y A (δ ij − d i d j ) C ij  3  C  (4.3)  ]  d i d j + Y C (δ ij − d i d j )  H ijkl = 8πμ a 3Y H ε ijl d l d k  (4.4) (4.5)  in which X A , Y A , X C , Y C , Y H are only functions of eccentricity. The internal resistance torque in the joints, Y j is a bending torque given by  (  )  Y jb = − k b θ j − θ 0j nˆ PB  (4.6)  where k b = EI 2a is the bending constant, θ 0j the equilibrium value and nˆ PB the unit vector normal to the plane of bending. By some manipulation of the governing equations, the final forms of the translational and rotational equations of motion for spheroid i are r r r r&i = M ij1ω j + M 2 F jc + M 3 (4.7) r r r N il1 ωl = N 2j F jc + N 3Ti c + N i4 (4.8) The scalar matrices in the above equations are  (  ~ ~ M ij1 = −b ji + A−1 Ak b jk  )  (4.9) (4.10) (4.11)  −1  M =A M 3 = A−1 A jU ∞j 2  [  ) ]  (  ~ 3 ~ ~ ~ d ij A j blj + d ij A j A −1 Ak blk 2 ar 1 ~ N 2j = d ij I − A j 4ar 3  N il1 = δ ij Ci −  [  82  ]  (4.12) (4.13)  N3 =  N = 4 i  ~ 3d ij ar  2  [A U j  ∞ j  (  1 8ar 3  )]  (4.14) ∞  − Ai A−1 AkU k∞ + H j E j + C j Ω j +  1 8ar 3  (4.15)  N  where A = ∑ Ai and ar = a b . The tilde operator is used to perform a cross product i =1  operation on two vectors as matrix multiplications. I is the identity matrix and exists ~ only for segments in contact with a solid surface. The matrix d represents the fiber center of mass based on connectivity vectors, extended from the center of mass of ~ segments to the neighboring joints. The matrix b is the position of the center of the segments relative to the fiber center of mass. All of the parameters in the above equations are in dimensionless form. The length, time and force scales are chosen as 2b , πμb 3 k b and k b b .  4.2.2.1 Fiber in an Unbounded Field  r r When a fiber is in an unbounded flow, there are no contact forces, i.e. F c and T c are zero. So, the governing equations reduce to r r r&i = M ij1ω j + M 3 r N il1 ωl = N i4  (4.16) (4.17)  The rotational evolution of fiber can be directly calculated from Eq.4.17. By substitution r of ω into Eq.4.16, the translational motion is then found. By integration from these equations, the updated position and rotation of each segment can be determined.  4.2.2.2 Fiber in Contact with a Wall  When a fiber is in contact with a solid surface, defined by S ( x, y ) = 0 , the segments’ velocities and contact forces should be calculated simultaneously. Figure 4.7 shows segment i of a fiber which has come in contact with a cylinder surface. The contact unit vectors are given by  83  r ∇S nˆi = r ∇S  ,  tˆi .nˆ i = 0  (4.18)  nˆ i segment i  tˆi Y S ( x, y ) = 0  X  Figure 4.7: Schematic of a fiber in contact with a solid surface. The main constraint here is that the fiber cannot pass through the wall. This imposes a r constraint on the normal velocity of the segment that touches the wall, i.e. u i .nˆ i = 0 . The solution procedure is therefore different from what we presented in section 4.2.2.1. We first multiply the translational equation of motion of those segments in contact with the wall by their associated normal unit vectors and rewrite the systems of Eq.4.7 and Eq.4.8 in a matrix form as follows: r r r r r&i .nˆi = M ij1ω j + M 2 F jc + M 3 .nˆi = 0 ⎫ ⎛ω j ⎞ ⎧ j = 1: N ⎪ r r (4.19) ⎬ ⇒ [A]⎜⎜ c ⎟⎟ = [B], ⎨ 1 r 2 c 4 N ij ω j − N j F j = N i ⎩i = contacts ⎝ Fi ⎠ ⎪⎭ r r r The unknowns are now ω and F c ; and F c contains both the normal and tangential  (  )  component of the reaction forces. Two cases arise here: (1) Fiber staples to the wall and remains stationary; (2) Fiber slides off the wall. In the case of sliding, the components of the reaction force are related by the coefficient r r of friction i.e. F t = μ f F n . In contrast, if the fiber staples to the wall the two reaction forces are independent and should only satisfy the inequality  84  rt  ∑F  i  r < μ f ∑ Fi n , where  the summation operator is performed over all segments in contact with the wall. The procedure solution is as follows: r (1) First, assume that the fiber remains stationary. An extra constraint of u i .tˆi = 0 is  simultaneously solved with systems of Eq.4.19. This assures all components of velocity are zero. (2) Check the contact forces if the stationary condition is satisfied,  rt  ∑F  i  r < μ f ∑ Fi n .  (3) If it is satisfied at the given time step, the fiber remains stationary. (4) If the stationary condition is not satisfied, the fiber slides. Then, (4-1) Assume that the fiber slides towards the right of the contact point. r r (4-2) The components of contact forcers are related by F t = μ f F n (4-3) Find the velocities and forces from the systems of Eq.4.19. r (4-4) Check the velocity if the right-sliding condition is satisfied u i .tˆi < 0 .  (4-5) If the last condition is satisfied, the velocities and forces are acceptable; otherwise; (4-6) The fiber should slide off towards the left of the contact point. (4-7) Find new forces and velocities. (4-8) Finally, by integration of the velocities the updated position can be determined. The fiber model consists of rigid spheroids connected by perfect ball and socket joints, making the entire fiber inextensible. As a consequence, when more than two segments of a fiber are in contact with the wall, the spheroid equations of motion are statically indeterminate. We overcame this indeterminacy by breaking the system down into two subsystems, as shown in Figure 4.8. We then assume in each subsystem that the segment on the wall remains stationary, and that the rest of the segments may hang around these fixed segments, i.e. segment (i − 1) for the left leg and segment (i + 1) for the right leg. The remaining segments between the left and right legs are also assumed to remain stationary.  85  i  y  i+1  i-1  x 1  N  (a) y  y  i+1  i-1  x 1  N  (b)  (c)  Figure 4.8: A statically indeterminate stationary system which is split into two subsystems as shown in (b) and (c). In each subsystem, the segments in contact with the wall are assumed to be stationary, but the rest may hang around them. 4.2.2.3 Fiber Model Validation  To validate the fiber model, the results of the fiber code are compared against different benchmarks. For a rigid spheroid in a simple shear flow, the period of rotation is compared with the analytical solution proposed by Jeffery [1], who showed that the dimensionless period of rotation, Tγ& , solely depends on the fiber aspect ratio given by: ⎛ 1⎞ Tγ& = 2π ⎜⎜ rP + ⎟⎟ (4.20) rP ⎠ ⎝ We ran the flow field simulation in FLUENT to obtain the flow between two parallel  planes with translational periodic lateral walls. A fiber possessing different aspect ratios was then introduced into the flow. The center of the fiber was located at a point in the fluid where its translational velocity was zero. Figure 4.9 shows a comparison of the period of rotation between the fiber model simulations and the analytical solution by Jeffery. The generally good agreement between the two graphs confirms that the fiber code is able to simulate rigid spheroids properly. The systematically lower rotation period 86  values relative to Jeffery’s theory has been noted by others [13], who defined an equivalent aspect ratio to match the simulation to the theory.  Figure 4.9: Comparison between Jeffery’s orbital periods and the fiber model simulations in a simple shear flow. 10 segments are used for the simulations. Pulp fibers vary in length and stiffness, which allows them to assume a different shape while they move and rotate within the flow [25]. Thus, to evaluate the flexibility segment of the fiber code, we compared the fiber simulations with those of a beam supporting a transverse loading. The free end deflection of a beam with a uniform transverse loading w , shown in Figure 4.10, is δ = wL4 8EI . Figure 4.11 plots the ratio of free end  deflection from the small deflection theory to that predicted by simulations. As the number of spheroids increases, the simulations asymptotically approach the theoretical result.  87  y  y L  L  w  w δ  y ( x) = −  (  wx 2 2 x − 4 xL + 6 L2 24 EI  x  x  )  (a)  (b)  Figure 4.10: A cantilever beam under a uniform transverse loading and its counterpart in the discrete form.  Figure 4.11: Ratio of the free end deflection from the small deflection theory to that predicted by the fiber model simulations.  Next, we estimated the deflection of a cantilevered fiber in a two-dimensional laminar channel flow (Figure 4.12) [20].  88  x y H  Um δ  Deflected fiber  Figure 4.12: A suspended fiber is a channel flow with a parabolic velocity profile.  Using the first approximation of Cox’s equations for tangential and normal forces on the fiber, the deflection of the free end with the given parabolic velocity profile is found to be [25]  δ=  μU m L5 2π (33H − 26 L ) 45 (EI )H 2 ln(D L )  (4.21)  Here, H is the width of the channel, L and D are the fiber length and diameter, μ is the viscosity of the fluid, and EI is the fiber stiffness. The same flow was computed and a fiber was placed in the simulated flow. A comparison between the fiber model simulations and those from Eq. 22 is presented in Table 1 for the specific case U m = 0.05 m s , H = 4mm , L = 2mm and D = 50μm .  10 spheroids were used to  represent the fiber in both simulations. The fibre end deflection is consistently overpredicted. A likely reason for this overprediction is that fibre segmentation in this model implies linearity. Consider fibers that are 40μm in diameter. The drag on two such attached fibers each of length 200μm will be exactly 2 times that of a single fiber 200μm long. In contrast, according to Cox’s theory, a fibre 400μm long (and 40μm in diameter) has only 1.6 times the drag of a 200μm fibre.  89  Table 4.1 Comparison between small deflection theory and the fiber model simulations of the free-end deflection of a cantilevered fiber in a channel flow. δ Eq.4.21 (m ) δ Code (m ) EI (N .m 2 ) Diff (%)  10 −12  2.45 × 10 −4  3.11 × 10 −4  27  10 × 10 −12  2.45 × 10 −5  3.19 × 10 −5  30  4.2.2.3 Fiber Friction  To validate the friction part of the fiber code, we considered the fiber interaction with both planar and curvilinear surfaces. First, a straight fiber was considered on a planar surface over which a uniform flow velocity was present. We were able to predict the onset of fibre sliding velocity to within 10%. To test the fiber friction on a curved surface we compared the fiber model with the Capstan equation either in the classical or modified form [27]. The simple form of the Capstan equation provides the ratio of tension on the tauter side to the slacker side of a rope wrapping a cylinder, as a function of the coefficient of friction and the wrap angle (Figure 4.13). Though the equation is derived for a continuous medium, we simulate such a system with a finite number of segments in our discrete approach. If we assume that the fiber wrapped around the cylinder is stationary, owing to the inextensibility of the fiber we have a statically indeterminate system, as mentioned already. Generally speaking, we were able to match the fiber model simulations with the Capstan equation for low values of the coefficient of friction and the wrap angle. To remove indeterminacy at higher values of friction and wrap, the Capstan equation could be matched by slightly inclining the last fibre segment to the cylinder.  Δθ  T1  Δθ  T1  T2  T2  T2 T1 = exp( f c Δθ )  Figure 4.13: A belt winding around a circular object and its counterpart in discrete form.  90  4.3 Results and Discussion  Numerous independent variables, which are tabulated in Table 4.2, govern the interaction of a fiber moving past a cylinder. These variables include the fluid properties ( ρ , μ ), the flow speed ( U ∞ ), the cylinder diameter ( D ), the fiber properties ( L , d , EI ), the fiber initial offset and orientation ( χ , ϑ ), and the frictional property of the fiber/cylinder interface ( f c ). As shown in Table 4.3, these variables may be non-dimensionalized, yielding 7 Π –groups.  Table 4.2 Salient variables that affect the interaction of a slender flexible particle with a rigid cylinder in a flow. Quantity  Density  ρ  Viscosity  μ  Flow Properties  Free stream velocity  U∞  Obstruction Properties  Cylinder diameter  D  Fiber length  L  Fiber diameter  d  Fiber stiffness  EI  Fluid Properties  Particle Properties Contact Properties Spatial Orientation  Coefficient of friction  fc  Initial offset  χ  Initial orientation  ϑ  91  Units ⎡ kg ⎤ ⎢⎣ m 3 ⎥⎦ ⎡ kg ⎤ ⎢⎣ m.s ⎥⎦ ⎡m⎤ ⎢⎣ s ⎥⎦  [m] [m] [m]  [N .m ] 2  [1] [m] [1]  Table 4.3 Dimensionless parameters for a slender particle in the flow of a cylinder Dimensionless parameter Fluid & Flow & Obstruction Properties  Reynolds number Fiber aspect ratio  Particle Properties  Internal restoring torque to External viscous torque  Particle to Object Properties  Fiber length / Filament diameter  Contact Properties  Coefficient of friction  Spatial Distribution  Initial offset / Filament diameter Initial orientation  Π –groups Π1 =  ρU ∞ D μ  L d EI Π3 = μU ∞ L3  Π2 =  L D Π5 = fc  Π4 =  Π6 =  χ D  Π7 = θ  Figure 4.14: Typical configurations of a fiber in the proximity of a cylinder. The angle θ and offset χ are measured far upstream of the cylinder. Papermaking fibers are suspended in water ( ρ ≈ 1000 kg m 3 , μ ≈ 10 −3 N .s m 2 ). In the frame of reference of the forming fabric, the free stream velocity will be 0 < U ∞ < 1.5 m s . Fabric filament diameters are approximately D = 150μm . Typical  92  fiber lengths are 0.5mm < L < 3mm . Typical fiber diameters are 15μm < D < 40μm . It is shown in the work of Kerekes et al. [26] that the fiber flexural rigidity varies in the range of 10 −12 N .m 2 < EI < 10 −10 N .m 2 ; in recent work the lower limit of EI was found to be  10 −13 N .m 2 [28]. The coefficient of friction is 0.2 < f c < 0.6 [29]. The initial offset and orientation are arbitrary. Hereafter, the simulations are performed based on the mentioned range of variables. The motion of fibers within the flow and their interaction with the cylinder are then characterized based on the dimensionless groups. As examples of typical fiber motion in the vicinity of a cylinder, we present two figures. In Figure 4.15, the fiber is comparatively rigid (larger value of П3), and its motion may be described as follows: i) the fiber is carried with the flow towards the cylinder surface without bending, appreciably, ii) it impacts the wall and starts to rotate. For fibers located almost symmetrically and for small initial angles of the fiber, the rotation process is slower, iii) it continues pivoting around the contact point, with or without sliding, iv) it finally slides off the cylinder surface. Upon leaving the surface, due to the high viscous forces at the rear of cylinder, the fiber can be drifted slightly back towards the centerline of the cylinder (i.e. for a fiber leaving the cylinder on the left, this would involve a motion towards the right). It then continues in purely translational motion without further appreciable bending or rotation. Another phenomenon shown in Figure 4.15b is that of a fiber slanted downwards on the left, whose center is shifted to the right of the cylinder center, which rotates clockwise and slides off from the right of the cylinder. Obviously, for smaller rightward offsets, a left-slanted fiber would slide off to the left of the cylinder.  93  initial position Initial position  Figure 4.15: Motion of a nearly rigid fiber near a cylinder. The fiber is shown at unequal time increments. Π 1 = 1 , Π 2 = 40 , Π 3 = 1.25 , Π 4 = 8 , Π 5 = 0.2 , Π 7 = 5 o . (a) Π 6 = 0 (b) Π 6 = 1 . In Figure 4.16, the fiber flexibility plays a significant role. Ten spheroids are used to represent the fiber. In contrast with a rigid fiber, which reaches an equilibrium “hung-up” position on the cylinder only when the fiber is located symmetrically on the cylinder, flexible fibers may become trapped on the cylinder, even if they have been significantly offset from the cylinder center. As explained above, the motion of a fiber in the vicinity of a cylinder is a function of 7 dimensionless groups. Owing to the large number of Π groups that govern the motion, the complete state space is extremely rich and difficult to convey graphically. In the sections that follow, subsets of this space will be explored and described.  94  initial position Equilibrium position  initial position Equilibrium position  (a)  (b)  initial position Equilibrium position  (c)  Figure 4.16: Fiber hang-up on a cylinder for different values of flexibility with f c = 0.2 , and Re = 1 . The fiber stiffness in figures (a), (b), and (c) is respectively EI = [80,25,5]× 10 −12 N .m 2 . 4.3.1 Subset of Space State with Π 3 =  EI χ , Π5 = fc , Π 6 = as Variables 3 D μU ∞ L  Figure 4.17 shows neutral stability curves for fiber hang-up, with Π 6 = χ D ,  Π 3 = EI μU ∞ L3 , and Π 5 = f c as variables. In this figure, the remaining Π -groups were held constant: Π 1 = 1 , Π 2 = 40 , Π 4 = 8 , Π 7 = 0 o . For any particular value of Π 5 (one curve in the figure), fibers with configurations to the left of each curve became “hung-up” on the cylinder, and fibers to the right of each curve did not get “hung-up” on the cylinder. So, for example, a fiber with Π 6 = 0.2 and Π 3 = 0.1 will be hung-up on the cylinder if Π 5 = 0.4 but not if Π 5 = 0.2 . The behaviours shown in Figure 4.17 are physically sensible. For higher values of friction coefficient, fibers can be offset further and still get hung-up. These larger regions of hang-up are shown in Figure 4.17 as the shift in the graphs to the right. Owing to symmetry, all fibers with Π 7 = 0 will be hung up if Π 6 = 0 . Thus, all curves should asymptotically approach the y-axis, as they do. Highly flexible fibers (smaller values of  Π 3 ), will get hung-up for larger fiber offsets than will less flexible fibers. This is why all curves are concave upwards.  95  Not-hung-up  Hung-up  Figure 4.17: The hung-up regions of the space state as a function of Π 3 , Π 5 , and Π 6 . 4.3.2 Subset of Space State with Π 3 =  EI , Π 5 = f c , Π 7 = θ as Variables μU ∞ L3  In this case, we kept all dimensionless groups fixed except Π 3 , Π 5 and Π 7 . As illustrative examples, fiber motion and hang-up on the cylinder for different initial orientation and different stiffness are shown in Figure 4.18. The coefficient of friction is 0.2 and Re = 1 . Fiber and cylinder dimensions are similar to those considered in Figure 4.17. As fiber bending increases with reducing stiffness, flexible fibers, will become trapped by the cylinder at a higher initial orientation angle than a rigid ones.  (a)  (b)  (c)  Figure 4.18: Fiber hung-up on a cylinder for different values of flexibility with f c = 0.2 , and Re = 1 . The fiber stiffness from left to right is EI = [140,40,2]× 10 −12 N .m 2 , and θ = [13o ,16 o ,40 o ]  96  We ran the simulations for different coefficients of friction and, as expected, increasing the coefficient of friction caused the region of hang-up to increase (Figure 4.19). The shape of the curves in Figure 4.19 is similar to that in Figure 4.17 because these curves too must asymptotically approach the y-axis for the same reason (symmetry) as do the curves in Figure 4.17.  Π7 = θ o  Figure 4.19: The hung-up region for as a function of Π 3 , Π 5 , and Π 7 . It is also instructive to study the way that a fiber slides off of the cylinder. A rigid fiber with an initial orientation of 14˚ is shown in Figure 4.20a. It is carried by the upstream flow until impact with the cylinder. After its contact with the cylinder, it pivots around the contact point while sliding, and then continues in translational motion, without significant rotation. In contrast, the flexible fiber in Figure 4.20b follows the streamlines more closely. After its contact with the cylinder, it bends and then slides off. After detaching from the cylinder, it is dragged slightly towards the rear of the cylinder, and becomes more aligned with the flow when it reaches the rear of the cylinder. Generally, a flexible fiber follows the flow streamlines more than a rigid one. This tendency to more closely follow the flow can be attributed to the fact that flexible fibers possess lower effective Stokes number values than rigid fibers of the same overall length and diameter.  97  (a)  (b)  Figure 4.20: Snapshots of a fiber sliding on a cylinder (a) rigid fiber (b) flexible fiber. Flow Reynolds number is 1, and EI = [80,2]× 10 −12 N .m 2 from left to right. 4.3.3 Subset of Space State with Π 1 =  ρU ∞ D EI χ , Π3 = , and Π 6 = as 3 D μ μU ∞ L  Variables  To consider the effect of Reynolds number on the trapping behavior, we decreased the Reynolds number from 1 to 0.5. As shown in Figure 4.21, the effect of this reduction in Reynolds number was to increase (slightly) the amount of fibre offset at which fiber trapping could still be observed. The difference between the two Reynolds number cases may be attributable to the shorter range of viscous action at the higher Re. This shorter range implies higher local velocity at a fixed distance from the cylinder, which would cause offset fibres to experience greater drag.  98  Figure 4.21: The hung-up region as a function of Π 1 , Π 3 , and Π 6 .  4.3.4 Subset of Space State with Π 1 =  ρU ∞ D EI , Π3 = , and Π 7 = θ as μ μU ∞ L3  Variables  The shape of the curves in this subset of state space is similar to that described in 4.3.3, except that now the x-axis represents the initial orientation of the fiber. For this case too, increasing the Reynolds number creates a smaller trapping region.  Π7 = θ o Figure 4.22: The hung-up region as a function of Π 1 , Π 3 , and Π 7 .  99  4.3.4 Subset of Space State with, Π 2 =  EI L χ , Π3 = , and Π 6 = as Variables 3 D D μU ∞ L  In this section we consider the influence of the ratio of fiber length to cylinder diameter on hang-up. Figure 4.23 shows these trapping regions for an aspect ratio of 20 and 28. As expected, longer fibers possess a larger hung-up region. It is interesting to note that for a rigid fiber, the shift of the graph is not as significant, while a flexible fiber responds to the change in fiber length more drastically. The flexible fiber effect is readily understood. With L D = 28 , even with χ D = 0.9 , an approximately 14.1D length of fiber trails behind the cylinder on one side and a 12.3D length trails on the other side. The normal force acting on the cylinder, which is roughly proportional to the drag on the two exposed lengths of the fiber (26.4D) generates a frictional force sufficient to support the differential drag between the two sides (roughly proportional to 1.8D). The ratio between these two forces is roughly 15. In contrast, with L D = 20 and χ D = 0.9 , an approximately 10.1D length of fiber trails behind the cylinder on one side and a 8.3D length trails on the other side. In this case, the ratio of normal force to differential drag will be about 10. The rigid fiber case is more difficult to understand as the rigid fibers tend to slide off the cylinder in response to more subtle flow differences.  Figure 4.23: The hung-up region as a function of Π 2 , Π 3 , and Π 6 .  100  4.4 Summary and Conclusion  In the forming section of a papermachine, pulp slurry drains through a forming fabric, causing individual fibers in the slurry to be trapped by the fabric. The interaction of fibers amongst themselves, with the surrounding water, and in relation to the fabric, ultimately contributes to the final spatial orientation of fibers in the final fiber mat. Carrying out comprehensive modeling of such a multiphase flow would be highly complicated, and thus, as an initial step we simplified the process and developed a numerical scheme based on the following assumptions. First, we considered only a highly diluted pulp suspension, where the fiber/fiber interactions can be neglected. Therefore, only single fibers needed to be tracked. Furthermore, the three-dimensional geometry of the fabrics was replaced with a single straightened strand. We further assumed that the fiber motion was constrained to the plane perpendicular to the axis of the strand. These assumptions resulted in our simplification of the tremendously complicated multiphase flow problem into our canonical problem of “single fiber motion in the flow of a single cylinder.” Fibers are free to move, rotate, and bend within these flows. In representing flexible fibers, we adopted the particle level method listed in [13], and modified it to fit our purposes. Each fiber was represented by a chain of rigid spheroids, connected by a ball and socket joint. The torsional springs at the joints serve to give flexibility to the fibers. The interaction of an inextensible fiber with a cylinder is statically indeterminate. We resolved this indeterminacy by breaking the system down into two subsystems (one on either side of the cylinder centerline) and solved the equations of motion of the two systems simultaneously. In real forming fabric flows, the Reynolds number may reach up to 40. Unfortunately, nowhere in the literature is there information on the drag force acting on spheroids, with arbitrary angles to the flow, at Reynolds numbers between 1 and 40. Consequently, we employed the creeping flow approximation, which is only valid for Re less than about 1. We validated the fiber model simulations by comparing our results with different benchmarks. In general, our results matched the benchmarks reasonably well, but due to the limitations of our fiber model and the benchmarks, in some cases the difference between the simulations and the benchmarks was up to 30%. 101  We found seven dimensionless groups that govern the motion of the fiber within the flow, as well as the fiber’s interaction with the cylinder. They are Π 1 = ρU ∞ D μ (flow Reynolds number), Π 2 = L d (fiber length-diameter ratio), Π 3 = EI μU ∞ L3 (ratio of the fiber bending stiffness restoring force to the drag force), Π 4 = L D (ratio of fiber length to cylinder diameter), Π 5 = f c (friction coefficient between fiber and cylinder),  Π 6 = χ D (non-dimensional fiber offset relative to the cylinder center), and Π 7 = θ o (fiber initial angular orientation relative to the cylinder). We defined the hung-up region of the fiber based on tests involving different combinations of these groups. In general, increasing the fiber flexibility, the fiber length, and the coefficient of friction resulted in a higher likelihood of hang-up. Conversely, increasing the Reynolds number decreased the likelihood of hang-up. Many of these effects are explained in terms of the complex interplay between the hydrodynamics of the flow and the dynamics of the fiber/cylinder interaction.  102  4.5 References  [1] G.B. Jeffery. The motion of ellipsoidal particles immersed in a viscous Fluid. Proc.  Roy. Soc. London Ser. A, 102:161-179, 1922. [2] F.P. Bretheron. The motion of rigid particles in a shear flow at low Reynolds number.  Journal of Fluid Mechanics, 14:284-304, 1962. [3] R.G. Cox. The motion of long slender bodies in a viscous fluid, part 1, general theory.  Journal of Fluid Mechanics, 44:791-810, 1970. [4] R.E. Khayat, RG. Cox. Inertia effects on the motion of long slender bodies. Journal of  Fluid Mechanics, 209:435-462, 1989. [5] D. Qi. Direct simulation of flexible cylindrical fiber suspension in finite Reynolds number flows. Journal of Chemical Physics, 125:114901-10, 2006. [6] J.M. Stockie, S.I. Green. Simulating the motion of flexible pulp fibers using the immersed boundary method. Journal of Computational Physics, 147:147-165, 1998. [7] L. Zhu, C.S. Peskin. Simulation of a flapping flexible filament in a flowing soap film by the Immersed Boundary Method. Journal of Computational Physics, 179:452-468, 2002. [8] E.J. Hinch. The distortion of a flexible inextensible thread in a shearing flow. Journal  of Fluid Mechanics, 74:317-333, 1976. [9] L. Gradon, P. Grzybowksi, W. Pilacinski. Analysis of motion and deposition of fibrous particle on single filter element. Chemical Engineering Science, 43:1263-1259, 1988. [10] L. Gradon, A. Podgorski. Flexible fibrous particle behaviour in the carrier gas flow around cylindrical obstacle. Chemical Engineering Science, 45:3435-3441, 1990. [11] Y.A. Lawryshyn. Statics and dynamics of pulp fibers, Ph.D. Thesis, University of Toronto, 1997. [12] S. Yamamoto, T. Matsuoka. A method for dynamic simulation of rigid and flexible fibers in a flow field. J. Chem. Phys., 98:644-650, 1993. [13] R.F. Ross, D. Klingenberg. Dynamic simulation of flexible fibers composed of linked rigid bodies. J. Chem. Phys., 106:2949-2960, 1997. [14] P. Skjetne, R.F. Ross, D. Klingenberg. Simulation of single fiber dynamics. J. Chem.  Phys., 107:2949-2960, 1997. 103  [15] F. Folgar, C. Tucker. Orientation behavior of fibers in concentrated suspensions. Journal of Reinforced plastics and Composites, 4:98-119, 1984. [16] C.A. Stover, D.L. Koch, C. Cohen. Observations of fiber orientation in simple shear flow of semi-dilute suspensions. Journal of Fluid Mechanics, 238:277-296, 1992. [17] C.F. Schmid, L.H. Switzer, D.J. Klingenberg. Simulations of fiber flocculation: Effects of fiber properties and interfiber friction. J. Rheol, .44: 781-809, 2000. [18] L.H. Switzer. Simulating systems of flexible fibers, Ph.D. Thesis, University of Wisconsin-Madison, 2002. [19] S. Dong. Modeling of fiber motion in pulp and paper equipment. Ph.D. Thesis, University of British Columbia, 2002. [20] S.B. Lindstrom. Modeling and simulation of paper structure development. Ph.D. Thesis, Mid Sweden University, 2008. [21] M.M. Zdravkovich. Flow Around Circular Cylinders, Vol 1: Fundamentals. 1997;  Vol 2: Applications. 2003; Oxford University Press, New York. [22] S. Kim, S.J. Karrila. Microhydrodynamics, Principle and Selected Applications. Dover Publications, Inc. New York, 2005. [23] D. Sucker, H. Brauer. Fluiddynamik bei quer angestromten Zylindern. Warme-  Stoffubertragung, 8:149-158, 1975. [24] D.J. Tritton. Experiments on the flow past a circular cylinder at low Reynolds numbers. Journal of Fluid Mechanics, 6:547-567, 1959. [25] O.L. Forgacs, S.G. Mason. The flexibility of wood-pulp fibers. Tappi, 41:695-704, 1958. [26] P.A Tam Doo, R.J. Kerekes. A method to measure wet fiber flexibility. Tappi, 64:113-116, 1981. [27] J.H. Jung, N. Pan, T.J. Kang. Capstan equation including bending rigidity and nonlinear frictional behavior. Mechanism and Machine Theory, 43:661-675, 2008. [28] D. Yan, K. Li. Measurement of wet fiber flexibility by confocal laser scanning microscopy. Journal of Material Science, 43:2869-2878, 2008. [29] S.R. Andersson, A. Rasmuson A. Dry and wet friction of single pulp and synthetic fibres. Journal of pulp and paper science, 23:J5-J11.1997.  104  Chapter 5 Drag and Lift Coefficients of Inclined Finite Circular Cylinders at Moderate Reynolds Numbers4 5.1 Introduction  In papermaking, pulp fibers commonly move relative to the surrounding water at speeds up to U ∞ = 1.5 m s . Typical pulp fibers have a diameter of D = 20 − 40μm and are  300 − 3000μm long, and may be oriented in arbitrary directions relative to the flow. Pulp fines are cut pulp fibers, and are typically 20 − 40μm in diameter and 40 − 300μm in length. The flow around pulp fibers and fines is therefore approximately that around a circular cylinder of aspect ratio 2 − 100 with arbitrary yaw angle, and with a Reynolds number, Re =  ρU ∞ D ( ρ = fluid density, μ = fluid dynamic viscosity), of up to 60. μ  Analytical solutions exist for the creeping flow ( Re << 1 ) around finite cylinders with arbitrary orientation to the flow [2]. However, if Re is not much smaller than 1, no exact analytical solution for the flow exists, although approximate theories do exist [3, 4, 5]. As the Reynolds number increases above 1, the flow field about a long cylinder normal to the flow changes from a small wake, to a larger wake with steady regions of recirculating flow, to a laminar unsteady wake for Re > 47 [6]. The shed vortices may be either parallel or oblique to the cylinder axis [7]. For both infinite and finite length cylinders, the onset of vortex shedding is found to be strongly dependent on the end configuration used [8], and is also a strong function of the boundary layer thickness on walls to which the cylinder is attached [9]. At low Reynolds numbers, 60 < Re < 200 , Slaouti and Gerrard have studied the effect of various end constructions on the wake structure behind a circular cylinder with aspect ratio about 4  A version of this chapter has been published. A. Vakil and S.I. Green. (2009) Drag and Lift Coefficients of Inclined Finite Circular Cylinders at Moderate Reynolds Numbers. Computers and Fluids, 38:17711781.  105  30 [10]. They observed that vortices originating from the free end bend to spread along the span. These vortices are strongly affected by the constraints imposed by the end configuration. They postulated that oblique shedding may arise from a difference in the two end effects. Norberg has investigated the influence of aspect ratio on the Strouhal number for vortex shedding from a circular cylinder perpendicular to a flow, for a range of Reynolds numbers from about 50 to 4 × 10 4 [11]. He believes that the aspect ratio effects and /or different end conditions provide a convincing explanation for the data scatter reported in the literature. His measurements indicate that the onset of vortex shedding is more or less independent of the aspect ratio in the range from about L D = 40 to L D = 1000 , whereas for smaller aspect ratios the onset of laminar vortex shedding was delayed: L D = 100 ⇒ Re c = 47.5 , L D = 30 ⇒ Re c = 49.5 and L D = 20 ⇒ Re c = 51.5 . It is  shown  that  critical  Reynolds  numbers  fit  well  into  the  relation  −2  ⎛L⎞ Re c = 47.4 + 1.8 × 10 ⎜ ⎟ [11]. Szepessy and Bearman found that aspect ratio also has ⎝D⎠ 3  an effect on the shedding frequency. For example, St = 0.19 for Re < 4.5 × 10 4 and L D = 6.7 reduced to St = 0.17 at similar Reynolds numbers but L D = 1 [12]. Norberg  wrote a review of previous simulations of fluctuating lift and has found that most simulations have covered the range 45 < Re < 500 and 100 < Re < 10 5 for twodimensional and three-dimensional simulations, respectively [13,14]. In the recent direct numerical simulations by Inoue and Sakuragi [15], some new types of periodic vortex shedding and transition between different types have been classified for cylinders of finite lengths at low Reynolds numbers 40 < Re < 300 . The dependence of the critical Reynolds number and the Strouhal number on L D in their simulations were in close agreement with the previous work mentioned above. Schouveiler and Provansal used an LDV system to study unsteady wake formation behind cylinders in a wind tunnel with 1 ≤ L D ≤ 5 [16]. They also found that the critical Reynolds number decreases with  increasing aspect ratio. When a circular cylinder is yawed to the flow, the cross section parallel to the flow becomes elliptical. The flow behavior around the yawed cylinder cannot be deduced  106  from the flow around elliptical cylinders normal to the flow, as the streamlines bend as they approach the stagnation line, at least for higher velocity flows [17]. Ramberg investigated the flow over yawed cylinders in the Reynolds-number range 160 < Re < 10 3 . He highlighted the sensitivity of the results to the cylinder end conditions [18]. Sears showed that the governing equations in a laminar boundary layer in the normal and spanwise directions are uncoupled, and thus the component of the force normal to a yawed cylinder may be calculated through “cosine laws”; the “Independence Principle” [19]. This simple approach has been used widely due to its practical application but it has a restricted range of validity. Thus, the empirical correlation of Roshko [20] for a cylinder normal to the flow has been extended to yawed cylinders by simply replacing the Reynolds number with the one defined based on the normal component of velocity [21, 22]. Although there is abundant previous work on finite and/or infinite cylinders either normal or oblique to the free stream velocity, the authors could find no previous research on yawed cylinders with finite aspect ratio in the Reynolds number range 1 − 40 . This particular Reynolds number range is of interest because, as stated previously, it is the range of relevance to papermaking fibers in water. In this article we describe computer simulations of the lift and drag of finite aspect ratio cylinders oriented to the flow at angles between 0 ≤ α ≤ 90 for 1 ≤ Re ≤ 40 . The following section of this article describes the numerical methods, meshing and mesh independency checks and code validation. In the third section the drag and lift curves versus the angle of inclination for different Reynolds numbers and different aspect ratios are presented. Best curve fits to the data are obtained. Conclusions and future work round out the article.  107  5.2 Numerical Method 5.2.1 Meshing and Boundary Conditions  Cylinders with four different aspect ratios ( L D = 2, 5, 10, and 20), with angles of inclination 0 ≤ α ≤ 90 to the flow were considered. Each cylinder was placed in the domain shown in Figure 5.1. The domain dimensions varied slightly between cases, but its rough dimensions are as given in Figure 5.1. The domain extends from 10 L upstream to 15 L downstream of the cylinder, where L is the cylinder length. The lateral slip walls are extended far enough away from the cylinder to provide a uniform far field flow. As an example of the uniformity of the far field velocity, the computed velocity profile is shown in Figure 5.2 for the case L D = 10 , α = 50o , Re = 40 and in Figure 5.3 for the case L D = 10 , α = 50o , Re = 1 . As one approaches the lateral boundaries of the domain,  typically for Y > 2 L , the velocity falls within 0.1% of the free stream velocity. Note that the velocity profiles illustrate the dominance of the viscous forces for Re = 1 , and the  inertial forces for Re = 40 : near the cylinder the flow is accelerated for Re = 40 , whereas it is decelerated for Re = 1 .  Figure 5.1: Computational domain indicating the inlet, outlet, and slip walls.  108  Figure 5.2: Velocity profile along the dashed line at Re = 40 , showing the negligible influence of the wall U∞  Velocity extracted along the dashed line.  Figure 5.3: Velocity profile along the dashed line at Re = 1 , showing the negligible influence of the wall. To reduce the computational time a hybrid mesh is created for the domain, as shown in Figure 5.4. In the areas close to the cylinder, tetrahedral elements are used,  109  whereas wedge elements are used far from the cylinder. Figure 5.5 depicts a close-up view of the unstructured mesh around the cylinders.  U∞  9L  2L  14 L  Figure 5.4: Hybrid computational mesh used in this research.  Figure 5.5: Unstructured mesh used close to the cylinder. The meshed geometry was then imported into the CFD package FluentTM. Owing to the low to moderate Reynolds numbers of interest here, the laminar viscous model in Fluent was used, resulting in a direct numerical simulation of the flow. Velocity inlet and pressure outlet boundary conditions were used to define the inlet and outlet. Since  110  unsteady vortex shedding normally begins when Re ≥ 50 for cylinders of aspect ratio L D ≤ 20 , almost all simulations were run under steady flow conditions. We ran some  simulations both steady and unsteady for the case L D = 10 , α = 40 o , Re = 40 . The time-averaged lift and drag coefficients for this case are given in Table 5.1. No Karman vortex shedding was observed for the unsteady simulations. The maximum discrepancy between unsteady to steady simulations is within 2% for this worst case. At lower Re the discrepancy is smaller. Table 5.1 Comparison of the lift and drag coefficient between different solvers, for L D = 10 , α = 40o at Re = 40 Steady Solver Unsteady Solver Difference % Lift Coefficient  0.382  0.374  2.1  Drag Coefficient  1.143  1.119  2.1  5.2.2 Mesh Independence  To ensure that the obtained data are mesh independent, several flows were simulated with different mesh densities. Two typical cases are shown in Figure 5.6 and Figure 5.7. In general, the number of finite volumes required to reach mesh independence increases as Re increases, but for all flows simulated the drag and lift computed with 1.6 million  volumes were within 1% of those computed with a finer mesh. Therefore, the number of volumes used in all other cases was between 1.5 and 2 million.  111  Figure 5.6: Drag coefficient versus number of mesh volumes for L D = 2 , Re = 40 , α = 90 o . Note the highly expanded ordinate axis.  Figure 5.7: Drag coefficient versus number of mesh volumes for L D = 5 , Re = 1 , α = 10o . Note the highly expanded ordinate axis.  5.2.3 Discretization Error  To estimate the discretization error, simulations for several cases are run on different set of grids of increasing fineness. We followed the procedure recommended in [23] for the estimation of the error in the drag and lift coefficients. As we are using integral quantities  112  for error estimation, it is appropriate to choose average global representative cells h1 , h2 , and h3 as ⎡1 h=⎢ ⎣N  1  3 (ΔVi )⎤⎥ ∑ i =1 ⎦ N  Where ΔVi is the volume of the i th cell and N is the total number of cells used for the computations. The procedure requires taking the difference between the extrapolated quantities from course to fine mesh sizes, and then normalizing them by an appropriate grid refinement factor, based on the calculated apparent order of the method. Following this procedure for the meshes used here, the numerical uncertainty in the fine-grid solution was found to be 1% and 1.5% for drag and lift coefficient, respectively. 5.2.4 Code Validation  The drag coefficient for a long cylinder perpendicular to the flow is reported in the literature and tabulated in Table 5.2. A comparison between our simulations and previous simulations and experiments on two-dimensional cylinders is shown in Figure 5.8. To quantify the difference between simulations and both experimental and previous simulation data, the simulation drag coefficients were plotted versus L D (for fixed Re ). An exponential curve fit the simulations well, and the asymptotic limit of the exponential at each value of Re is tabulated in Table 5.2 along with experimental [24] and other numerical simulation results [25]. The rms difference between simulations and experiment is 10% (probably within the large experimental error), while this difference reduces to 1% in comparison with previous numerical simulations.  113  Figure 5.8: Drag coefficient versus Reynolds number ( 0 p Re ≤ 40 ) for a cylinder normal to the flow with L D as a parameter.  Table 5.2 Comparison of drag coefficient between asymptotic value of exponential curve fit to simulations and previous two-dimensional simulations [25] and experiments [24]. The experimental results were read from a logarithmic graph with limited resolution. Simulation Re Ref [25] Difference (%) Ref [24] Difference (%) CD ( L D → ∞ ) 1  12.06  -  -  10.88 ± 0.1  10  5  4.24  -  -  4.44 ± 0.1  -5  10  2.92  2.87  1.7  2.8 ± 0.1  4  15  2.38  2.36  0.8  2.4 ± 0.1  -1  30  1.75  1.73  1.1  1.8 ± 0.1  -3  40  1.55  1.54  0.65  1.6 ± 0.1  -3  At low Reynolds number there is a strong dependence of C D on L D , consistent with r the creeping flow theory of Cox [26], in which the total force F acting upon a circular r cylinder of finite length at rest in a uniform flow of velocity U is given by ⎛L⎞ Fi = μ ⎜ ⎟ AijU j ⎝2⎠  114  where Aij is a tensor given by ( i = 1 refers to the direction along the cylinder axis) 4π ⎧ ⎪ 3 ⎪ ln 2 L D − + ln (2 ) 2 ⎪ 8π ⎪ Aij = ⎨ ⎪ ln 2 L − 1 + ln (2 ) D 2 ⎪ ⎪ 0 ⎪ ⎩  (  (  ,  )  )  i = j =1  ,  i = j = 2,3  ,  i≠ j  For Re = 1 (beyond the range of validity of the creeping flow approximations), Cox’s theory is compared with the simulations in Figure 5.9. The simulations parallel the predictions of Cox’s theory. In summary, the simulations have been validated by the experimental and numerical results on infinite aspect ratio cylinders perpendicular to the flow and also show similar behaviour with L D to the creeping flow theory of Cox.  Figure 5.9: Comparison of drag coefficient between simulations and Cox, Re = 1 .  115  5.3 Results and Discussions  In this section we present the drag and lift coefficients as graphs and best curve fits to the simulation results. The force parallel to the flow (drag) and perpendicular to the flow (lift) were non-dimensionalized as follows CD =  FD  , CL =  FL  1 1 ρU ∞2 LD ρU ∞2 LD 2 2 where L and D are the cylinder length and diameter, respectively, and U ∞ is the free stream velocity.  5.3.1 Flow Visualization  It is instructive to draw some streamlines of the flow for representative inclination angles, aspect ratios and Reynolds numbers. These streamlines help elucidate the observed behaviors in the lift and drag curves. Figure 5.10 shows the streamline patterns for L D = 20 at different inclination angles. The flow Reynolds number is 40.  ( α = 90o )  ( α = 70 o )  116  ( α = 20o ) ( α = 50o ) Figure 5.10: Streamlines around a cylinder with L D = 20 at different inclination angles. The Reynolds number is 40. For a cylinder with L D = 20 normal to the flow ( α = 90 o ), the wake is very small near the tips, and broadens at the cylinder centre. The wake length at the centre of the cylinder is consistent with the experimental measurements of Taneda [27]. At a Reynolds number of 41, for a two-dimensional cylinder, Taneda's experiments show a separation bubble of length 2.2 D . Our simulations, at a Reynolds number of 40 and for a finite cylinder with L D = 20 , show a separation bubble at the middle of the cylinder with a length of 1.8D . The small discrepancy in separation bubble length is perhaps due to an inhibition of the wake owing to the cylinder ends. As the inclination angle decreases, the length of the recirculation zone decreases until, for α ≤ 50o , the wake disappears entirely and the flow is attached everywhere along the cylinder surface. The peak in the lift curve (shown later in Figure 5.14) occurs at about the same angle as the disappearance of the wake. Further decreases in the inclination angle result in a decreased lift coefficient. Similar flow behaviours are observed at smaller aspect ratios (Figure 5.11).  117  ( α = 90 o )  ( α = 40 o )  Figure 5.11: Streamlines around a cylinder with L D = 2 at different inclinations angles. The Reynolds number is 40. At lower Reynolds numbers, e.g. Re = 1 , the flow contains no regions of flow separation (Figure 5.12). The streamlines expand as they approach the cylinder, wrap around the cylinder, and contract downstream, similar to a potential flow.  ( α = 90 o , L D = 20 )  ( α = 40o , L D = 20 )  118  ( α = 40 o , L D = 2 ) ( α = 90 o , L D = 2 ) Figure 5.12: Streamlines around a cylinder with different aspect ratios and inclination angles. Re = 1. At a Reynolds number of 40 for a cylinder normal to the flow with L D = 20 , the pressure at the front stagnation line is nearly constant, with C P = 1.03 (see Figure 5.13a). By way of comparison, 2D simulations by Sen et al. [28] show a stagnation point pressure coefficient of C P = 1.16 at Re = 40 . The difference between these values might be due to the three-dimensionality of the finite cylinder flow, or differences in the numerical dissipation on the front stagnation streamline. The base pressure coefficient at the midspan in our simulations is CP = −0.48 (see Figure 5.13b), which compares favourably with the values reported by Sen et al. [28] ( CP = −0.48 ) and Baranyi and Lewis [25] ( CP = −0.51 ). As tabulated in Table 2, at a Reynolds number of 40 the asymptotic value of C D , obtained by fitting an exponential curve through a plot of C D versus L D , is 1.55. For this Reynolds number, in 2D simulations Baranyi & Lewis [25] report CD = 1.54 and Sen et al. [28] report CD = 1.51 . In both cases the agreement may be considered good, in view of the 1% numerical error and the asymptotic function fitting error.  119  (a) (b) Figure 5.13: Nondimensional pressure distribution along a cylinder with L D = 20 for different inclination angles at Re = 40 . (a) Front pressure (b) Rear pressure. As the inclination angle decreases, the front side pressure decreases, but remains almost constant along the cylinder span. With decreasing inclination angle the rear side pressure rises (Figure 5.13b); more rapidly near s=-0.5 than near s=0.5. These pressure behaviours are consistent with the observed variations in wake shape. 5.3.2 Drag and Lift Coefficient  The left column in Figure 5.14 shows the drag coefficient versus the angle of inclination for different aspect ratios, with Reynolds number as a parameter. To make the graphs more readable, we have only shown curves for Re = 1 and Re = 5 . Graphs at higher Re are similar, but show a reduced C D . For a given L D and Re , the drag coefficient increases monotonically as the inclination angle increases. Owing to symmetry, the C D versus α graph must have zero slope at both α = 0 o and α = 90 o and therefore must have an inflection point somewhere between. In all cases, this inflection point occurs at  α ≈ 45o .  120  Figure 5.14: Forces on a finite cylinder with ( L D = 2, 5, 10, 20) versus the angle of inclination with Re as a parameter. a) Drag coefficient b) Lift coefficient  121  The right side of Figure 5.14 is similar to the left side, but in this case the lift coefficient is plotted against the inclination angle. The lift coefficient increases as the inclination angle increases, reaches a maximum around α ≈ 45o , and then it decreases to 0 for a cylinder normal to the flow. The peak value of the lift coefficient is diminished as the aspect ratio decreases and the Reynolds number increases. It is instructive to consider the contribution of the skin friction and pressure drag to the total drag on the cylinder. Figure 5.15a shows the ratio of pressure drag to viscous drag as a function of inclination angle at Re = 1 . At this low Reynolds number, the viscous drag somewhat exceeds the pressure drag, more so when the cylinder axis becomes more parallel to the flow. As the Reynolds number increases, the curve shapes are relatively unchanged, although the presence of the wake of course leads to higher pressure drag (Figure 5.15b).  (a) (b) Figure 5.15: Ratio of pressure drag to viscous drag as a function of inclination angle at two different aspect ratios at (a) Re = 1 , (b) Re = 40 . The variation of pressure along the cylinder span generates a net torque around the center of the cylinder. The nondimensionalized torque is plotted against inclination angle in Figure 16. Owing to symmetry, a cylinder normal or parallel to the flow experiences no torque about the centroid. For low Reynolds numbers viscous forces are dominant. There is as a result little top/bottom asymmetry in the flow, and therefore the developed torque is insignificant, as seen in Figure 5.16a. Increasing the flow Reynolds number increases the contribution from the pressure forces, which results in a higher developed torque (Figure 5.16b).  122  (a)  (b)  Figure 5.16: Nondimensionalized torque as a function of the inclination angle for different aspect ratios. (a) Re = 1 , (b) Re = 40 .  5.3.3 Normalized Drag Coefficient  In the previous section the drag coefficient data versus the inclination angle were plotted on a graph only for Re = 1 and Re = 5 . To show the rest of the data on a legible graph, we normalized each graph with the y-axis intercept of the graph, i.e. the drag coefficient of the cylinder normal to the flow. The normalized drag coefficients  CD  CD , ⊥  are shown in  Figure 5.17 for different cases. The trends of all graphs are almost the same, showing a monotonic increase of the normalized drag with the inclination angle. The variation in normalized drag is greatest at higher Reynolds numbers. The variation in normalized drag is also greatest at higher aspect ratios. Best curve fits to these graphs are given in section 3.4. The dependence of C D ,⊥ on Re and L D is shown in Fig. 8.  123  Figure 5.17: Non-dimensional drag coefficient ( L D = 2, 5, 10, 20) versus the angle of inclination with Re as a parameter. 5.3.4 Ratio of Lift to Drag Coefficient  Figure 5.18 shows the ratio of C L C D as a function of α , with L D and Re as parameters. Because the lift coefficient curves show a maximum at angles α ≈ 45o , the ratio of lift to drag coefficient is also a maximum around α = 45o . Furthermore, since C L versus α is nearly symmetric with respect to α = 45o , whereas C D increases with increasing α , the C L C D curve is asymmetric with respect to α = 45o , being higher near α = 0 o than near α = 90 o . The asymmetry of the graphs becomes larger as L D rises.  124  Figure 5.18: Lift-to-drag ratio on a finite cylinder with ( L D = 2, 5, 10, 20) versus the angle of inclination with Re as a parameter. 5.3.5 Curve Fitting to Data  Curve fits to the lift and drag coefficients are provided in Table 5.6 and Table 5.7. The coefficients have dependencies on (i) inclination angle (ii) Reynolds number (iii) aspect ratio. The root mean square (RMS) difference between simulations and curve fits are 13% and 12% for drag coefficient and lift coefficient, respectively. For a given aspect ratio, the maximum difference and RMS % difference calculated for all Reynolds numbers and inclination angles are tabulated in Table 3. The largest discrepancy in the lift coefficients happens when Re = 40 for α = 10o or α = 80o . At smaller Re and inclination angles 10o < α < 80o this discrepancy is limited to less than 10%. The largest discrepancy in drag coefficient occurs for Re = 40 and for α < 30 o .  125  Table 5.3 Maximum and RMS difference between simulations and curve fit. drag coefficient lift coefficient L D Maximum Diff. (%) RMS Diff. (%) Maximum Diff. (%) RMS Diff. (%) 2  29  15  34  11  5  21  11  26  10  10  23  12  34  12  20  27  13  49  13  As an independent test of the accuracy of our curve fits, we ran the simulations for the randomly selected case α = 45o , L D = 7 and Re = 7 , which had not been used as data for the curve fitting. The values of C L and C D are given in Table 5.4. There is a good agreement between the curve fit and simulation results, consistent with the results in Table 5.3. Table 5.4 Comparison of C L and C D between simulations and curve fits for a specific simulation case that had not been used to derive the curve fit ( α = 45o , L D = 7 , Re = 7 ) Simulation Curve fit Difference % CL  0.65  0.63  3.1  CD  2.78  2.93  7.3  According to the Independence Principle, FN =  1 ρC DU N2 LD , where U N is the velocity 2  component normal to the cylinder axis, i.e. U N = U ∞ sin α . Table 5.5 lists the ratio of cylinder drag force predicted by simulations to that from the Independence Principle, for cylinders with various values of α , L D , and Re . In this table the simulation drag is computed from the best fit curves. The Independence Principle drag is computed with CD  based on the Reynolds number of the normal component of velocity  ( Re = ρU N D / μ ). The Independence Principle is accurate to within about 5% for  α ≥ 70o , overestimates the drag by 15% or more for α = 45o , and can be more than 50%  126  in error for α ≤ 30o . The error of the Independence Principle is greater at larger Reynolds numbers and at larger values of L D .  Table 5.5 Ratio of cylinder drag force predicted by simulations to that from the Independence Principle, for cylinders with various values of α, L/D, and Re. α = 10o α = 30o α = 45o α = 45o α = 45o α = 70o α = 70o Cylinder Geometry  FSim FIP  Re = 30  Re = 30  Re = 7  Re = 30  Re = 30  Re = 2  Re = 30  L D = 15  L D = 15  L D=7  L D=7  L D = 15  L D=2  L D = 15  0.60  0.61  0.87  0.80  0.74  1.00  0.94  Table 5. 6 Lift coefficient curve fit ⎛ L ⎞ ⎛L ⎞ ⎛L ⎞ C L ⎜ α , , Re ⎟ = A2 ⎜ , Re ⎟ sin(2α ) + A4 ⎜ , Re ⎟ sin(4α ) ⎝ D ⎠ ⎝D ⎠ ⎝D ⎠  L⎞ δ ⎛ A2 ⎜ Re, ⎟ = 1 + δ 2 D ⎠ Re ⎝  A4 = δ 3 ln(Re) + δ 4 A2  ⎛L⎞ ⎟ = λ1 + λ2 exp(−λ3 L D) D ⎝ ⎠ ⎛L⎞ δ 2 ⎜ ⎟ = λ4 + λ5 exp(−λ6 L D) ⎝D⎠  δ1 ⎜  ⎛L⎞ ⎟ = λ7 + λ8 exp(−λ9 L D) ⎝D⎠ ⎛L⎞ δ 4 ⎜ ⎟ = λ10 + λ11 exp(−λ12 L D) ⎝D⎠  δ3⎜  127  λ1 = 1.9276 λ2 = −1.9124 λ3 = 0.4270 λ4 = 0.3870 λ5 = −0.4420 λ6 = 0.5050  λ7 = 0.0330 λ8 = −0.1366 λ9 = 0.3912 λ10 = 0.0515 λ11 = −0.0588 λ12 = 0.0382  Table 5.7 Drag coefficient curve fit  CD ⎛ L ⎞ ⎛L ⎞ ⎛L ⎞ ⎜α , , Re ⎟ = A2 ⎜ , Re ⎟ cos(2α ) + A0 ⎜ , Re ⎟ C D ,⊥ ⎝ D ⎠ ⎝D ⎠ ⎝D ⎠  (  )  (  )  (  )  L⎞ ⎛ A2 ⎜ Re, ⎟ = β 1 ln(Re ) + β 2 D⎠ ⎝  ⎛L⎞ ⎟ = γ 1 + γ 2 exp − γ 3 L D ⎝D⎠ ⎛L⎞ β 2 ⎜ ⎟ = γ 4 + γ 5 exp − γ 6 L D ⎝D⎠  L⎞ ⎛ A0 ⎜ Re, ⎟ = β 3 ln(Re ) + β 4 D⎠ ⎝  ⎛L⎞ ⎟ = γ 7 + γ 8 exp − γ 9 L D ⎝D⎠ ⎛L⎞ β 4 ⎜ ⎟ = γ 10 + γ 11 exp − γ 12 L D ⎝D⎠  β1 ⎜  β3 ⎜  (  )  γ 1 = −0.0414 γ 2 = 0.07010 γ 3 = 0.9342 γ 4 = −0.2377 γ 5 = 0.2685 γ 6 = 0.1697 γ 7 = −0.0470 γ 8 = 0.0790 γ 9 = 0.8194 γ 10 = 0.7609 γ 11 = 0.2577 γ 12 = 0.1583  in which  ⎛L ⎞ C D ,⊥ ⎜ , Re ⎟ = κ 1 Reκ 2 ⎝D ⎠  ⎛ L ⎞ η1 +η2 ⎟= ⎝D⎠ L D η L ⎛ ⎞ κ 2 ⎜ ⎟ = 3 +η4 ⎝D⎠ L D  κ1 ⎜  128  η1 = 18.7565 η 2 = 10.8586 η 3 = −0.3233 η 4 = −0.6024  5.4. Concluding Remarks  In papermaking, fibers, which are roughly finite aspect ratio cylinders, move at arbitrary angles with respect to the surrounding fluid. Under real papermaking conditions the Reynolds number varies up to 60. CFD simulations of cylinders with different aspect ratios, L D = 2, 5, 10, and 20 oriented to the flow at angles between 0 o ≤ α ≤ 90 o and for 1 ≤ Re ≤ 40 were performed. The results were in very good agreement with the limited  experimental data and simulations that exist for these conditions (respectively, rms errors of 3% and 1%). For fixed L D and Re , the drag coefficient increases monotonically with α , whereas the lift coefficient is approximately parabolic, achieving a maximum near  α = 45o . The curves of C L and C D against α are similar as one varies L D and Re . Both C L and C D decrease with increasing Re . C L increases rapidly with L D for L D ≤ 10 , and thereafter is nearly independent of L D . C D decreases rapidly with L D  for L D ≤ 10 , and thereafter is likewise nearly independent of L D . Curve fits of C L (α , L D , Re ) and C D (α , L D , Re ) are provided. As a test of the “Independence  Principle”, various simulations were compared against the results of the principle. The Independence Principle is reasonably accurate for cylinders nearly perpendicular to the flow (error < 5% for α ≥ 70o ), and overestimates the drag force by increasing amounts as α is diminished (error > 15% for α = 45o and error > 50% for α ≤ 30o ). The results presented in this article have filled an important gap in the literature on drag on finite cylinders for industrial flows such as in papermaking.  5.5. Acknowledgments  The authors are grateful for the financial support provided by AstenJohnson Inc. and the Natural Sciences and Engineering Research Council of Canada.  129  5.6. References  [1] M.M. Zdravkovich. Flow Around Circular Cylinders, Vol 2: Applications, Oxford University Press, New York, 2003. [2] S. Kim S and S.J. Karrila. Microhydrodynamics, Principle and Selected Applications, Dover Publications, Inc. New York, 2005. [3] S.C.R. Dennis and J.D.A. Walker. Calculation of the steady flow past a sphere at low and moderate Reynolds numbers, Journal of Fluid Mechanics, 44:771-789, 1971. [4] S. Weinbaum, M.S. Kolansky, M.J. Gluckman and R. Pfeffer. An approximate theory for incompressible viscous flow past two-dimensional bluff bodies in the intermediate Reynolds number regime O(1) < Re < O(100 ) . Journal of Fluid Mechanics, 77:129-152, 1975. [5] M.S. Kolansky, S. Weinbaum and R.Pfeffer. An approximate theory for the streaming motion past axisymmetric bodies at Reynolds numbers of from 1 to approximately 100. Journal of Fluid Mechanics, 82:583-604, 1976. [6] L. Schouveiler, M. Provansal. Periodic wakes of low aspect ratio cylinders with free hemispherical ends. Journal of Fluids and Structures, 15:565-573, 2001. [7] C.H.K. Williamson. Vortex dynamics in the cylinder wake. Annu. Rev. Fluid. Mech. 28:477-539, 1996. [8] F.H. Shair, A.S. Grove, E. Peterson and A. Acrivos. The effect of confining walls on a steady wake behind a circular cylinder, Journal of Fluid Mechanics, 17:546-550, 1965. [9] D. Sumner, J.L. Heseltine. Tip vortex structure for a circular cylinder with a free end, Journal of Wind Engineering and Industrial Aerodynamics, 96:1185-1196, 2008.  [10] A. Slaouti, J.H. Gerrard. An Experimental investigation of end effects on the wake of a circular cylinder towed through water at low Reynolds numbers. Journal of Fluid Mechanics, 112:297-314, 1981.  [11] C. Norberg. An experimental investigation of the flow around a circular cylinder: influence of aspect ratio, Journal of Fluid Mechanics, 258:287-346, 1994. [12] S. Szepessy, P.W. Bearman. Aspect ratio and end plate effects on vortex shedding from a circular cylinder, Journal of Fluid Mechanics, 234:191-217, 1992. [13] C. Norberg. Flow around a circular cylinder: aspect of fluctuating lift, Journal of Fluids and Structures, 15: 459-469, 2001.  130  [14] C. Norberg. Fluctuating lift on a circular cylinder: review and new measurements. Journal of Fluids and Structures. 17:57-96, 2003.  [15] Q. Inoue, A. Sakuragi. Vortex shedding from a circular cylinder of finite length at low Reynolds numbers. Physics of Fluids, 20:033601-033601-12, 2008. [16] L. Schouveiler, M. Provansal. Periodic wakes of low aspect ratio cylinders with free hemispherical ends. Journal of Fluids and Structures,15:565-573, 2001. [17] M. Shirakashi, Y. Ishida, S.Wakiya. Higher Velocity Resonance of Circular Cylinder in Cross Flow. Trans. ASME J. Fluids Eng, 107:392–396, 1985. [18] S.E. Ramberg. The Effect of yaw and finite length upon the vortex wakes of stationary and vibrating circular cylinders. Journal of Fluid Mechanics, 128:81-107, 1983. [19] W.R. Sears. The Boundary layer of yawed cylinder. J. Aero. Sci, 15:49-52, 1948. [20] A. Roshko. On the development of turbulent flow wakes from vortex streets, NACA report, 1191, 1954. [21] C.W. Van Atta. Experiments on vortex shedding form yawed circular cylinders. AIAA J, 6:931-933, 1968.  [22] A.R. Hanson, Vortex shedding from yawed cylinders. AIAA J, 4:738-740, 1966. [23] P.J. Roache. Quantification of Uncertainty in Computational Fluid Dynamics. Annu. Rev. Fluid. Mech. 29:123-160, 1997.  [24] K.O.L.F. Jayaweera and B.J. Mason. The behaviour of freely falling cylinders and cones in a viscous fluid, Journal of Fluid Mechanics, 22:709-720, 1965. [25] L. Baranyi, R.I. Lewis. Comparison of a grid-based CFD method and vortex dynamics predictions of low Reynolds number cylinder flows, The Aeronautical Journal, 110: 63-71, 2006. [26] R.G. Cox. The motion of long slender bodies in a viscous fluid, part 1, general theory, Journal of Fluid Mechanics, 44:791-810, 1970. [27] S. Taneda. Experimental investigation of the wake behind cylinders and plates at low Reynolds numbers, J. Phys. Soc. Japan, 11:302-307, 1956. [28] S. Sen, S. Mittal, G. Biswas. Steady separated flow past a circular cylinder at low Reynolds numbers, Journal of Fluid Mechanics, 620 : 89-119, 2009.  131  Chapter 6  Closing Remarks 6.1 Introduction  Three-dimensional flows through real forming fabrics were dealt with first in this dissertation. The motion of a single flexible fiber in the flow of a circular cylinder was then considered as the interaction between one fiber and one filament of fabrics. The drag and lift coefficients of inclined finite circular cylinders are documented in a range relevant to papermaking conditions. In this chapter, a summary of findings and conclusion of each chapter are separately presented. Some of the frontiers that offer exciting prospects for future research are then stated.  6.2 Summary of Chapter 2  A novel method to build accurate 3D models of forming fabrics was first presented [1]. CFD simulations of pressure-driven flow through the 3D fabric models were then conducted. The striking agreement between the predicted flow rates of the simulations and those from experimental measurements (less then 1%) verified the CAD modeling of forming fabrics and validated the flow simulations through them. For a number of different pressure drops across fabrics, the bulk flow rates upstream of the fabrics were calculated. The pressure drop coefficient was then plotted as a function of Reynolds number, based on the filament diameter and upstream bulk velocity. It was found that: 1) Best curve fits to the computational data were power laws in the form C p = a Re b , with coefficient of determinations close to 1.  2) At the moderate Reynolds number of the studies, 20 ≤ Re ≤ 140 , all calculated powers on the Reynolds number fell in the range − 1 ≤ b ≤ 0 , where b = −1 corresponds to the creeping flow and b = 0 to the fully turbulent flow regime. 132  The flow non-uniformity on the filament scale was considered responsible for small fibers distribution on fabrics. It was found that: 1) On a representative plane, a quarter of the filament diameter in the upstream of the paper-side filaments, the ratio of maximum to minimum velocities was about 3; up to 3 times more small fibers will initially accumulate on one area of the fabric than on other areas in the forming section. The flow non-uniformity was averaged over areas comparable to the cross-sectional area of a fiber. It was found that: 1) The non-uniformity was smoothed out, as expected. 2) A less pronounced weave-based non-uniformity of larger scale emerged which is believed to influence even long fiber motion. 3) The accelerating motion of fibers towards the fabric depended highly on their spatial location. 4) Small fibers oriented in the CMD would be more affected by the flow nonuniformity than fibers in the MD.  6.2.1 Further Comments  1) At high Reynolds numbers, the pressure drop correlates with the dynamic pressure , ρV 2 , and consequently the pressure coefficient with C p ∝ Re 0 . In contrast, at low Reynolds number, the pressure drop correlates with viscous effects, μV L , which leads to C p ∝ Re −1 . Therefore, considering the sum of the viscous and inertial resistance, the pressure drop may be expressed as Δp = aμV + bρV 2 or C p = a Re −1 + b Re 0 . In this dissertation, pressure drop-flow rate relationships were of the form C p = a Re b . This agrees with some published work [2]. A more general form, however, is C p = a Re −1 + b Re c [3].  133  2) The velocity averaging was performed on a plane parallel to the MD-CMD plane. Therefore, the conclusions hold true only for fibers landing parallel to fabrics. If fibers were at an angle to fabrics, a more rigorous 3D averaging should be performed.  6.3 Summary of Chapter 3  We ran simulations of the flow through 3D models of forming fabrics, built based on the procedure described in chapter 2 [1]. A specific fabric and a variant thereof were used to study two quantities on real forming fabrics; the influence of machine-side filaments and jet-to-wire speed ratio on the flow through a forming fabric. It was found that [5]: 1) The machine-side filaments have a small influence on the flow field upstream of the paper-side filaments, and consequently on the initial fiber mass distribution of pulp fibers. 2) The pressure drop of a fabric was found to be approximately the sum of that through the paper and machine side layers. 3) The rush or drag of the impingement jet slightly changed the air permeability results from 58˚ to 90˚ impingement. 4) The rush or drag had a more pronounced effect on the MD shear stress, which would tend to orient fibres in the machine direction.  6.3.1 Further Comments  Some major assumptions behind this work of forming fabric modeling should be iterated: 1) The free surface interaction between the suspension phase and the surrounding air was neglected. 2) Fibers in the suspension were neglected and therefore a pure water flow through forming fabrics was simulated. 3) The inlet walls of the simulation box were tilted to consider the rush and drag of the impingement jet.  134  The zero-shear wall boundary conditions on inlet walls are prone to produce biased results. Therefore, caution should be exercised, given the assumptions that were made to simulate the flow, to prevent any misinterpretation of the results.  6.4 Summary of Chapter 4  We assumed a highly diluted pulp suspension, where the fiber/fiber interactions can be neglected. Therefore, only single fibers were tracked in flow fields. A circular cylinder represented one filament of fabrics. Fibers were then assumed to move only on a plane perpendicular to the cylinder axis. This made a problem of the “single fiber motion in the flow field of a cylinder” [6]. Each fiber was represented by a chain of rigid spheroids connected by ball and socket joints, originally developed by Klingenberg et al [7]. A linear dependence of hydrodynamic drag force and torque on translational and rotational velocities was employed in the simulations. The fiber model was validated against different benchmarks. The inextensibility of the fiber modeling led to some system indeterminacy. This problem was treated by breaking the system down into two subsystems to close the equations of motion. It was found that 1) There were seven dimensionless groups that govern the motion of a flexible fiber in the flow and its interaction with the cylinder. They were Π 1 = ρU ∞ D μ (flowcylinder property), Π 2 = L d and Π 3 = EI μU ∞ L3 (fiber and fiber-flow property), Π 4 = L D (fiber-cylinder geometrical property), Π 5 = f c (fiber-cylinder property), and  Π6 = χ D , Π7 = θ .  2) The hang-up region of fiber was defined based on some combinations of these dimensionless groups. In general, increasing fiber flexibility, fiber length, and coefficient of friction resulted in a higher likelihood of the hung-up region. Conversely, increasing the Reynolds number decreased the likelihood of a fiber becoming hang-up. 3) Finally, the flexible fibers were found to follow the flow streamlines more closely than the rigid ones.  135  6.4.1 Further Comments  The lack of inertia, inextensibility in the fiber model, and low Reynolds number were some major assumptions in the fiber model simulations. To prevent the fiber from passing through the wall, the normal component of velocity was set to zero. Another way of doing this could be introducing a repulsive force coming into action to satisfy the “no-pass-through” condition [8]. Away from fabrics, the relative velocity between fiber and the surrounding water is small. Therefore, the low Re assumption for fibers is valid, whereas for a fiber in contact with a wall, the Reynolds number may reach up to 40. Therefore, using the low Re theory we can acceptably predict the trend of the desired graphs qualitatively but not quantitatively. The velocity-force relations used are valid for a fiber in the free stream velocity; there is no presence of wall effects. Similar relations were also used when fiber was in the vicinity of a wall. This may have led to constant offsets in our results, compared to if true relations were used instead.  6.5 Summary of Chapter 5  CFD simulations of cylinders with different aspect ratios, L D = 2, 5, 10, and 20 oriented to the flow at angles between 0 o ≤ α ≤ 90 o and for 1 ≤ Re ≤ 40 were performed, at conditions relevant to fibers in papermaking conditions [9]. The results were in very good agreement with the limited experimental data and simulations that exist for these conditions. It was found that: 1) For fixed L D and Re , the drag coefficient increased monotonically with α , whereas the lift coefficient was approximately parabolic, achieving a maximum near  α = 45o . The curves of C L and C D against α were similar as one varies L D and Re . Both C L and C D decreased with increasing Re . 2) C L increased rapidly with L D for L D ≤ 10 , and thereafter was nearly independent of L D . C D decreased rapidly with L D for L D ≤ 10 , and thereafter was likewise nearly independent of L D .  136  3) Curve fits of C L (α , L D , Re ) and C D (α , L D , Re ) were provided. 4) As a test of the “Independence Principle”, various simulations were compared against the results of the principle. 5) The Independence Principle was reasonably accurate for cylinders nearly perpendicular to the flow (error < 5% for α ≥ 70o ), and overestimated the drag force by increasing amounts as α was diminished (error > 15% for α = 45o and error > 50% for  α ≤ 30o ). 6) The non-dimensional torque coefficient is zero at inclination angles α = 0 o and  α = 90 o . It should be noted that while α = 0 o is an unstable condition, a cylinder oriented with α = 90 o is stable. Hence, if a cylinder aligned with the incoming flow is disturbed, it rotates to a position normal to the flow. That is why fibers with arbitrary initial orientation in a tank at equilibrium settle normal to the flow.  6.6 Recommendations and Future Work  Some potential future works suggested by this thesis are summarized below.  Chapter 2  1) Simulations of the fiber motion in the 3D flow fields of fabrics to investigate the initial mass distribution of fibers. This may be compared with the hypothesis of the linear proportionality of the initial mass distribution to the velocity fields [1]. 2) An experimental study could be used to validate the computational results. A PIV (Particle Image Velocitmetry) technique could be employed to visualize the flow upstream of the fabric. The global and local pressure measurements would provide compelling evidence for the numerical validity.  Chapter 3  1) This dissertation simulated the rush and drag of jet impingement by tilting the inlet volume of the simulations [5]. A more comprehensive study might take into account the free surface interaction of the pulp suspension and the surrounding air.  137  2) The air permeability results were presented as power laws. A study may correlate the constants of the power laws to the pores structures and the flow properties [3,4].  Chapter 4  1) An experimental study is highly recommended to validate the computational results of the hung-up regions [6]. 2) The fiber model should be extended to an extensible fiber model to remove the system indeterminacy for cases where several segments of a flexible fiber are in contact with a wall. 3) A row of cylinders is an acceptable representation for the paper-side layer of fabrics. Therefore, the fiber motion in the flow through rows of cylinders may give an insight into the influence of both topography and hydrodynamics on wire mark [10,11] . 4) After the extension of the fiber model to 3D, it would be useful to study the motion of a fiber in the flow of a cylinder where the fiber does not lie on a plane perpendicular to the cylinder axis. This adds two more dimensionless parameters (spatial orientations) which complete the hang-up criteria under study [6]. 5) In the attempt to model sheet formation, the fiber/fiber interaction should be incorporated into the fiber model. 6) Replacing the drag coefficients with more generalized formulas would simulate fibers in real papermaking conditions where Re >1.  Chapter 5  1) The drag/lift/torque coefficients were presented for a translational flow without rotation. The next step would be CFD simulations of a rotating cylinder in a quiescent flow to obtain the rotational resistance. 2) The superposition of results, rotational and translational resistance, at Re>1 is questionable. Therefore, CFD simulations of a rotating cylinder in a translational flow need to be carried out to test the range of the superposition validity. 3) CFD simulations of a cylinder bounded by a wall provide the no-slip wall effect on the drag and lift coefficients of a cylinder in the free stream flow.  138  6.7 References  [1] A. Vakil, A. Olyaei, and S.I. Green. Three-Dimensional Geometry and Flow Field Modeling of Forming Fabrics. Nordic Pulp and Paper Research Journal, 24:338-346, 2009. [2] W.M. Lu, K.L. Tung, and K.J. Hwang. Fluid flow through basic weaves of monofilament filter cloth. Textile Research Journal, 66(5), 311-323, 1996. [3] W.L. Ingmanson, S.T. Han, H.D. Wilder and W.T. Myers. Resistance of wire screen to flow of water. Tappi Journal, 44(1):47-54, 1961. [4] J.C. Armour and J.N. Cannon. Fluid flow through woven screens. American Institute of Chemical Engineers, 14(3):415-420, 1968.  [5] Ali Vakil, Arash Olyaei, and Sheldon Green. The Influence of Machine-Side Filaments and Jet-to-Wire Speed Ratio on the Flow Through a Forming Fabric. Submitted for publication, 2009. [6] A. Vakil and S.I. Green. Flexible fiber motion in the flow field of a cylinder. submitted for publication, 2010. [7] R.F. Ross, D. Klingenberg. Dynamic simulation of flexible fibers composed of linked rigid bodies. J. Chem. Phys., 106:2949-2960, 1997. [8] S.B. Lindstrom. Modeling and simulation of paper structure development. Ph.D. Thesis, Mid Sweden University, 2008. [9] A. Vakil and S.I. Green. (2009) Drag and Lift Coefficients of Inclined Finite Circular Cylinders at Moderate Reynolds Numbers. Computers and Fluids, 38:1771-1781. [10] T. Helle. Analysis of wire mark in printing paper. Journal of Pulp and Paper Science. 14(4): J91-J95, 1988.  [11] T. Helle. How forming fabric design affects drainage and release. Pulp and Paper Canada. 79(11):91-98, 1978.  139  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0069055/manifest

Comment

Related Items